Page MenuHomec4science
Files F65651996

matrix_pod.py
No OneTemporary
Actions

File Metadata

Created
Wed, Jun 5, 07:29

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.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
test = 2
method = "RationalPade"
#method = "RationalInterpolant"
#method = "RBCentered"
#method = "RBDistributed"
verb = 0
N = 100
solver = MEB(verbosity = verb)
solver.npar = 1
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 == "RationalPade":
params = {'N':10, 'M':9, 'S':11, 'POD':True}
approx = RP(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RationalInterpolant":
params = {'N':10, 'M':9, 'S':11, 'POD':True, 'polybasis':"CHEBYSHEV",
'sampler':QS(murange, "CHEBYSHEV")}
approx = RI(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RBCentered":
params = {'R':10, 'S':11, 'POD':True}
approx = RBC(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RBDistributed":
params = {'R':10, 'S':11, 'POD':True, 'sampler':QS(murange, "CHEBYSHEV")}
approx = RBD(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
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