Page MenuHomec4science
Files F60754001

laplace_base_problem_engine.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 09:52

laplace_base_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 fenics as fen
from rrompy.hfengines.base.problem_engine_base import ProblemEngineBase
from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal
from rrompy.solver.fenics import (fenZERO, fenONE, L2InverseNormMatrix,
H1NormMatrix, Hminus1NormMatrix)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.parameter import checkParameter
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LaplaceBaseProblemEngine']
class LaplaceBaseProblemEngine(ProblemEngineBase):
"""
Solver for generic Laplace problems.
- \nabla \cdot (a \nabla u) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h 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.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing inner product
over V'.
dualityMatrix: Scipy sparse matrix representing duality V-V'.
energyNormPartialDualMatrix: Scipy sparse matrix representing dual
inner product between Riesz representers V-V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
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.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
_energyDualNormCompress = None
def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.mu0 = checkParameter(mu0)
self.npar = self.mu0.shape[1]
self.omega = np.abs(self.mu0(0, 0)) if self.npar > 0 else 0.
self.diffusivity = fenONE
self.forcingTerm = fenZERO
self.DirichletDatum = fenZERO
self.NeumannDatum = fenZERO
self.RobinDatumG = fenZERO
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
ProblemEngineBase.V.fset(self, V)
self.dsToBeSet = True
@property
def diffusivity(self):
"""Value of a."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
self.resetAs()
if not isinstance(diffusivity, (list, tuple,)):
diffusivity = [diffusivity, fenZERO]
self._diffusivity = diffusivity
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm, fenZERO]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
if not isinstance(DirichletDatum, (list, tuple,)):
DirichletDatum = [DirichletDatum, fenZERO]
self._DirichletDatum = DirichletDatum
@property
def NeumannDatum(self):
"""Value of g1."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
self.resetbs()
if not isinstance(NeumannDatum, (list, tuple,)):
NeumannDatum = [NeumannDatum, fenZERO]
self._NeumannDatum = NeumannDatum
@property
def RobinDatumG(self):
"""Value of g2."""
return self._RobinDatumG
@RobinDatumG.setter
def RobinDatumG(self, RobinDatumG):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG, fenZERO]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
vbMng(self, "INIT", "Initializing boundary measures.", 20)
mesh = self.V.mesh()
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = mesh,
subdomain_data = boundary_markers)
self.dsToBeSet = False
vbMng(self, "DEL", "Done assembling boundary measures.", 20)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self.energyNormMatrix = H1NormMatrix(self.V, np.abs(self.omega)**2)
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self.energyNormDualMatrix = Hminus1NormMatrix(
self.V, np.abs(self.omega)**2,
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildDualityPairingForm(self):
"""Build sparse matrix (in CSR format) representative of duality."""
vbMng(self, "INIT", "Assembling duality matrix.", 20)
self.dualityMatrix = L2InverseNormMatrix(
self.V, solverType = self._solver,
solverArgs = self._solverArgs,
compressRank = self._dualityCompress)
vbMng(self, "DEL", "Done assembling duality matrix.", 20)
def buildEnergyNormPartialDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy partial dual matrix.", 20)
self.energyNormPartialDualMatrix = Hminus1NormMatrix(
self.V, np.abs(self.omega)**2,
compressRank = self._energyDualNormCompress,
duality = False)
vbMng(self, "DEL", "Done assembling energy partial dual matrix.", 20)
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)
if derI <= 0 and self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
termNames = ["diffusivity", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip([aRe, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip([aIm, hIm],
[x + "Imag" for x in termNames]))
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))
self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
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):
if bs[j] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
termNames, terms = [], []
if j == 0:
u0Re, u0Im = self.DirichletDatum
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
termNames += ["forcingTerm", "NeumannDatum", "RobinDatumG"]
terms += [[fRe, fIm], [g1Re, g1Im], [g2Re, g2Im]]
else:
u0Re, u0Im = fenZERO, fenZERO
fRe, fIm = fenZERO, fenZERO
g1Re, g1Im = fenZERO, fenZERO
g2Re, g2Im = fenZERO, fenZERO
if len(termNames) > 0:
parsRe = self.iterReduceQuadratureDegree(zip(
[term[0] for term in terms],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[term[1] for term in terms],
[x + "Imag" for x in termNames]))
else:
parsRe, parsIm = {}, {}
L0Re = (fen.dot(fRe, self.v) * fen.dx
+ fen.dot(g1Re, self.v) * self.ds(0)
+ fen.dot(g2Re, self.v) * self.ds(1))
L0Im = (fen.dot(fIm, self.v) * fen.dx
+ fen.dot(g1Im, self.v) * self.ds(0)
+ fen.dot(g2Im, self.v) * self.ds(1))
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
b = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
if homogeneized:
Ader = self.A(0, hashI(j, self.npar))
b -= Ader.dot(self.liftedDirichletDatum)
if homogeneized:
self.bsH[j] = b
else:
self.bs[j] = b
vbMng(self, "DEL", "Done assembling forcing term.", 20)
return self._assembleb(mu, der, derI, homogeneized)

Event Timeline

c4science · Help