Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/fenics_engines/helmholtz_problem_engine_augmented.py
74dce5ca631dmaster

Last Change
Raw File

helmholtz_problem_engine_augmented.py
View Options

# Copyright (C) 2018-2020 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 pad
from scipy.sparse import eye, bmat, block_diag
from collections.abc import Iterable
from .helmholtz_problem_engine import (HelmholtzProblemEngine,
ScatteringProblemEngine)
from rrompy.solver.fenics import (augmentedH1NormMatrix,
augmentedHminus1NormMatrix)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import parameterMap as pMap
__all__ = ['HelmholtzProblemEngineAugmented',
'ScatteringProblemEngineAugmented']
class HelmholtzProblemEngineAugmented(HelmholtzProblemEngine):
"""
Solver for generic Helmholtz problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega * n**2 * v = f in \Omega
omega * u = v in \overline{\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.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
refractionIndex: Value of n.
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).
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.parameterMap = pMap(1., self.npar)
@property
def spacedim(self):
if (hasattr(self, "bs") and isinstance(self.bs, Iterable)
and self.bs[0] is not None): return len(self.bs[0])
return 2 * super().spacedim
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self._energyNormMatrix = augmentedH1NormMatrix(self.V)
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self._energyNormDualMatrix = augmentedHminus1NormMatrix(self.V,
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
ANone = any([A is None for A in self.As])
if not ANone: return
self.nAs = 2
super().buildA()
I = eye(self.spacedim // 2)
self.As[0] = block_diag((self.As[0], I), format = "csr")
self.As[1] = bmat([[None, self.As[1]], [- I, None]], format = "csr")
def buildb(self):
"""Build terms of operator of linear system."""
bNone = any([b is None for b in self.bs])
if not bNone: return
self.nbs = 1
dim = self.spacedim // 2
super().buildb()
self.bs[0] = pad(self.bs[0], (0, dim), "constant")
def plot(self, u, warping = None, is_state = False, name = "u",
save = None, what = 'all', forceNewFile = True,
saveFormat = "eps", saveDPI = 100, show = True, colorMap = "jet",
fenplotArgs = {}, **figspecs):
uh = u[: self.spacedim // 2] if is_state or self.isCEye else u
return super().plot(uh, warping, is_state, name, save, what,
forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
def outParaview(self, u, warping = None, is_state = False, name = "u",
filename = "out", time = 0., what = 'all',
forceNewFile = True, folder = False, filePW = None):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaview(u[: self.spacedim // 2], warping, is_state,
name, filename, time, what, forceNewFile,
folder, filePW)
def outParaviewTimeDomain(self, u, omega, warping = None, is_state = False,
timeFinal = None, periodResolution = 20,
name = "u", filename = "out",
forceNewFile = True, folder = False):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaviewTimeDomain(u[: self.spacedim // 2], omega,
warping, is_state, timeFinal,
periodResolution, name, filename,
forceNewFile, folder)
class ScatteringProblemEngineAugmented(ScatteringProblemEngine):
"""
Solver for scattering problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega * n**2 * v = f in \Omega
omega * u = v in \overline{\Omega}
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu +- i v = 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.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
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).
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.nAs = 2
self._weight0 = 1.
@property
def spacedim(self):
if (hasattr(self, "bs") and isinstance(self.bs, Iterable)
and self.bs[0] is not None): return len(self.bs[0])
return 2 * super().spacedim
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self._energyNormMatrix = augmentedH1NormMatrix(self.V)
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self._energyNormDualMatrix = augmentedHminus1NormMatrix(self.V,
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
ANone = any([A is None for A in self.As])
if not ANone: return
self.nAs = 3
super().buildA()
self._nAs = 2
I = eye(self.spacedim // 2)
self.As[0] = bmat([[self.As[0], self._weight0 * self.As[1]],
[None, I]], format = "csr")
self.As[1] = bmat([[(1. - self._weight0) * self.As[1], self.As[2]],
[- I, None]], format = "csr")
self.thAs.pop()
self.As.pop()
def buildb(self):
"""Build terms of operator of linear system."""
bNone = any([b is None for b in self.bs])
if not bNone: return
self.nbs = 1
dim = self.spacedim // 2
super().buildb()
self.bs[0] = pad(self.bs[0], (0, dim), "constant")
def plot(self, u, warping = None, is_state = False, name = "u",
save = None, what = 'all', forceNewFile = True,
saveFormat = "eps", saveDPI = 100, show = True, colorMap = "jet",
fenplotArgs = {}, **figspecs):
uh = u[: self.spacedim // 2] if is_state or self.isCEye else u
return super().plot(uh, warping, is_state, name, save, what,
forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
def outParaview(self, u, warping = None, is_state = False, name = "u",
filename = "out", time = 0., what = 'all',
forceNewFile = True, folder = False, filePW = None):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaview(u[: self.spacedim // 2], warping, is_state,
name, filename, time, what, forceNewFile,
folder, filePW)
def outParaviewTimeDomain(self, u, omega, warping = None, is_state = False,
timeFinal = None, periodResolution = 20,
name = "u", filename = "out",
forceNewFile = True, folder = False):
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
return super().outParaviewTimeDomain(u[: self.spacedim // 2], omega,
warping, is_state, timeFinal,
periodResolution, name, filename,
forceNewFile, folder)
c4science · Help