Page MenuHomec4science
Files F60464753

square_pod.py
No OneTemporary
Actions

File Metadata

Created
Tue, Apr 30, 10:44

square_pod.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pole_matching import \
RationalInterpolantPoleMatching as RIPM
from rrompy.reduction_methods.pole_matching import \
ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 7
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
MN = 6
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 2
matchingWeight = 10.
samples = "centered"
samples = "standard"
samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples == "pole matching":
radialM = ""
radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [6.5 ** .5, .9 ** .5]
murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]]
aEff = 1.1
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
if ignore_forcing: solver.nbs = 1
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV"
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RIPM
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RBPM
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'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "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 != "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:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), 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":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 20, [2., 1.], clip = clip,
relative = False, nobeta = True)

Event Timeline

c4science · Help