Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/tests/2_hfengines/helmholtz_square_bubble_domain_problem_engine.py
d6451a36626adevelop

Last Change
Raw File

helmholtz_square_bubble_domain_problem_engine.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/>.
#
import numpy as np
from numpy.polynomial.polynomial import polyfit as fit
import fenics as fen
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import fenZERO
from rrompy.hfengines.fenics_engines import HelmholtzProblemEngine
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
from rrompy.parameter import parameterMap as pMap
class HelmholtzSquareBubbleDomainProblemEngine(HelmholtzProblemEngine):
"""
Solver for square bubble Helmholtz problems with parametric domain heigth.
- \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi]
u = 0 on \Gamma_mu = \partial\Omega_mu
with exact solution square bubble times plane wave.
"""
def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [1.],
*args, **kwargs):
super().__init__(mu0, *args, **kwargs)
self._affinePoly = False
self.nAs, self.nbs = 2, 15
self.kappa = kappa
self.theta = theta
mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(np.pi, np.pi),
3 * n, 3 * n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.parameterMap = pMap(1.)
def buildA(self):
"""Build terms of operator of linear system."""
if self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
self.thAs[0] = self.getMonomialSingleWeight([0])
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[1] is None:
vbMng(self, "INIT", "Assembling operator term A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2)
k2n2Re = k2Re * n2Re - k2Im * n2Im
k2n2Im = k2Re * n2Im + k2Im * n2Re
parsRe = self.iterReduceQuadratureDegree(zip([k2n2Re],
["kappaSquaredRefractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([k2n2Im],
["kappaSquaredRefractionIndexSquaredImag"]))
a2Re = (fen.dot(self.u.dx(0), self.v.dx(0))
- k2n2Re * fen.dot(self.u, self.v)) * fen.dx
a2Im = - k2n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0))
self.thAs[1] = self.getMonomialSingleWeight([2])
vbMng(self, "DEL", "Done assembling operator term.", 20)
def buildb(self):
"""Build terms of operator of linear system."""
if self.thbs[0] is None:
self.thbs = self.getMonomialWeights(self.nbs)
bDEIMCoeffs = None
for j in range(self.nbs):
if self.bs[j] is None:
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
if bDEIMCoeffs is None:
muDEIM = np.linspace(.5, 4., self.nbs)
for jj, muD in enumerate(muDEIM):
pi = np.pi
c, s = np.cos(self.theta), np.sin(self.theta)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
muR, muI = np.real(muD), np.imag(muD)
mu2R, mu2I = np.real(muD ** 2.), np.imag(muD ** 2.)
C = 16. / pi ** 4.
bR = C * (2 * (x * (pi - x) + y * (pi - y))
+ (self.kappa * s) ** 2. * (mu2R - 1.)
* x * (pi - x) * y * (pi - y))
bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y
* (pi - y) + s * x * (pi - x)
* (pi - 2 * y))
+ (self.kappa * s) ** 2. * mu2I
* x * (pi - x) * y * (pi - y))
wR = (fen.cos(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
wI = (fen.sin(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
fRe, fIm = bR * wR + bI * wI, bI * wR - bR * wI
fRe = mu2R * fRe - mu2I * fIm + fenZERO
fIm = mu2R * fIm + mu2I * fRe + fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([fRe],
["forcingTerm{}Real".format(jj)]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm],
["forcingTerm{}Imag".format(jj)]))
LR = fen.dot(fRe, self.v) * fen.dx
LI = fen.dot(fIm, self.v) * fen.dx
DBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
bjj = (fenics2Vector(LR, parsRe, DBC0, 1)
+ 1.j * fenics2Vector(LI, parsIm, DBC0, 1))
if jj == 0:
bDEIM = np.empty((self.nbs, len(bjj)),
dtype = np.complex)
bDEIM[jj] = bjj
bDEIMCoeffs = fit(muDEIM, bDEIM, self.nbs - 1)
self.bs[j] = bDEIMCoeffs[j]
vbMng(self, "DEL", "Done assembling forcing term.", 20)
c4science · Help