Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/greedy/squareTransmissionNonHomog.py
51cd9a0a331afeature_greedy_coarseni…

Last Change
Raw File

squareTransmissionNonHomog.py
View Options

import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade
verb = 7
homog = True
#homog = False
k0s = np.power(np.linspace(4, 15, 200), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'muBounds':[kl, kr], 'nTrainingPoints':500,
'greedyTol':1e-3, 'nTestPoints':10}
solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4.,
n = 20, verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb, homogeneized = homog)
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
MassInv = approx.massInv
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.normApp(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (solver.norm(MassInv.solve(approx.getRes(k0s[j],
homogeneized = homog)))
/ solver.norm(MassInv.solve(approx.getRHS(k0s[j],
homogeneized = homog))))
err[j] = (approx.normErr(k0s[j], homogeneized = homog)
/ approx.normHF(k0s[j], homogeneized = homog))
resApp = approx.errorEstimator(k0s, homogeneized = homog)
polesApp = approx.getPoles()
polesApp = polesApp[np.abs(np.imag(polesApp)) < 1e-3]
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.semilogy(np.real(polesApp),
2.*np.max(norm)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.semilogy(np.real(polesApp),
2.*np.max(resApp)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.semilogy(np.real(polesApp),
2.*np.max(err)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
c4science · Help