Page MenuHomec4science
Files F63068273

FenicsHelmholtzEngine.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 17, 13:58

FenicsHelmholtzEngine.py
View Options

#!/usr/bin/python
from copy import copy
import warnings
import numpy as np
import sympy as sp
import sympy.printing as sprint
import scipy.sparse as scsp
import fenics as fen
from RROMPyTypes import DictAny, Np1D, Np2D, GenExpr, FenExpr, BfSExpr, Tuple, List
fenZERO = fen.Constant(0)
fenZEROC = [fenZERO, fenZERO]
class DomainError(Exception):
"""
Domain error exception.
Args:
value: Human readable string describing the exception.
Attributes:
value: Human readable string describing the exception.
"""
def __init__(self, value:str):
self.value = value
def __str__(self):
return self.value
def CustomExpressionParser(value:GenExpr, degree:int)\
-> Tuple[FenExpr, FenExpr]:
"""
Numerical and Fenics expressions parser.
Args:
value: Expression to be parsed. Accepts 2-tuples composed of real and
imaginary parts. Available elementary formats are numbers, strings
and Fenics Expression's.
degree: Degree of Fenics FE interpolant of expression.
Returns:
Fenics FE interpolant of expression.
"""
try:
if isinstance(value, (list, tuple)):
if len(value) == 1:
return CustomExpressionParser(value[0])
elif len(value) != 2:
raise Exception("Parsing error")
if isinstance(value[0], (int, float, complex)):
if any([isinstance(y, complex) for y in value]):
raise Exception("Parsing error")
valueRe = fen.Constant(value[0])
valueIm = fen.Constant(value[1])
elif isinstance(value[0], str):
x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]",
real=True)
xyDict = {"x": x, "y": y, "z": z}
valueReParsed = value[0].replace("x[0]", "x")\
.replace("x[1]", "y")\
.replace("x[2]", "z")
valueImParsed = value[1].replace("x[0]", "x")\
.replace("x[1]", "y")\
.replace("x[2]", "z")
valueReSym = sp.sympify(valueReParsed, locals=xyDict)\
.subs([(x, x0), (y, x1), (z, x2)])
valueImSym = sp.sympify(valueImParsed, locals=xyDict)\
.subs([(x, x0), (y, x1), (z, x2)])
valueReStr = sprint.ccode(valueReSym).rpartition("\n")[-1]
valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1]
valueRe = fen.Expression(valueReStr, degree = degree)
valueIm = fen.Expression(valueImStr, degree = degree)
else:
valueRe = value[0]
valueIm = value[1]
else:
if isinstance(value, (int, float, complex)):
valueRe = fen.Constant(np.real(value))
valueIm = fen.Constant(np.imag(value))
elif isinstance(value, str):
x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]",
real=True)
xyDict = {"x": x, "y": y, "z": z}
valueParsed = value.replace("x[0]", "x").replace("x[1]", "y")\
.replace("x[2]", "z")
valueSym = sp.sympify(valueParsed, locals=xyDict)\
.subs([(x, x0), (y, x1), (z, x2)])
valueReStr = sprint.ccode(sp.re(valueSym)).rpartition("\n")[-1]
valueImStr = sprint.ccode(sp.im(valueSym)).rpartition("\n")[-1]
valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1]
valueRe = fen.Expression(valueReStr, degree = degree)
valueIm = fen.Expression(valueImStr, degree = degree)
else:
valueRe = value
valueIm = fenZERO
except:
raise Exception("Parsing error")
return valueRe, valueIm
class FenicsHelmholtzEngine:
"""
Fenics-based solver for Helmholtz problems.
- \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Args:
mesh: Domain of Helmholtz problem.
FEDegree: FE degree.
wavenumber2: Value of k^2.
diffusivity(optional): Value of a. Defaults to identically 1.
refractionIndex(optional): Value of n. Defaults to identically 1.
forcingTerm(optional): Value of f. Defaults to identically 0.
DirichletBoundary(optional): Function flagging Dirichlet boundary as
True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
accepted. Defaults to False everywhere.
NeumannBoundary(optional): Function flagging Neumann boundary as True,
in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
RobinBoundary(optional): Function flagging Robin boundary as True, in
Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
DirichletDatum(optional): Value of u0. Defaults to identically 0.
NeumannDatum(optional): Value of g1. Defaults to identically 0.
RobinDatum(optional): Value of h. Defaults to identically 0.
RobinDatum2(optional): Value of g2. Defaults to identically 0.
Attributes:
FEDegree: FE degree.
wavenumber2: Copy of processed wavenumber2 parameter.
diffusivity: Copy of processed diffusivity parameter.
refractionIndex: Copy of processed refractionIndex parameter.
forcingTerm: Copy of processed forcingTerm parameter.
DirichletBoundary: Copy of processed DirichletBoundary parameter.
NeumannBoundary: Copy of processed NeumannBoundary parameter.
RobinBoundary: Copy of processed RobinBoundary parameter.
DirichletDatum: Copy of processed DirichletDatum parameter.
NeumannDatum: Copy of processed NeumannDatum parameter.
RobinDatum: Copy of processed RobinDatum parameter.
RobinDatum2: Copy of processed RobinDatum2 parameter.
LU: Whether to use LU solver for computation of derivatives.
V: Real FE space.
u: Helmholtz problem solution as numpy complex vector.
v1: Generic trial functions for variational form evaluation.
v2: Generic test functions for variational form evaluation.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
nRe,nIm: Real and imaginary parts of n.
n2Re,n2Im: Real and imaginary parts of n^2.
fRe,fIm: Real and imaginary parts of f.
u0Re,u0Im: Real and imaginary parts of u0.
g1Re,g1Im: Real and imaginary parts of g1.
hRe,hIm: Real and imaginary parts of h.
g2Re,g2Im: Real and imaginary parts of g2.
DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
parts of Dirichlet data.
A*: Scipy sparse array representation (in CSC format) of A*.
b*: Numpy array representation of b*.
"""
def __init__(self, mesh:fen.mesh,
FEDegree:int,
wavenumber:complex,
diffusivity : GenExpr = 1,
refractionIndex : GenExpr = 1,
forcingTerm : GenExpr = fenZEROC,
DirichletBoundary : BfSExpr = None,
NeumannBoundary : BfSExpr = None,
RobinBoundary : BfSExpr = None,
DirichletDatum : GenExpr = fenZEROC,
NeumannDatum : GenExpr = fenZEROC,
RobinDatum : GenExpr = fenZEROC,
RobinDatum2 : GenExpr = fenZEROC):
self.FEDegree = FEDegree
# process boundaries
boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary]
for i in range(3):
if isinstance(boundariesList[i], str):
boundariesList[i] = boundariesList[i].upper()
if boundariesList[i] == "NONE": boundariesList[i] = None
DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList
if boundariesList.count(None) == 3:
raise DomainError("At least one boundary must be prescribed.")
if boundariesList.count("ALL") + boundariesList.count("REST") >= 2:
raise DomainError("Only one keyword 'ALL'/'REST' can be used.")
def DirichletB(x, on_boundary):
if DirichletBoundary is None:
return False
elif DirichletBoundary == "ALL":
return on_boundary
elif DirichletBoundary == "REST":
return (on_boundary
and not self.NeumannBoundary(x, on_boundary)
and not self.RobinBoundary(x, on_boundary))
else:
return DirichletBoundary(x, on_boundary)
def NeumannB(x, on_boundary):
if NeumannBoundary is None:
return False
elif NeumannBoundary == "ALL":
return on_boundary
elif NeumannBoundary == "REST":
return (on_boundary
and not self.DirichletBoundary(x, on_boundary)
and not self.RobinBoundary(x, on_boundary))
else:
return NeumannBoundary(x, on_boundary)
def RobinB(x, on_boundary):
if RobinBoundary is None:
return False
elif RobinBoundary == "ALL":
return on_boundary
elif RobinBoundary == "REST":
return (on_boundary
and not self.DirichletBoundary(x, on_boundary)
and not self.NeumannBoundary(x, on_boundary))
else:
return RobinBoundary(x, on_boundary)
self.DirichletBoundary = DirichletB
self.NeumannBoundary = NeumannB
self.RobinBoundary = RobinB
# create boundary measure
class NBoundary(fen.SubDomain):
def inside(self, x, on_boundary):
return NeumannB(x, on_boundary)
class RBoundary(fen.SubDomain):
def inside(self, x, on_boundary):
return RobinB(x, on_boundary)
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NBoundary().mark(boundary_markers, 0)
RBoundary().mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain=mesh,
subdomain_data=boundary_markers)
self.V = fen.FunctionSpace(mesh, "P", FEDegree)
self.v1 = fen.TrialFunction(self.V)
self.v2 = fen.TestFunction(self.V)
self.wavenumber2 = wavenumber ** 2
self.diffusivity = diffusivity
self.refractionIndex = refractionIndex
self.forcingTerm = forcingTerm
self.DirichletDatum = DirichletDatum
self.NeumannDatum = NeumannDatum
self.RobinDatum = RobinDatum
self.RobinDatum2 = RobinDatum2
functional = lambda self, u: 0.
def Vdim(self) -> int:
"""Dimension of finite element space."""
return self.V.dim()
def energyNormMatrix(self, w:float) -> Np2D:
"""
Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
Args:
w: Weight.
Returns:
Sparse matrix in CSR format.
"""
normMatFen = fen.assemble((
fen.dot(fen.grad(self.v1), fen.grad(self.v2))
+ w * fen.dot(self.v1, self.v2)
) * fen.dx)
normMat = fen.as_backend_type(normMatFen).mat()
return scsp.csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
def problemData(self) -> DictAny:
"""List of HF problem data."""
dataDict = {}
dataDict["forcingTerm"] = self.forcingTerm
dataDict["DirichletDatum"] = self.DirichletDatum
dataDict["NeumannDatum"] = self.NeumannDatum
dataDict["RobinDatum2"] = self.RobinDatum2
return dataDict
def setProblemData(self, dataDict:DictAny, k2:complex):
"""
Set problem data.
Args:
dataDict: Dictionary of problem data.
k2: Parameter value.
"""
self.wavenumber2 = k2
self.forcingTerm = dataDict["forcingTerm"]
self.DirichletDatum = dataDict["DirichletDatum"]
self.NeumannDatum = dataDict["NeumannDatum"]
self.RobinDatum2 = dataDict["RobinDatum2"]
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
if hasattr(self, "b0"): del self.b0
def liftDirichletData(self) -> Np1D:
"""Lift Dirichlet datum."""
solLRe = fen.interpolate(self.u0Re, self.V)
solLIm = fen.interpolate(self.u0Im, self.V)
return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector())
@property
def diffusivity(self):
"""Value of a. Its assignment changes aRe and aIm."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
if hasattr(self, "A0"): del self.A0
self.aRe, self.aIm = CustomExpressionParser(diffusivity,
degree = self.FEDegree)
self._diffusivity = [self.aRe, self.aIm]
@property
def refractionIndex(self):
"""Value of n. Its assignment changes nRe, nIm, n2Re and n2Im."""
return self._refractionIndex
@refractionIndex.setter
def refractionIndex(self, refractionIndex):
if hasattr(self, "A1"): del self.A1
self.nRe, self.nIm = CustomExpressionParser(refractionIndex,
degree = self.FEDegree)
self.n2Re = self.nRe*self.nRe - self.nIm*self.nIm
self.n2Im = 2 * self.nRe * self.nIm
self._refractionIndex = [self.nRe, self.nIm]
@property
def forcingTerm(self):
"""Value of f. Its assignment changes fRe and fIm."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
if hasattr(self, "b0"): del self.b0
self.fRe, self.fIm = CustomExpressionParser(forcingTerm,
degree = self.FEDegree)
self._forcingTerm = [self.fRe, self.fIm]
@property
def DirichletDatum(self):
"""
Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
DirichletBCIm.
"""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
if hasattr(self, "b0"): del self.b0
self.u0Re, self.u0Im = CustomExpressionParser(DirichletDatum,
degree = self.FEDegree)
self.DirichletBCRe = fen.DirichletBC(self.V, self.u0Re,
self.DirichletBoundary)
self.DirichletBCIm = fen.DirichletBC(self.V, self.u0Im,
self.DirichletBoundary)
self.DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
self._DirichletDatum = [self.u0Re, self.u0Im]
@property
def NeumannDatum(self):
"""Value of g1. Its assignment changes g1Re and g1Im."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
if hasattr(self, "b0"): del self.b0
self.g1Re, self.g1Im = CustomExpressionParser(NeumannDatum,
degree = self.FEDegree)
self._NeumannDatum = [self.g1Re, self.g1Im]
@property
def RobinDatum(self):
"""Value of h. Its assignment changes hRe and hIm."""
return self._RobinDatum
@RobinDatum.setter
def RobinDatum(self, RobinDatum):
if hasattr(self, "A0"): del self.A0
self.hRe, self.hIm = CustomExpressionParser(RobinDatum,
degree = self.FEDegree)
self._RobinDatum = [self.hRe, self.hIm]
@property
def RobinDatum2(self):
"""Value of g2. Its assignment changes g2Re and g2Im."""
return self._RobinDatum2
@RobinDatum2.setter
def RobinDatum2(self, RobinDatum2):
if hasattr(self, "b0"): del self.b0
self.g2Re, self.g2Im = CustomExpressionParser(RobinDatum2,
degree = self.FEDegree)
self._RobinDatum2 = [self.g2Re, self.g2Im]
def assembleA(self):
"""Assemble matrix blocks of linear system."""
if not hasattr(self, "A0"):
a0Re = (self.aRe * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
+ self.hRe * fen.dot(self.v1, self.v2) * self.ds(1))
a0Im = (self.aIm * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
+ self.hIm * fen.dot(self.v1, self.v2) * self.ds(1))
A0Re = fen.assemble(a0Re)
A0Im = fen.assemble(a0Im)
self.DirichletBC0.apply(A0Re)
self.DirichletBC0.zero(A0Im)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
A0Imr, A0Imc, A0Imv = fen.as_backend_type(
A0Im).mat().getValuesCSR()
self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ReMat.size))
if not hasattr(self, "A1"):
a1Re = - self.n2Re * fen.dot(self.v1, self.v2) * fen.dx
a1Im = - self.n2Im * fen.dot(self.v1, self.v2) * fen.dx
A1Re = fen.assemble(a1Re)
A1Im = fen.assemble(a1Im)
self.DirichletBC0.zero(A1Re)
self.DirichletBC0.zero(A1Im)
A1Rer, A1Rec, A1Rev = fen.as_backend_type(A1Re).mat()\
.getValuesCSR()
A1Imr, A1Imc, A1Imv = fen.as_backend_type(A1Im).mat()\
.getValuesCSR()
self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = self.A0.shape)
+ 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
shape = self.A0.shape))
def assembleb(self):
"""Assemble RHS blocks of linear system."""
if not hasattr(self, "b0"):
L0Re = (fen.dot(self.fRe, self.v2) * fen.dx
+ fen.dot(self.g1Re, self.v2) * self.ds(0)
+ fen.dot(self.g2Re, self.v2) * self.ds(1))
L0Im = (fen.dot(self.fIm, self.v2) * fen.dx
+ fen.dot(self.g1Im, self.v2) * self.ds(0)
+ fen.dot(self.g2Im, self.v2) * self.ds(1))
b0Re = fen.assemble(L0Re)
b0Im = fen.assemble(L0Im)
self.DirichletBCRe.apply(b0Re)
self.DirichletBCIm.apply(b0Im)
self.b0 = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
parDegree = 1
def As(self) -> List[Np2D]:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1]
def bs(self) -> List[Np1D]:
"""RHS blocks of linear system."""
self.assembleb()
return [self.b0]
def A(self) -> Np2D:
"""Assemble matrix of linear system."""
self.assembleA()
return self.A0 + self.wavenumber2 * self.A1
def b(self) -> Np1D:
"""Assemble RHS of linear system."""
self.assembleb()
return self.b0
def solve(self) -> Np1D:
"""
Find solution of linear system.
Returns:
Solution vector.
"""
return scsp.linalg.spsolve(self.A(), self.b())
def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]:
"""
Get blocks of linear system.
Returns:
Blocks of LS (\sum_{j=0}^J k^j A_j = \sum_{j=0}^K k^j b_j).
"""
return self.As(), self.bs()
class FenicsHelmholtzScatteringEngine(FenicsHelmholtzEngine):
"""
Fenics-based solver for Helmholtz scattering problems with Robin ABC.
- \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu - i k u = g2 on \Gamma_R
Args:
mesh: Domain of Helmholtz problem.
FEDegree: FE degree.
wavenumber: Value of k.
diffusivity(optional): Value of a. Defaults to identically 1.
refractionIndex(optional): Value of n. Defaults to identically 1.
forcingTerm(optional): Value of f. Defaults to identically 0.
DirichletBoundary(optional): Function flagging Dirichlet boundary as
True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
accepted. Defaults to False everywhere.
NeumannBoundary(optional): Function flagging Neumann boundary as True,
in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
RobinBoundary(optional): Function flagging Robin boundary as True, in
Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
DirichletDatum(optional): Value of u0. Defaults to identically 0.
NeumannDatum(optional): Value of g1. Defaults to identically 0.
RobinDatum2(optional): Value of g2. Defaults to identically 0.
Attributes:
FEDegree: FE degree.
wavenumber: Copy of processed wavenumber parameter.
diffusivity: Copy of processed diffusivity parameter.
refractionIndex: Copy of processed refractionIndex parameter.
forcingTerm: Copy of processed forcingTerm parameter.
DirichletBoundary: Copy of processed DirichletBoundary parameter.
NeumannBoundary: Copy of processed NeumannBoundary parameter.
RobinBoundary: Copy of processed RobinBoundary parameter.
DirichletDatum: Copy of processed DirichletDatum parameter.
NeumannDatum: Copy of processed NeumannDatum parameter.
RobinDatum: Value of h, i.e. (- i * k).
RobinDatum2: Copy of processed RobinDatum2 parameter.
V: Real FE space.
u: Helmholtz problem solution as numpy complex vector.
v1: Generic trial functions for variational form evaluation.
v2: Generic test functions for variational form evaluation.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
nRe,nIm: Real and imaginary parts of n.
n2Re,n2Im: Real and imaginary parts of n^2.
fRe,fIm: Real and imaginary parts of f.
u0Re,u0Im: Real and imaginary parts of u0.
g1Re,g1Im: Real and imaginary parts of g1.
hRe,hIm: Real and imaginary parts of h.
h: Value of h.
g2Re,g2Im: Real and imaginary parts of g2.
DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
parts of Dirichlet data.
A*: Scipy sparse array representation (in CSC format) of A*.
b*: Numpy array representation of b*.
"""
def __init__(self, mesh:fen.mesh,
FEDegree:int,
wavenumber:complex,
diffusivity : GenExpr = 1,
refractionIndex : GenExpr = 1,
forcingTerm : GenExpr = fenZEROC,
DirichletBoundary : BfSExpr = None,
NeumannBoundary : BfSExpr = None,
RobinBoundary : BfSExpr = None,
DirichletDatum : GenExpr = fenZEROC,
NeumannDatum : GenExpr = fenZEROC,
RobinDatum2 : GenExpr = fenZEROC):
self.wavenumber = wavenumber
FenicsHelmholtzEngine.__init__(self, mesh = mesh,
FEDegree = FEDegree,
wavenumber = wavenumber,
diffusivity = diffusivity,
refractionIndex = refractionIndex,
forcingTerm = forcingTerm,
DirichletBoundary = DirichletBoundary,
NeumannBoundary = NeumannBoundary,
RobinBoundary = RobinBoundary,
DirichletDatum = DirichletDatum,
NeumannDatum = NeumannDatum,
RobinDatum = -1.j * wavenumber,
RobinDatum2 = RobinDatum2)
def energyNormMatrix(self, w:float) -> Np2D:
"""
Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
Args:
w: Weight.
Returns:
Sparse matrix in CSR format.
"""
return FenicsHelmholtzEngine.energyNormMatrix(self, w ** 2.)
def setProblemData(self, dataDict:DictAny, k:complex):
"""
Set problem data.
Args:
dataDict: Dictionary of problem data.
k: Parameter value.
"""
self.wavenumber = k
self.forcingTerm = dataDict["forcingTerm"]
self.DirichletDatum = dataDict["DirichletDatum"]
self.NeumannDatum = dataDict["NeumannDatum"]
self.RobinDatum2 = dataDict["RobinDatum2"]
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
if hasattr(self, "A2"): del self.A2
if hasattr(self, "b0"): del self.b0
@property
def wavenumber2(self):
"""Value of k^2."""
return self._wavenumber2
@wavenumber2.setter
def wavenumber2(self, wavenumber2):
self._wavenumber2 = wavenumber2
self._wavenumber = self.wavenumber2 ** .5
if (not hasattr(self, "h")
or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)):
self.RobinDatum = -1.j * self.wavenumber
@property
def wavenumber(self):
"""Value of k."""
return self._wavenumber
@wavenumber.setter
def wavenumber(self, wavenumber):
self.wavenumber2 = wavenumber ** 2.
@property
def RobinDatum(self):
"""
Value of h. Its assignment changes hRe and hIm (and maybe kRe, kIm, k
and z).
"""
return super().RobinDatum
@RobinDatum.setter
def RobinDatum(self, RobinDatum):
self.h = RobinDatum
self.hRe, self.hIm = np.real(self.h), np.imag(self.h)
self._RobinDatum = [self.hRe, self.hIm]
if (not hasattr(self, "wavenumber")
or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)):
self.wavenumber = 1.j * self.h
def assembleA(self):
"""Assemble matrix blocks of linear system."""
if not hasattr(self, "A0"):
a0Re = self.aRe * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
a0Im = self.aIm * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
A0Re = fen.assemble(a0Re)
A0Im = fen.assemble(a0Im)
self.DirichletBC0.apply(A0Re)
self.DirichletBC0.zero(A0Im)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
A0Imr, A0Imc, A0Imv = fen.as_backend_type(
A0Im).mat().getValuesCSR()
self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ReMat.size))
if hasattr(self, "A1"):
aux = copy(self.A1)
else:
aux = None
if not hasattr(self, "A2"):
if aux is not None:
del self.A1
FenicsHelmholtzEngine.assembleA(self)
self.A2 = copy(self.A1)
if aux is None:
a1 = fen.dot(self.v1, self.v2) * self.ds(1)
A1 = fen.assemble(a1)
self.DirichletBC0.zero(A1)
A1r, A1c, A1v = fen.as_backend_type(A1).mat().getValuesCSR()
self.A1 = - 1.j * scsp.csr_matrix((A1v, A1c, A1r),
shape = self.A0.shape)
else:
self.A1 = aux
parDegree = 2
def As(self) -> List[Np2D]:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1, self.A2]
def A(self) -> Np2D:
"""Assemble matrix of linear system."""
self.assembleA()
return self.A0 + self.wavenumber * self.A1 + self.wavenumber2 * self.A2
class FenicsHelmholtzScatteringAugmentedEngine(FenicsHelmholtzScatteringEngine):
"""
Fenics-based solver for Helmholtz scattering problems with Robin ABC.
- \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu - i k u = g2 on \Gamma_R
Linear dependence on k is achieved by introducing an additional variable.
Args:
mesh: Domain of Helmholtz problem.
FEDegree: FE degree.
wavenumber: Value of k.
diffusivity(optional): Value of a. Defaults to identically 1.
refractionIndex(optional): Value of n. Defaults to identically 1.
forcingTerm(optional): Value of f. Defaults to identically 0.
DirichletBoundary(optional): Function flagging Dirichlet boundary as
True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
accepted. Defaults to False everywhere.
NeumannBoundary(optional): Function flagging Neumann boundary as True,
in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
RobinBoundary(optional): Function flagging Robin boundary as True, in
Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
Defaults to False everywhere.
DirichletDatum(optional): Value of u0. Defaults to identically 0.
NeumannDatum(optional): Value of g1. Defaults to identically 0.
RobinDatum2(optional): Value of g2. Defaults to identically 0.
constraintType(optional): Type of augmentation. Keywords 'IDENTITY' and
'MASS' are accepted. Defaults to 'IDENTITY'.
Attributes:
FEDegree: FE degree.
wavenumber: Copy of processed wavenumber parameter.
diffusivity: Copy of processed diffusivity parameter.
refractionIndex: Copy of processed refractionIndex parameter.
forcingTerm: Copy of processed forcingTerm parameter.
DirichletBoundary: Copy of processed DirichletBoundary parameter.
NeumannBoundary: Copy of processed NeumannBoundary parameter.
RobinBoundary: Copy of processed RobinBoundary parameter.
DirichletDatum: Copy of processed DirichletDatum parameter.
NeumannDatum: Copy of processed NeumannDatum parameter.
RobinDatum: Value of h, i.e. (- i * k).
RobinDatum2: Copy of processed RobinDatum2 parameter.
constraintType: Type of augmentation.
V: Real FE space.
u: Helmholtz problem solution as numpy complex vector.
v1: Generic trial functions for variational form evaluation.
v2: Generic test functions for variational form evaluation.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
nRe,nIm: Real and imaginary parts of n.
n2Re,n2Im: Real and imaginary parts of n^2.
fRe,fIm: Real and imaginary parts of f.
u0Re,u0Im: Real and imaginary parts of u0.
g1Re,g1Im: Real and imaginary parts of g1.
hRe,hIm: Real and imaginary parts of h.
g2Re,g2Im: Real and imaginary parts of g2.
DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
parts of Dirichlet data.
A*: Scipy sparse array representation (in CSC format) of A*.
b*: Numpy array representation of b*.
"""
def __init__(self, mesh:fen.mesh,
FEDegree:int,
wavenumber:complex,
diffusivity : GenExpr = 1,
refractionIndex : GenExpr = 1,
forcingTerm : GenExpr = fenZEROC,
DirichletBoundary : BfSExpr = None,
NeumannBoundary : BfSExpr = None,
RobinBoundary : BfSExpr = None,
DirichletDatum : GenExpr = fenZEROC,
NeumannDatum : GenExpr = fenZEROC,
RobinDatum2 : GenExpr = fenZEROC,
constraintType : str = "IDENTITY"):
self.constraintType = constraintType
FenicsHelmholtzScatteringEngine.__init__(self, mesh = mesh,
FEDegree = FEDegree,
wavenumber = wavenumber,
diffusivity = diffusivity,
refractionIndex = refractionIndex,
forcingTerm = forcingTerm,
DirichletBoundary = DirichletBoundary,
NeumannBoundary = NeumannBoundary,
RobinBoundary = RobinBoundary,
DirichletDatum = DirichletDatum,
NeumannDatum = NeumannDatum,
RobinDatum2 = RobinDatum2)
def energyNormMatrix(self, w:float) -> Np2D:
"""
Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
Args:
w: Weight.
Returns:
Sparse matrix in CSR format.
"""
normMatFen = FenicsHelmholtzScatteringEngine.energyNormMatrix(self, w)
return scsp.block_diag((normMatFen, 1/w * normMatFen))
def liftDirichletData(self):
"""Lift Dirichlet datum."""
lift1 = FenicsHelmholtzScatteringEngine.liftDirichletData(self)
return np.pad(lift1, (0, len(lift1)), 'constant')
def setProblemData(self, dataDict:dict, k:complex):
"""
Set problem data.
Args:
dataDict: Dictionary of problem data.
k: Parameter value.
"""
self.wavenumber = k
self.forcingTerm = dataDict["forcingTerm"]
self.DirichletDatum = dataDict["DirichletDatum"]
self.NeumannDatum = dataDict["NeumannDatum"]
self.RobinDatum2 = dataDict["RobinDatum2"]
if hasattr(self, "A0"): del self.A0
if hasattr(self, "A1"): del self.A1
if hasattr(self, "b0"): del self.b0
@property
def constraintType(self):
"""Value of constraintType."""
return self._constraintType
@constraintType.setter
def constraintType(self, constraintType):
try:
constraintType = constraintType.upper().strip().replace(" ","")
if constraintType not in ["IDENTITY", "MASS"]: raise
self._constraintType = constraintType
except:
warnings.warn(("Prescribed constraintType not recognized. "
"Overriding to 'IDENTITY'."))
self.constraintType = "IDENTITY"
def assembleA(self):
"""Assemble matrix blocks of linear system."""
builtA0 = hasattr(self, "A0")
builtA1 = hasattr(self, "A1")
if not (builtA0 and builtA1):
if builtA1 and self.constraintType == "IDENTITY":
self.A2 = None
FenicsHelmholtzScatteringEngine.assembleA(self)
if self.constraintType == "IDENTITY":
if not builtA0:
block = scsp.eye(self.A0.shape[0])
else:
block = scsp.eye(self.A1.shape[0])
else: #MASS
block = - self.A2
I = list(self.DirichletBC0.get_boundary_values().keys())
warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning)
block[I, I] = 1.
warnings.simplefilter("default", scsp.SparseEfficiencyWarning)
if not builtA0:
self.A0 = scsp.block_diag((self.A0, block))
if not builtA1:
self.A1 = scsp.bmat([[self.A1, self.A2], [- block, None]])
if hasattr(self, 'A2'): del self.A2
def assembleb(self):
"""Assemble RHS of linear system."""
if not hasattr(self, "b0"):
FenicsHelmholtzScatteringEngine.assembleb(self)
self.b0 = np.pad(self.b0, (0, self.b0.shape[0]), 'constant')
parDegree = 1
def As(self) -> List[Np2D]:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1]
def A(self) -> Np2D:
"""Assemble matrix of linear system."""
self.assembleA()
return self.A0 + self.wavenumber * self.A1
def solve(self) -> Tuple[Np1D, Np1D]:
"""
Find solution of linear system.
Returns:
Solution vectors.
"""
u = FenicsHelmholtzScatteringEngine.solve(self)
lu2 = int(len(u) / 2)
return np.array([u[: lu2], u[lu2 :]])

Event Timeline

c4science · Help