Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/fenics/helmholtz_problem_engine.py
c6af74908c2cmaster

Last Change
Raw File

helmholtz_problem_engine.py
View Options

#!/usr/bin/python
# 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 scipy.sparse as scsp
import fenics as fen
from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine, fenZERO
__all__ = ['HelmholtzProblemEngine']
class HelmholtzProblemEngine(HelmholtzBaseProblemEngine):
"""
Solver for Helmholtz 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 + h u = g2 on \Gamma_R
Attributes:
omega: Value of omega.
diffusivity: Value of a.
refractionIndex: Value of n.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumH: Value of h.
RobinDatumG: Value of g2.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
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.
b0: Numpy array representation of b0.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
"""
_RobinDatumH = [fenZERO, fenZERO]
def __init__(self):
super().__init__(self)
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.DirichletBoundary = "ALL"
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
if hasattr(self, "A0"): del self.A0
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
def assembleA(self):
"""Assemble matrix blocks of linear system."""
self.autoSetDS()
if not hasattr(self, "A0"):
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ hRe * fen.dot(self.u, self.v) * self.ds(1))
a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ hIm * fen.dot(self.u, self.v) * self.ds(1))
A0Re = fen.assemble(a0Re)
A0Im = fen.assemble(a0Im)
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.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if not hasattr(self, "A1"):
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
A1Re = fen.assemble(a1Re)
A1Im = fen.assemble(a1Im)
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.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size)
+ 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
shape = A1ImMat.size))
self.As = [self.A0, self.A1]
c4science · Help