Page MenuHomec4science
Files F61663507

greedy.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 8, 04:27

greedy.py
View Options

import numpy as np
from diapason_engine import DiapasonEngine, DiapasonEngineDamped
from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RI
from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RB
from rrompy.solver.fenics import L2NormMatrix
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
verb = 5
timed = False
algo = "RI"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
if timed: verb = 0
dampingEta = 0 * 1e4 / 2. / np.pi
k0s = np.linspace(2.5e2, 7.5e3, 100)
#k0s = np.linspace(2.5e3, 1.5e4, 100)
#k0s = np.linspace(5.0e4, 1.0e5, 100)
k0s = np.linspace(2.0e5, 2.5e5, 100)
kl, kr = min(k0s), max(k0s)
theta = 20. * np.pi / 180.
phi = 10. * np.pi / 180.
c = 3.e2
rho = 8e3 * (2. * np.pi) ** 2.
E = 1.93e11
nu = .3
T = 1e6
###
if np.isclose(dampingEta, 0.):
solver = DiapasonEngine(kappa = np.mean(np.power(k0s, 2.)) ** .5, c = c,
rho = rho, E = E, nu = nu, T = T, theta = theta,
phi = phi, meshNo = 1, degree_threshold = 8,
verbosity = 0)
else:
solver = DiapasonEngineDamped(kappa = np.mean(k0s), c = c, rho = rho,
E = E, nu = nu, T = T, theta = theta,
phi = phi, dampingEta = dampingEta,
meshNo = 1, degree_threshold = 8,
verbosity = 0)
params = {'sampler':QS([kl, kr], "UNIFORM"),#, solver.rescalingExp),
'nTestPoints':500, 'greedyTol':1e-2, 'S':5, 'polybasis':polyBasis,
# 'errorEstimatorKind':'DIAGONAL'}
# 'errorEstimatorKind':'INTERPOLATORY'}
'errorEstimatorKind':'AFFINE'}
if algo == "RI":
approx = RI(solver, mu0 = solver.mu0, approxParameters = params,
verbosity = verb)
else:
params.pop("polybasis")
params.pop("errorEstimatorKind")
approx = RB(solver, mu0 = solver.mu0, 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)
polesApp = approx.getPoles()
print("Poles:\n", polesApp)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
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], duality = False))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], duality = False)))
err[j] = approx.normErr(k0s[j]) / norm[j]
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus.data),
1.05*np.max(norm)*np.ones_like(approx.mus.data, 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(np.real(approx.mus.data),
approx.greedyTol*np.ones_like(approx.mus.data, dtype = float),
'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesAppEff = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesAppEff), np.imag(polesAppEff), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()

Event Timeline

c4science · Help