Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FreeFem++/HelmholtzTaylorPoleIdentification.py
6e6c171cf128master

Last Change
Raw File

HelmholtzTaylorPoleIdentification.py
View Options

#!/usr/bin/env python3
import numpy as np
from context import utilities
from context import FreeFemHelmholtzEngine as HFEngine
from context import FreeFemHSEngine as HSEngine
from context import ROMApproximantTaylorPade as Pade
from context import ROMApproximantTaylorRB as RB
from FreeFem_snippets import SquareHomogeneousBubble
bdr, V, f = SquareHomogeneousBubble(kappa = 12 ** .5,
theta = np.pi / 3,
n = 25)
z0 = 12 + 1.j
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 = HFEngine(V, z0 ** .5, forcingTerm = f, DirichletBoundary = bdr)
plotter = HSEngine(solver.V)
approxP = Pade(solver, plotter, k0 = z0, w = np.real(z0**.5),
approxParameters = params)#, equilibration = True,
# cleanupParameters = cleanupParameters)
approxR = RB(solver, plotter, k0 = z0, w = np.real(z0**.5),
approxParameters = params)
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(True)
rE[j] = approxR.getPoles(True)
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]), '*r',
label="RB", markersize = 10)
axes[i1, i2].axhline(linewidth=1, color='k')
xmin, xmax = axes[i1, i2].get_xlim()
res = utilities.squareResonances(xmin, xmax, False)
axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4)
axes[i1, i2].grid()
axes[i1, i2].set_xlim(xmin, xmax)
axes[i1, i2].axis('equal')
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