Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/greedy/squareBubbleHomog.py
f972e1cdd3b2feature_greedy_coarseni…

Last Change
Raw File

squareBubbleHomog.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
from rrompy.utilities.base import squareResonances
from rrompy.utilities.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
errorEstimatorKind = "BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
k0s = np.power(np.linspace(95, 149, 250), .5)
#k0s = np.power(np.linspace(95, 129, 100), .5)
#k0s = np.power(np.linspace(95, 109, 100), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5))
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':10, 'polybasis':polyBasis,
'errorEstimatorKind':errorEstimatorKind, 'interactive':False}
if timed:
verb = 0
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20,
verbosity = verb)
solver.omega = np.real(k0)
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
normApp = np.zeros_like(k0s)
norm = np.zeros_like(k0s)
res = np.zeros_like(k0s)
err = np.zeros_like(k0s)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(polesexact,
2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.')
plt.semilogy(np.real(approx.mus),
4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(polesexact,
2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.semilogy(polesexact,
2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.plot(np.real(polesexact), np.imag(polesexact), 'm.')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
c4science · Help