Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python/HelmholtzSolver.py
7e2d435a8d1bmaster

Last Change
Raw File

HelmholtzSolver.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
import fenics as fen
import numpy as np
import sympy as sp
from context import FenicsHelmholtzEngine as HF
from context import FenicsHSEngine as HS
testNo = 1
if testNo == 1:
PI = np.pi
def boundary(x, on_boundary):
return on_boundary
nu = 12 ** .5
theta = PI / 3
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 = HF.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu + 1.j, forcingTerm = forcingTerm, FEDegree = 3,\
DirichletBoundary = boundary, DirichletDatum = 0)
uh = solver.solve()
plotter = HS.FenicsHSEngine(solver.V)
print(plotter.norm(uh, np.real(nu)))
plotter.plot(uh)
###########
elif testNo == 2:
def boundary(x, on_boundary):
return on_boundary
PI = np.pi
n1 = 4**.5
n2 = 1**.5
kappa = 3
theta = PI * 70 / 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)
# refraction index
n2Re = fen.Expression('x[1]<0 ? n1r : n2r',
n1r = n1.real, n2r = n2.real, degree=4)
n2Im = 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 = HF.FenicsHelmholtzEngine(mesh = mesh, wavenumber = kappa, refractionIndex = (n2Re, n2Im),\
forcingTerm = 0, FEDegree = 3, DirichletBoundary = boundary,\
DirichletDatum = (uexRe, uexIm))
uh = solver.solve()
plotter = HS.FenicsHSEngine(solver.V)
print(plotter.norm(uh, kappa))
plotter.plot(uh)
###########
elif testNo == 3:
import mshr
from matplotlib import pyplot as plt
PI = np.pi
R = 5
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
A = 10
kappa = 12**.5
theta = - PI * 90 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - A * sp.exp(1.j * phiex)
npoints = 50
scatterer = mshr.Polygon([fen.Point(-1, -.5),
fen.Point(1, -.5),
fen.Point(1, .5),
fen.Point(.8, .5),
fen.Point(.8, -.3),
fen.Point(-.8, -.3),
fen.Point(-.8, .5),
fen.Point(-1, .5),
])
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer, npoints)
plt.jet()
plt.figure()
fen.plot(mesh)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
solver = HF.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0, FEDegree = 3,\
DirichletBoundary = Dboundary, RobinBoundary = 'rest',\
DirichletDatum = DirichletTerm)
baseRe, baseIm = DirichletTerm
baseRe = fen.project(fen.Expression(baseRe, degree = 4), solver.V)
baseIm = fen.project(fen.Expression(baseIm, degree = 4), solver.V)
uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
uh = solver.solve()
plotter = HS.FenicsHSEngine(solver.V)
print(plotter.norm(uh, kappa))
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh)
plotter.plot(uh - uinc)
###########
elif testNo == 4:
import mshr
from matplotlib import pyplot as plt
PI = np.pi
R = 5
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
A = 10
kappa = 12**.5
theta = - PI * 90 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - A * sp.exp(1.j * phiex)
npoints = 40
scatterer = mshr.Polygon([fen.Point(-1, -.5),
fen.Point(1, -.5),
fen.Point(1, .5),
fen.Point(.8, .5),
fen.Point(.8, -.3),
fen.Point(-.8, -.3),
fen.Point(-.8, .5),
fen.Point(-1, .5),
])
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer, npoints)
plt.jet()
plt.figure()
fen.plot(mesh)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
solver = HF.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0, FEDegree = 3,\
DirichletBoundary = Dboundary, RobinBoundary = 'rest',\
DirichletDatum = DirichletTerm, constraintType = 'MASS')
baseRe, baseIm = DirichletTerm
baseRe = fen.project(fen.Expression(baseRe, degree = 4), solver.V)
baseIm = fen.project(fen.Expression(baseIm, degree = 4), solver.V)
uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
uinc = np.concatenate((uinc, kappa * uinc))
uh = solver.solve()
plotter = HS.FenicsHSAugmentedEngine(solver.V, 2)
print(plotter.norm(uh, kappa))
print(plotter.norm(uh - uinc, kappa))
plotter.plot(uh)
plotter.plot(uh - uinc)
c4science · Help