Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/2d/pod/matrix_passive_pod.py
913de11fda96master

Last Change
Raw File

matrix_passive_pod.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem.tridimensional import \
MatrixDynamicalPassive as MDP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.reduction_methods.pivoting import RationalInterpolantPivoted as RIP
from rrompy.reduction_methods.pivoting import RBDistributedPivoted as RBP
from rrompy.reduction_methods.pivoting import \
RationalInterpolantPivotedPoleMatching as RIPPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 3
show_sample = True
show_norm = True
Delta = 0
MN = 7
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 1
samples = "centered"
samples = "distributed"
samples = "pivoted"
samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples in ["distributed", "pivoted", "pole matching"]:
radial = 0
# radial = "gaussian"
# radial = "thinplate"
# radial = "multiquadric"
rW0 = 10.
radialWeight = [rW0]
if samples in ["pivoted", "pole matching"]:
radialM = 0
radialM = "gaussian"
# radialM = "thinplate"
# radialM = "multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
matchingWeight = 1.
if size == 1:
mu0 = [2.7e2, 20]
mutar = [3e2, 25]
murange = [[20., 10], [5.2e2, 30]]
elif size == 2:
mu0 = [2.7e2, 60]
mutar = [3e2, 75]
murange = [[20., 10], [5.2e2, 110]]
elif size == 3:
mu0 = [2.7e2, 160]
mutar = [3e2, 105]
murange = [[20., 10], [5.2e2, 310]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1]]]
n = 100
b = 10
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
params['radialBasis'] = radial
params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples in ["pivoted", "pole matching"]:
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV"
params['polybasisMarginal'] = "MONOMIAL"
params['radialBasisPivot'] = radial
params['radialBasisMarginal'] = radialM
params['radialBasisWeightsPivot'] = radialWeight
params['radialBasisWeightsMarginal'] = radialWeightM
if samples == "pivoted":
method = RIP
else:
params['matchingWeight'] = matchingWeight
method = RIPPM
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered":
params['S'] = R
method = RBD
elif samples == "pivoted":
params['R'] = SPivot[0]
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = SPivot[1] - 2
params['polybasisMarginal'] = "MONOMIAL"
params['radialBasisMarginal'] = radialM
params['radialBasisWeightsMarginal'] = radialWeightM
method = RBP
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples in ["pivoted", "pole matching"]:
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples not in ["pivoted", "pole matching"]:
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered":
approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
what = "REAL")
approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if hasattr(approx.trainedModel, "getQVal") and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [1., 1])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[1., 1.], relative = False,
nobeta = True)
print(1.j * approx.getPoles([None, 50.]))
c4science · Help