Page MenuHomec4science
Files F90053442

HelmholtzTaylorPoleIdentification.py
No OneTemporary
Actions

File Metadata

Created
Mon, Oct 28, 21:36

HelmholtzTaylorPoleIdentification.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
import numpy as np
from context import utilities
from context import FenicsHelmholtzEngine as HFEngine
from context import FenicsHSEngine as HSEngine
from context import ROMApproximantTaylorPade as Pade
from context import ROMApproximantTaylorRB as RB
from FEniCS_snippets import SquareHomogeneousBubble
boundary, mesh, forcingTerm = 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(mesh = mesh, wavenumber = z0**.5, forcingTerm = forcingTerm,
FEDegree = 3, DirichletBoundary = boundary,
DirichletDatum = 0)
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')

Event Timeline

c4science · Help