Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/examples/pod/TaylorPoles.py
51cd9a0a331adevelop

Last Change
Raw File

TaylorPoles.py
View Options

import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
from rrompy.utilities.base import squareResonances
verb = 0
k0 = (12+0.j) ** .5
Nmin, Nmax = 2, 10
Nvals = np.arange(Nmin, Nmax + 1, 2)
params = {'N':Nmin, 'M':0, 'Emax':Nmax, 'POD':True, 'sampleType':'Arnoldi'}
#, 'robustTol':1e-14}
#boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x)
# - np.real(z0))
#cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True}
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 25, verbosity = verb)
solver.omega = np.real(k0)
approxP = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb)#,
# equilibration = True, cleanupParameters = cleanupParameters)
approxR = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
rP, rE = [None] * len(Nvals), [None] * len(Nvals)
verbose = 1
for j, N in enumerate(Nvals):
if verbose > 0:
print('N = E = {}'.format(N))
approxP.approxParameters = {'N':N, 'E':N}
approxR.approxParameters = {'R':N, 'E':N}
if verbose > 1:
print(approxP.approxParameters)
print(approxR.approxParameters)
rP[j] = approxP.getPoles()
rE[j] = approxR.getPoles()
if verbose > 2:
print(rP)
print(rE)
from matplotlib import pyplot as plt
plotRows = int(np.ceil(len(Nvals) / 3))
fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows))
for j, N in enumerate(Nvals):
i1, i2 = int(np.floor(j / 3)), j % 3
axes[i1, i2].set_title('N = E = {}'.format(N))
axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb',
label="Pade'", markersize = 8)
axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), 'Pr',
label="RB", markersize = 8)
axes[i1, i2].axhline(linewidth=1, color='k')
xmin, xmax = axes[i1, i2].get_xlim()
height = (xmax - xmin) / 2.
res = np.power(squareResonances(xmin**2., xmax**2., False), .5)
axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4)
axes[i1, i2].plot(np.real(k0), np.imag(k0), 'om', markersize = 5)
axes[i1, i2].plot(np.real(k0) * np.ones(2),
1.5 * height * np.arange(-1, 3, 2), '--m')
axes[i1, i2].grid()
axes[i1, i2].set_xlim(xmin, xmax)
axes[i1, i2].set_ylim(- height, height)
p = axes[i1, i2].legend()
plt.tight_layout()
for j in range((len(Nvals) - 1) % 3 + 1, 3):
axes[plotRows - 1, j].axis('off')
c4science · Help