Page MenuHomec4science

HelmholtzPadeTaylorApproximant.py
No OneTemporary

File Metadata

Created
Tue, Aug 20, 08:22

HelmholtzPadeTaylorApproximant.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Example homogeneous Dirichlet forcing wave
from __future__ import print_function
import numpy as np
from context import FreeFemHelmholtzEngine as HFEngine
from context import FreeFemHelmholtzScatteringEngine as HFSEngine
from context import FreeFemHelmholtzScatteringAugmentedEngine as HFSAEngine
from context import FreeFemHSEngine as HSEngine
from context import FreeFemHSAugmentedEngine as HSAEngine
from context import ROMApproximantTaylorPade as Pade
PI = np.pi
testNo = 4
if testNo == 1:
params = {'N':3, 'M':3, 'E':3, 'sampleType':'Krylov', 'POD':True}
z0 = 12+.5j
ztar = 11
V = ("int n = 40;\n"
"real nu = 12^.5, theta = pi / 3;\n"
"mesh Th = square(n, n, [pi * x, pi * y]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);")
f = ("16 / pi^4 * exp(-1i * nu*(cos(theta) * x + sin(theta) * y)) * ( 2i *"
" nu *cos(theta) * (2 * x * y^2 - 2 * pi * x * y - pi * y^2 + pi^2 * "
"y) + 2i * nu * sin(theta) * (2 *x^2 * y - pi * x^2 - 2 * pi * x * y "
"+ pi^2 * x) - (2 * y^2 - 2 * pi * y + 2 * x^2 - 2 * pi * x))")
solver = HFEngine(V, z0**.5, forcingTerm = f,
DirichletBoundary = '1,2,3,4')
plotter = HSEngine(solver.V)
approx = Pade(solver, plotter, k0 = z0, w = np.real(z0**.5),
approxParameters = params)
approx.plotApp(ztar, name = 'u_Pade''')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(appErr, solNorm, appErr/solNorm)
print(approx.getPoles(True))
############
elif testNo == 2:
params = {'N':7, 'M':7, 'E':7, 'sampleType':'Arnoldi', 'POD':True}
ztar = 3.95 ** 2
kappa = 4
theta = np.pi * 45 / 180
k1 = kappa * 1 * np.cos(theta)
if 2 * kappa > k1:
k2 = (4. * kappa**2. - k1**2.)**.5
else:
k2 = 1.j * (k1**2. - 4. * kappa**2.)**.5
R = (kappa * 1 * np.sin(theta) - k2) / (kappa * 1 * np.sin(theta) + k2)
T = R + 1
V = ("int n = 50;\n"
"real kappa = {}, theta = {};\n"
"mesh Th = square(n, n, [pi * (x - .5), pi * (y - .5)]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);").format(kappa, theta)
n = "(y >= 0 ? 2. : 1.)"
u0 = ("(y >= 0 ? {0}*exp(1.i*({1}*x+{2}*y)) : exp(1.i*{3}*({4}*x+{5}*y)) "
"+ {6}*exp(1.i*{3}*({4}*x-{5}*y)))").format(T, k1, k2, kappa,
np.cos(theta),
np.sin(theta), R)
solver = HFEngine(V, kappa, refractionIndex = n,
DirichletBoundary = "1,2,3,4", DirichletDatum = u0)
plotter = HSEngine(solver.V)
approx = Pade(solver, plotter, k0 = kappa ** 2,
w = kappa, approxParameters = params,
plotDer = 'ALL')
approx.plotApp(ztar, name = 'u_Pade''')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(appErr, solNorm, appErr/solNorm)
print(approx.getPoles(True))
############
elif testNo == 3:
V = ("int n = 50;\n"
"real nu = 2, theta = pi / 6;\n"
"mesh Th = square(n, n, [pi * x, pi * y]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);")
f = ("4 / pi^2 * exp(-1i * nu * (cos(theta) * x + sin(theta) * y)) * "
"(2i * nu * cos(theta) * (pi - 2 * x) + 2)")
params = {'N':8, 'M':7, 'E':8, 'sampleType':'Krylov', 'POD':True}
kappa = 2.
solver = HFSEngine(V, kappa, forcingTerm = f, DirichletBoundary = "2,4",
RobinBoundary = "1,3")
plotter = HSEngine(solver.V)
approx = Pade(solver, plotter, k0 = kappa,
approxParameters = params,
plotDer = 'ALL')
approx.setupApprox()
ktar = 1.8 - .3j
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
print(appErr, solNorm, appErr/solNorm)
print(approx.getPoles(True))
############
elif testNo == 4:
V = ("int n = 50;\n"
"real nu = 2, theta = pi / 6;\n"
"mesh Th = square(n, n, [pi * x, pi * y]);\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);")
f = ("4 / pi^2 * exp(-1i * nu * (cos(theta) * x + sin(theta) * y)) * "
"(2i * nu * cos(theta) * (pi - 2 * x) + 2)")
params = {'N':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
kappa = 2.
solver = HFSAEngine(V, kappa, forcingTerm = f, DirichletBoundary = "2,4",
RobinBoundary = "1,3", constraintType = "IDENTITY")
plotter = HSAEngine(solver.V, 2)
approx = Pade(solver, plotter, k0 = kappa,
approxParameters = params,
plotDer = 'ALL')
approx.setupApprox()
ktar = 1.8 - .3j
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr = approx.approxError(ktar, kappa)
solNorm = approx.HFNorm(ktar, kappa)
print(appErr, solNorm, np.divide(appErr, solNorm))
print(approx.getPoles(True))

Event Timeline