Page MenuHomec4science
Files F62240114

scattering1_pod.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 11, 21:09

scattering1_pod.py
View Options

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.linear_problem.tridimensional import Scattering1d as S1D
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 50
size = 2
show_sample = True
show_norm = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
Delta = 0
MN = 5
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
S = [int(np.ceil(R ** (1. / 3.)))] * 3
PODTol = 1e-8
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "distributed":
radial = 0
# radial = "gaussian"
# radial = "thinplate"
# radial = "multiquadric"
rW0 = 5.
radialWeight = [rW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [4., np.pi, 0.]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, -.5], [6., 1.1 * np.pi, .5]]
if size == 2:
mu0 = [4., np.pi, .375]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, 0.], [6., 1.1 * np.pi, .75]]
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[0][2] + bEff*murange[1][2]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
n = 100
solver = S1D(mu0 = mu0, n = n, verbosity = verb)
rescaling = [lambda x: x, lambda x: x, lambda x: x]
rescalingInv = [lambda x: x, lambda x: x, lambda x: x]
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_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':S, 'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import fenics as fen
x = fen.SpatialCoordinate(solver.V.mesh())
warps = [x * mutar[1] / solver._L - x, x * solver._L / mutar[1] - x]
approx.plotApprox(mutar, warps, name = 'u_app',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, warps, name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, warps, name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
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)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(approx.trainedModel.data.mus(0), approx.trainedModel.data.mus(1),
approx.trainedModel.data.mus(2), '.')
ax.set_xlim3d(murangeEff[0][0], murangeEff[1][0])
ax.set_ylim3d(murangeEff[0][1], murangeEff[1][1])
ax.set_zlim3d(murangeEff[0][2], murangeEff[1][2])
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
if algo == "rational" and approx.N > 0:
from plot_zero_set_3 import plotZeroSet3
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], None], [1., 1., 1.],
clip = clip)
plotZeroSet3(murange, murangeEff, approx, mu0, 100, [None, None, None],
[1., 1., 1.], clip = clip, imagTol = 1e-2)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, np.pi, 0.]))

Event Timeline

c4science · Help