Page MenuHomec4science
Files F62204366

helmholtz_square_domain_problem_engine.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 11, 14:55

helmholtz_square_domain_problem_engine.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import scipy.sparse as scsp
import fenics as fen
from rrompy.utilities.base.types import (ScOp, List, Tuple, paramVal, Np1D,
FenExpr)
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
__all__ = ['HelmholtzSquareDomainProblemEngine']
class HelmholtzSquareDomainProblemEngine(HelmholtzProblemEngine):
"""
Solver for square Helmholtz problems with parametric laplacian.
- \dxx u - mu_2**2 \dyy u - mu_1**2 * u = f(mu_2) in \Omega = [0,\pi]^2
u = 0 on \partial\Omega
"""
def __init__(self, kappa:float, theta:float, n:int,
mu0 : paramVal = [12. ** .5, 1.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.nAs, self.nbs = 6, 11 * 12 // 2
self.npar = 2
self.rescalingExp = [2., 1.]
pi = np.pi
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
self.forcingTerm = [fen.cos(kappa * (c * x + s / self.mu0(0, 1) * y)),
fen.sin(kappa * (c * x + s / self.mu0(0, 1) * y))]
self.forcingTermDer = kappa * s * y
def getExtraFactorB(self, mu : paramVal = [],
derI : int = 0) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
mu = self.checkParameter(mu)
def getPowMinusj(x, power):
powR = x ** power
powI = fenZERO
if power % 2 == 1:
powR, powI = powI, powR
if power % 4 > 1:
powR, powI = - powR, - powI
return powR, powI
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
if derI == 0: return fenONE, fenZERO
coeffs = np.zeros(derI + 1)
coeffs[1] = - 2.
for j in range(2, derI + 1):
coeffs[1 :] = (-2. / j * coeffs[1 :]
- (3 - (1 + 2 * np.arange(derI)) / j) * coeffs[: -1])
for j in range(derI):
powR, powI = getPowMinusj(self.forcingTermDer, derI - j)
mupBase = coeffs[j + 1] * mu(0, 1) ** (- 3 * derI + 2 * j)
mupR, mupI = np.real(mupBase), np.imag(mupBase)
if j == 0:
exprR = mupR * powR - mupI * powI
exprI = mupI * powR + mupR * powI
else:
exprR += mupR * powR - mupI * powI
exprI += mupI * powR + mupR * powI
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return exprR, exprI
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
for j in range(2, 5):
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
A0Re = fen.assemble(a0Re)
DirichletBC0.apply(A0Re)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size,
dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
DirichletBC0.zero(A1Re)
DirichletBC0.zero(A1Im)
A1ReMat = fen.as_backend_type(A1Re).mat()
A1ImMat = fen.as_backend_type(A1Im).mat()
A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size)
+ 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
shape = A1ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 5 and self.As[5] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A5.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a5Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
A5Re = fen.assemble(a5Re)
DirichletBC0.zero(A5Re)
A5ReMat = fen.as_backend_type(A5Re).mat()
A5Rer, A5Rec, A5Rev = A5ReMat.getValuesCSR()
self.As[5] = scsp.csr_matrix((A5Rev, A5Rec, A5Rer),
shape = A5ReMat.size,
dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
nbsTot = self.nbsH if homogeneized else self.nbs
bs = self.bsH if homogeneized else self.bs
if homogeneized and self.mu0 != self.mu0BC:
self.liftDirichletData(self.mu0)
for j in range(derI, nbsTot):
derH = hashI(j, self.npar)
if bs[j] is None:
if np.sum(derH) != derH[-1]:
if homogeneized:
self.bsH[j] = self.checkbInBounds(-1)
else:
self.bs[j] = self.checkbInBounds(-1)
continue
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(j),
timestamp = self.timestamp)
if j == 0:
u0Re, u0Im = self.DirichletDatum
else:
u0Re, u0Im = fenZERO, fenZERO
if j < self.nbs:
fRe, fIm = self.forcingTerm
cRe, cIm = self.getExtraFactorB(self.mu0, derH[-1])
cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(j)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(j)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
DBCR.apply(b0Re)
DBCI.apply(b0Im)
b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
if homogeneized:
Ader = self.A(self.mu0, derH)
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return self._assembleb(mu, der, derI, homogeneized, self.mu0)

Event Timeline

c4science · Help