Page MenuHomec4science
Files F90355243

FenicsHelmholtzEngine.py
No OneTemporary
Actions

File Metadata

Created
Thu, Oct 31, 21:54

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
PI = np.pi
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:"expression",
degree:int) -> "Fenics Expression":
"""
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 copy(valueRe), copy(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:"Fenics mesh",
FEDegree:int,
wavenumber:complex,
diffusivity : "custom expression" = 1,
refractionIndex : "custom expression" = 1,
forcingTerm : "custom expression" = fenZEROC,
DirichletBoundary : "bool function or str" = None,
NeumannBoundary : "bool function or str" = None,
RobinBoundary : "bool function or str" = None,
DirichletDatum : "custom expression" = fenZEROC,
NeumannDatum : "custom expression" = fenZEROC,
RobinDatum : "custom expression" = fenZEROC,
RobinDatum2 : "custom expression" = 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 = copy(DirichletB)
self.NeumannBoundary = copy(NeumannB)
self.RobinBoundary = copy(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.FacetFunction("size_t", mesh)
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 energyNormMatrix(self, w:float) -> "CSC sparse matrix":
"""
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**2 * 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):
"""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:dict, 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):
"""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 = copy(diffusivity)
@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 = copy(refractionIndex)
@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 = copy(forcingTerm)
@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 = copy(DirichletDatum)
@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 = copy(NeumannDatum)
@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 = copy(RobinDatum)
@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 = copy(RobinDatum2)
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.array()[:] + 1.j * b0Im.array()[:],
dtype = np.complex)
parDegree = 1
def As(self) -> list:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1]
def bs(self) -> list:
"""RHS blocks of linear system."""
self.assembleb()
return [self.b0]
def A(self) -> "CSR scipy matrix":
"""Assemble matrix of linear system."""
self.assembleA()
return self.A0 + self.wavenumber2 * self.A1
def b(self) -> "numpy 1D array":
"""Assemble RHS of linear system."""
self.assembleb()
return self.b0
def solve(self) -> "numpy 1D array":
"""
Find solution of linear system.
Returns:
Solution vector.
"""
return scsp.linalg.spsolve(self.A(), self.b())
def getLSBlocks(self) -> (list, list):
"""
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.
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:"Fenics mesh",
FEDegree:int,
wavenumber:"custom complex",
diffusivity : "custom expression" = 1,
refractionIndex : "custom expression" = 1,
forcingTerm : "custom expression" = fenZEROC,
DirichletBoundary : "bool function or str" = None,
NeumannBoundary : "bool function or str" = None,
RobinBoundary : "bool function or str" = None,
DirichletDatum : "custom expression" = fenZEROC,
NeumannDatum : "custom expression" = fenZEROC,
RobinDatum2 : "custom expression" = 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 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, "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._RobinDatum = copy(RobinDatum)
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:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1, self.A2]
def A(self) -> "CSR scipy matrix":
"""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:"Fenics mesh",
FEDegree:int,
wavenumber:"custom complex",
diffusivity : "custom expression" = 1,
refractionIndex : "custom expression" = 1,
forcingTerm : "custom expression" = fenZEROC,
DirichletBoundary : "bool function or str" = None,
NeumannBoundary : "bool function or str" = None,
RobinBoundary : "bool function or str" = None,
DirichletDatum : "custom expression" = fenZEROC,
NeumannDatum : "custom expression" = fenZEROC,
RobinDatum2 : "custom expression" = 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) -> "CSC sparse matrix":
"""
Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
Args:
w: Weight.
Returns:
Sparse matrix in CSR format.
"""
normMatFen = FenicsHelmholtzEngine.energyNormMatrix(self, w)
return scsp.block_diag((normMatFen, 1/w * normMatFen))
def liftDirichletData(self):
"""Lift Dirichlet datum."""
lift1 = FenicsHelmholtzScatteringEngine.liftDirichletData(self)
return np.array([lift1, np.zeros_like(lift1)])
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:
"""Matrix blocks of linear system."""
self.assembleA()
return [self.A0, self.A1]
def A(self) -> "CSR scipy matrix":
"""Assemble matrix of linear system."""
self.assembleA()
return self.A0 + self.wavenumber * self.A1
def solve(self) -> ("numpy 1D array", "numpy 1D array"):
"""
Find solution of linear system.
Returns:
Solution vector.
"""
u = FenicsHelmholtzScatteringEngine.solve(self)
lu2 = int(len(u) / 2)
return np.array([u[: lu2], u[lu2 :]])

Event Timeline

c4science · Help