Page MenuHomec4science
Files F61265776

synthetic_pod.py
No OneTemporary
Actions

File Metadata

Created
Sun, May 5, 14:14

synthetic_pod.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
SyntheticBivariateEngine as SBE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 10
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
if samples == "standard":
radial = ""
radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0] * 2
if size == 1: # small
mu0 = [10. ** .5, 15. ** .5]
mutar = [12. ** .5, 14. ** .5]
murange = [[5. ** .5, 10. ** .5], [15 ** .5, 20 ** .5]]
if size == 2: # large
mu0 = [15. ** .5, 17.5 ** .5]
mutar = [18. ** .5, 22. ** .5]
murange = [[5. ** .5, 10. ** .5], [25 ** .5, 25 ** .5]]
if size == 3: # medium
mu0 = [17.5 ** .5, 15 ** .5]
mutar = [20. ** .5, 18. ** .5]
murange = [[10. ** .5, 10. ** .5], [25 ** .5, 20 ** .5]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
kappa = 20. ** .5
theta = - np.pi / 6.
n = 20
L = np.pi
solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
mu0 = mu0, verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':R,
'POD':True, 'PODTolerance':PODTol}
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "standard": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', 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 algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 2.], clip = clip, relative = False)

Event Timeline

c4science · Help