Page MenuHomec4science

HelmholtzVFLagrangeApproximant.py
No OneTemporary

File Metadata

Created
Tue, Aug 20, 00:32

HelmholtzVFLagrangeApproximant.py

#!/usr/bin/env python3
import numpy as np
from context import FenicsHelmholtzEngine as HFEngine
from context import FenicsHelmholtzScatteringEngine as HFSEngine
from context import FenicsHelmholtzScatteringAugmentedEngine as HFSAEngine
from context import FenicsHSEngine as HSEngine
from context import FenicsHSAugmentedEngine as HSAEngine
from context import ROMApproximantLagrangeVF as VF
testNo = 4
if testNo == 1:
params = {'N':4, 'M':3, 'S':5, 'polyBasis':'CHEBYSHEV', 'POD':True}
z0s = [10 + .5j, 14 + .5j]
z0 = np.mean(z0s)
ztar = 11
from FEniCS_snippets import SquareHomogeneousBubble
boundary, mesh, forcingTerm = SquareHomogeneousBubble(kappa = 12 ** .5,
theta = np.pi / 3,
n = 40)
solver = HFEngine(mesh = mesh, wavenumber = z0**.5,
forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = boundary, DirichletDatum = 0)
plotter = HSEngine(solver.V)
approx = VF(solver, plotter, mus = z0s, w = np.real(z0**.5),
approxParameters = params)
approx.plotApp(ztar, name = 'u_RB')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles VF:')
print(approx.getPoles())
############
elif testNo == 2:
params = {'N':9, 'M':8, 'S':10, 'polyBasis':'CHEBYSHEV', 'POD':True}
z0s = np.power([3.85 + .15j, 4.15 + .15j], 2.)
z0 = np.mean(z0s)
ztar = 4 ** 2.
from FEniCS_snippets import SquareTransmissionDirichlet
boundary, mesh, n, u0 = SquareTransmissionDirichlet(nT = 2, nB = 1,
theta = np.pi * 45 / 180,
kappa = 4., n = 50)
solver = HFEngine(mesh = mesh, wavenumber = z0**.5,
refractionIndex = n, FEDegree = 3,
DirichletBoundary = boundary, DirichletDatum = u0)
plotter = HSEngine(solver.V)
approx = VF(solver, plotter, mus = z0s, w = np.real(z0**.5),
approxParameters = params, plotSnap = 'ALL')
approx.plotApp(ztar, name = 'u_VF')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles VF:')
print(approx.getPoles())
############
elif testNo in [3, 4]:
params = {'N':40, 'M':39, 'S':45, 'polyBasis':'CHEBYSHEV', 'POD':True}
k0s = [0, 8]
k0 = np.mean(k0s)
ktar = 4.5
from FEniCS_snippets import SquareScatteringTB
bdrD, bdrN, mesh, forcingTerm = SquareScatteringTB(kappa = 4,
theta = np.pi / 2,
n = 40)
if testNo == 3:
solver = HFSEngine(mesh = mesh, wavenumber = k0, FEDegree = 3,
forcingTerm = forcingTerm, DirichletBoundary = bdrD,
RobinBoundary = bdrN)
plotter = HSEngine(solver.V)
else:
solver = HFSAEngine(mesh = mesh, wavenumber = k0, FEDegree = 3,
forcingTerm = forcingTerm,
DirichletBoundary = bdrD, RobinBoundary = bdrN)
plotter = HSAEngine(solver.V, 2)
approx = VF(solver, plotter, mus = k0s, approxParameters = params)
approx.plotApp(ktar, name = 'u_VF')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles VF:')
print(approx.getPoles())

Event Timeline