Page MenuHomec4science
Files F63805542

matrix_pod.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 22, 14:43

matrix_pod.py
View Options

import numpy as np
import scipy.sparse as sp
from matplotlib import pyplot as plt
from rrompy.hfengines.base import MatrixEngineBase as MEB
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as PadeT
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as PadeL
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RBT
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RBL
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
test = 2
method = "PadeTaylor"
method = "PadeLagrange"
method = "RBTaylor"
method = "RBLagrange"
verb = 0
N = 100
solver = MEB(verbosity = verb)
solver.nAs = 2
if test == 1:
solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- sp.eye(N)]
elif test == 2:
solver.setSolver("SOLVE")
fftB = np.fft.fft(np.eye(N)) * N**-.5
solver.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
- np.eye(N)]
np.random.seed(420)
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
mu0 = 10.25
mutar = 12.5
murange = [5.25, 15.25]
if method == "PadeTaylor":
params = {'N':10, 'M':9, 'E':10, 'sampleType':'Arnoldi', 'POD':True}
approx = PadeT(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "PadeLagrange":
params = {'N':10, 'M':9, 'S':11, 'POD':True, 'polybasis':"CHEBYSHEV"}
approx = PadeL(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.sampler = QS(murange, "CHEBYSHEV")
elif method == "RBTaylor":
params = {'R':10, 'E':10, 'sampleType':'Arnoldi', 'POD':True}
approx = RBT(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RBLagrange":
params = {'R':10, 'S':11, 'POD':True}
approx = RBL(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.sampler = QS(murange, "CHEBYSHEV")
approx.setupApprox()
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
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)))
polesTrue = np.arange(1, 1 + N)
polesTrue = polesTrue[polesTrue >= murange[0]]
polesTrue = polesTrue[polesTrue <= murange[1]]
polesApp = approx.getPoles()
mask = (np.real(polesApp) < murange[0]) | (np.real(polesApp) > murange[1])
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.plot(polesTrue, np.zeros_like(polesTrue), 'm.')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()

Event Timeline

c4science · Help