Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/pod/with_error_plot.py
6643b38fd8aemaster

Last Change
Raw File

with_error_plot.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import RationalMovingLeastSquares as RIM
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
verb = 100
sol = "single"
sol = "sweep"
algo = "RI"
algo = "RIM"
#algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
radialBasis = "GAUSSIAN"
#radialBasis = "THINPLATE"
#radialBasis = "MULTIQUADRIC"
#radialBasis = "NEARESTNEIGHBOR"
radialBasisDen = "GAUSSIAN"
radialBasisDen = "THINPLATE"
radialBasisDen = "MULTIQUADRIC"
radialBasisDen = "NEARESTNEIGHBOR"
k0sC = np.power([7 + 0.j, 55 + 0.j], .5)
k0 = np.mean(k0sC ** 2.) ** .5
ktar = 14 ** .5
n = 20
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
params = {'N':1, 'M':1, 'S':30, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0sC, "CHEBYSHEV", 2.)}
if algo == "RI":
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
elif algo == "RIM":
params["radialBasis"] = radialBasis
params["radialDirectionalWeights"] = [.75 * params["S"]]
params["radialBasisDen"] = radialBasisDen
params["nNearestNeighborDen"] = params["N"] + 1
approx = RIM(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
params.pop("N")
params.pop("M")
params.pop("polybasis")
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.setupApprox()
if sol == "single":
approx.plotSamples(what = "REAL")
approx.plotApprox(ktar, what = "REAL", name = "uApp")
approx.plotHF(ktar, what = "REAL", name = "uHF")
approx.plotErr(ktar, what = "REAL", name = "err")
approx.plotRes(ktar, what = "REAL", name = "res")
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
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)))
poles = approx.getPoles()
try:
print('Poles:', poles ** 2.)
except: pass
if sol == "sweep":
z0s = np.real(np.linspace(k0sC[0] ** 2., k0sC[1] ** 2., 100))
k0s = z0s ** .5
zl, zr = min(z0s), max(z0s)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
# normRHS = approx.normRHS(k0s)
norm = approx.normHF(k0s)
normApp = approx.normApprox(k0s)
# res = approx.normRes(k0s) / normRHS
# err = approx.normErr(k0s) / norm
plt.figure()
plt.semilogy(z0s, norm)
plt.semilogy(z0s, normApp, '--')
plt.semilogy(np.real(approx.mus.data) ** 2.,
1.05*np.max(norm)*np.ones_like(approx.mus.data, dtype = float),
'rx')
plt.xlim([zl, zr])
plt.grid()
plt.show()
plt.close()
# plt.figure()
# plt.semilogy(z0s, res)
# plt.xlim([zl, zr])
# plt.grid()
# plt.show()
# plt.close()
#
# plt.figure()
# plt.semilogy(z0s, err)
# plt.xlim([zl, zr])
# plt.grid()
# plt.show()
# plt.close()
#for j, k in enumerate(k0s):
# print(k ** 2., approx.getPoles(mu = k) ** 2., norm[j], normApp[j])
c4science · Help