Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/hfengines/linear_problem/scattering_problem_engine.py
0c9d322ca745feature_greedy_coarseni…

Last Change
Raw File

scattering_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/>.
#
from numpy import inf
import scipy.sparse as scsp
import fenics as fen
from rrompy.utilities.base.types import ScOp
from rrompy.utilities.base.fenics import fenZERO
from rrompy.utilities.base import verbosityDepth
from .helmholtz_problem_engine import HelmholtzProblemEngine
from rrompy.utilities.exception_manager import RROMPyWarning
__all__ = ['ScatteringProblemEngine']
class ScatteringProblemEngine(HelmholtzProblemEngine):
"""
Solver for scattering problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu +- i omega u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
signR: Sign in ABC.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
A1: Scipy sparse array representation (in CSC format) of A1.
A2: Scipy sparse array representation (in CSC format) of A2.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
nAs = 3
rescalingExp = 1.
signR = - 1.
def __init__(self, degree_threshold : int = inf, verbosity : int = 10,
timestamp : bool = True):
self.silenceWarnings = True
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
del self.silenceWarnings
@property
def RobinDatumH(self):
"""Value of h."""
return self.signR * self.omega
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
if not hasattr(self, "silenceWarnings"):
RROMPyWarning(("Scattering problems do not allow changes of h. "
"Ignoring assignment."))
return
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 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)
aRe, aIm = self.diffusivity
parsRe = self.iterReduceQuadratureDegree(zip([aRe],
["diffusivityReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([aIm],
["diffusivityImag"]))
a0Re = aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
a0Im = aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
DirichletBC0.apply(A0Re)
DirichletBC0.zero(A0Im)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0ImMat = fen.as_backend_type(A0Im).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 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)
a1 = fen.dot(self.u, self.v) * self.ds(1)
A1 = fen.assemble(a1)
DirichletBC0.zero(A1)
A1Mat = fen.as_backend_type(A1).mat()
A1r, A1c, A1v = A1Mat.getValuesCSR()
self.As[1] = self.signR * 1.j * scsp.csr_matrix((A1v, A1c, A1r),
shape = A1Mat.size)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
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"]))
a2Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a2Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
DirichletBC0.zero(A2Re)
DirichletBC0.zero(A2Im)
A2ReMat = fen.as_backend_type(A2Re).mat()
A2ImMat = fen.as_backend_type(A2Im).mat()
A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
self.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
shape = A2ReMat.size)
+ 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
shape = A2ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu * self.As[1] + mu**2. * self.As[2]
if der == 1:
return self.As[1] + 2 * mu * self.As[2]
return self.As[2]
c4science · Help