Page MenuHomec4science
Files F90320367

HelmholtzPadeTaylorApproximant.py
No OneTemporary
Actions

File Metadata

Created
Thu, Oct 31, 11:49

HelmholtzPadeTaylorApproximant.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Example homogeneous Dirichlet forcing wave
from __future__ import print_function
import fenics as fen
import numpy as np
import sympy as sp
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 ROMApproximantTaylorPade as Pade
PI = np.pi
testNo = 4
if testNo == 1:
params = {'N':4, 'M':3, 'E':4, 'sampleType':'Krylov', 'POD':True}
nu = 12**.5
theta = PI / 3
z0 = 12+.5j
ztar = 11
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
nx = ny = 40
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
solver = HFEngine(mesh = mesh, wavenumber = z0**.5, forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = 'all', DirichletDatum = 0)
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':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
ztar = 3.95 ** 2
n1 = 2**.5
n2 = 3**.5
kappa = 4
theta = PI * 75 / 180
d1, d2 = np.cos(theta), np.sin(theta)
K1 = kappa * n1 * d1
if kappa * n2 >= K1:
K2 = ((kappa*n2)**2 - K1**2)**.5
else:
K2 = 1.j * (K1**2 - (kappa*n2)**2)**.5
R = (kappa * n1 * d2 - K2) / (kappa * n1 * d2 + K2)
T = R + 1
x, y = sp.symbols('x[0] x[1]', real=True)
uex1 = T*sp.exp(1.j*(K1*x+K2*y))
uex2 = sp.exp(1.j*kappa*n1*(d1*x+d2*y)) + R*sp.exp(1.j*kappa*n1*(d1*x-d2*y))
# Exact solution
uexRe = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
sp.printing.ccode(sp.re(uex1)), sp.printing.ccode(sp.re(uex2))), degree=4)
uexIm = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
sp.printing.ccode(sp.im(uex1)), sp.printing.ccode(sp.im(uex2))), degree=4)
# wavenumber term
nRe = fen.Expression('x[1]<0 ? n1r : n2r',
n1r = n1.real, n2r = n2.real, degree=4)
nIm = fen.Expression('x[1]<0 ? n1i : n2i',
n1i = n1.imag, n2i = n2.imag, degree=4)
# Create mesh and define function space
nx = ny = 50
mesh = fen.RectangleMesh(fen.Point(-PI/2,-PI/2), fen.Point(PI/2,PI/2), nx, ny)
solver = HFEngine(mesh = mesh, wavenumber = kappa, refractionIndex = (nRe, nIm),
forcingTerm = 0, FEDegree = 3, DirichletBoundary = 'all',
DirichletDatum = (uexRe, uexIm))
plotter = HSEngine(solver.V)
approx = Pade(solver, plotter, k0 = kappa ** 2,
w = kappa, approxParameters = params,
plotDer = True)
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:
def Dboundary(x, on_boundary):
return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
A = 10
kappa = 2
theta = PI * 30 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 4/PI**2 * x * (PI - x)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
nx = ny = 50
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
params = {'N':8, 'M':7, 'E':8, 'sampleType':'Krylov', 'POD':True}
solver = HFSEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = 0)
plotter = HSEngine(solver.V)
approx = Pade(solver, plotter, k0 = kappa,
approxParameters = params,
plotDer = True)
approx.setupApprox()
print(np.roots(approx.Q[::-1]) + kappa)
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:
def Dboundary(x, on_boundary):
return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
A = 10
kappa = 2
theta = PI * 30 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 4/PI**2 * x * (PI - x)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
nx = ny = 20
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
params = {'N':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
solver = HFSAEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = 0,
constraintType = 'IDENTITY')
plotter = HSAEngine(solver.V, 2)
approx = Pade(solver, plotter, k0 = kappa,
approxParameters = params,
plotDer = True)
approx.setupApprox()
print(np.roots(approx.Q[::-1]) + kappa)
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, kappa), approx.HFNorm(ktar, kappa)
print(appErr, solNorm, np.divide(appErr, solNorm))
print(approx.getPoles(True))

Event Timeline

c4science · Help