Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_inverse_problem_engine.py
c11471199317master

Last Change
Raw File

helmholtz_square_bubble_domain_inverse_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.hfengines.fenics.helmholtz_problem_engine import (
HelmholtzProblemEngine, fenZERO)
__all__ = ['HelmholtzSquareBubbleDomainInverseProblemEngine']
class HelmholtzSquareBubbleDomainInverseProblemEngine(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.
"""
bterms = 2
def __init__(self, kappa:float, theta:float, n:int):
super().__init__()
self.omega = kappa
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
self.V = fen.FunctionSpace(mesh, "P", 3)
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 = kappa * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - 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]
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
a0Re = aRe * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
a0Im = aIm * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
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)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
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
a1Re = ((aRe * fen.dot(self.u.dx(0), self.v.dx(0))
- k2n2Re * fen.dot(self.u, self.v)) * fen.dx
+ hRe * fen.dot(self.u, self.v) * self.ds(1))
a1Im = ((aIm * fen.dot(self.u.dx(0), self.v.dx(0))
- k2n2Im * fen.dot(self.u, self.v)) * fen.dx
+ hIm * fen.dot(self.u, self.v) * self.ds(1))
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]
def assembleb(self):
"""Assemble RHS blocks of linear system."""
if not hasattr(self, "b1"):
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
L1Re = (fen.dot(fRe, self.v) * fen.dx
+ fen.dot(g1Re, self.v) * self.ds(0)
+ fen.dot(g2Re, self.v) * self.ds(1))
L1Im = (fen.dot(fIm, self.v) * fen.dx
+ fen.dot(g1Im, self.v) * self.ds(0)
+ fen.dot(g2Im, self.v) * self.ds(1))
b1Re = fen.assemble(L1Re)
b1Im = fen.assemble(L1Im)
fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary).apply(b1Re)
fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary).apply(b1Im)
self.b1 = np.array(b1Re[:] + 1.j * b1Im[:], dtype = np.complex)
if not hasattr(self, "b0"):
self.b0 = np.zeros_like(self.b1, dtype = np.complex)
self.bs = [self.b0, self.b1]
c4science · Help