Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/hfengines/scipy/helmholtz_square_bubble_domain_problem_engine.py
0c9a0bbd9f4edevelop

Last Change
Raw File

helmholtz_square_bubble_domain_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 scipy.sparse as scsp
import fenics as fen
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, FenExpr
from rrompy.hfengines.scipy.helmholtz_problem_engine import (
HelmholtzProblemEngine, fenZERO)
__all__ = ['HelmholtzSquareBubbleDomainProblemEngine']
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):
super().__init__()
self.omega = kappa
self.kappa = kappa
self.theta = theta
self.derMaxB = 20
self.bs = [None] * self.derMaxB
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
self.V = fen.FunctionSpace(mesh, "P", 3)
def getForcingTerm(self, mu:complex) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
if (not hasattr(self, "forcingTermMu")
or not np.isclose(mu, self.forcingTermMu)):
import sympy as sp
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 16/np.pi**4 * x * y * (x - np.pi) * (y - np.pi)
phiex = self.kappa * (x * np.cos(self.theta) + y * mu * np.sin(self.theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - self.kappa**2 * uex
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
sp.printing.ccode(sp.simplify(sp.im(fex)))]
self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm]
self.forcingTermMu = mu
return self.forcingTerm
def getExtraFactorB(self, mu:complex, der:int) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
import sympy as sp
from math import factorial as fact
y = sp.symbols('x[1]', real=True)
expr = mu**2.*np.power(-1.j * self.kappa * np.sin(self.theta) * y, der)
if der >= 1:
expr += 2.*mu*np.power(-1.j * self.kappa * np.sin(self.theta) * y,
der - 1) * der
if der >= 2:
expr += np.power(-1.j * self.kappa * np.sin(self.theta) * y,
der - 2) * der * (der - 1)
extraTerm = [sp.printing.ccode(sp.simplify(sp.re(expr / fact(der)))),
sp.printing.ccode(sp.simplify(sp.im(expr / fact(der))))]
return [fen.Expression(x, degree = 3) for x in extraTerm]
def rescaling(self, x:Np1D) -> Np1D:
"""Rescaling in parameter dependence."""
return x
def rescalingInv(self, x:Np1D) -> Np1D:
"""Inverse rescaling in parameter dependence."""
return x
def A(self, mu:complex, der : int = 0) -> Np2D:
"""Assemble (derivative of) operator of linear system."""
if der > 2 or der < 0:
d = self.V.dim()
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
self.autoSetDS()
if der <= 0 and not hasattr(self, "A0"):
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
A0Re = fen.assemble(a0Re)
DirichletBC0.apply(A0Re)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
self.A0 = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size,
dtype = np.complex)
if der <= 2 and not hasattr(self, "A2"):
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
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
A2Re = fen.assemble(a2Re)
A2Im = fen.assemble(a2Im)
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.A2 = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
shape = A2ReMat.size)
+ 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
shape = A2ImMat.size))
if der == 0:
return self.A0 + mu ** 2 * self.A2
if der == 1:
return 2. * mu * self.A2
return self.A2
def b(self, mu:complex, der : int = 0) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
if der < 0 or der >= self.derMaxB:
return np.zeros(self.V.dim(), dtype = np.complex)
if self.bs[der] is None:
fRe, fIm = self.getForcingTerm(mu)
cRe, cIm = self.getExtraFactorB(mu, der)
cfRe = cRe * fRe - cIm * fIm
cfIm = cRe * fIm + cIm * fRe
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
b0Re = fen.assemble(L0Re)
b0Im = fen.assemble(L0Im)
if der == 0:
DirichletBCRe = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DirichletBCIm = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
else:
DirichletBCRe = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
DirichletBCIm = DirichletBCRe
DirichletBCRe.apply(b0Re)
DirichletBCIm.apply(b0Im)
self.bs[der] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
return self.bs[der]
def affineBlocksA(self, mu : complex = 0.) -> Tuple[List[Np2D], callable]:
"""Assemble affine blocks of operator of linear system."""
def lambdas(x, j):
if j == 0: return np.ones(np.size(x))
if j in [1, 2]: return np.power(self.rescaling(x)
- self.rescaling(mu), j)
raise Exception("Wrong j value.")
A0 = self.A(mu, 0)
A1 = self.A(mu, 1)
return [A0, A1, self.A2], lambdas
def affineBlocksb(self, mu : complex = 0.) -> Tuple[List[Np1D], callable]:
"""Assemble affine blocks of RHS of linear system."""
def lambdas(x, j):
if j == 0: return np.ones(np.size(x))
if j in range(1, self.derMaxB): return np.power(self.rescaling(x)
- self.rescaling(mu),
j)
raise Exception("Wrong j value.")
for j in range(self.derMaxB):
self.b(mu, j)
return self.bs, lambdas
c4science · Help