Page MenuHomec4science
Files F60948395

cookie_engine_single.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 13:35

cookie_engine_single.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
import ufl
from rrompy.utilities.base.types import ScOp, List, paramVal
from rrompy.solver.fenics import fenONE, fenZERO
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenics2Sparse
__all__ = ['CookieEngineSingle']
class CookieEngineSingle(HelmholtzProblemEngine):
def __init__(self, kappa:float, theta:float, n:int, R : int = 1.,
L : int = 2., nX : int = 1, nY : int = 1,
mu0 : paramVal = [12. ** .5, 1.],
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.nAs = 5
self.npar = 2
self.rescalingExp = [2., 2.]
mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(L * nX, L * nY),
n * nX, n * nY)
self.V = fen.FunctionSpace(mesh, "P", 1)
x, y = fen.SpatialCoordinate(mesh)[:]
cxs = np.linspace(0, L * nX, 2 * nX + 1)[1::2]
cys = np.linspace(0, L * nY, 2 * nY + 1)[1::2]
self.cookieIn = fenZERO
for cx in cxs:
for cy in cys:
self.cookieIn += ufl.conditional(
ufl.le((x-cx)**2. + (y-cy)**2., R**2.),
fenONE, fenZERO)
self.cookieOut = fenONE - self.cookieIn
c, s = np.cos(theta), np.sin(theta)
self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
fen.sin(kappa * (c * x + s * y))]
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)
self.autoSetDS()
if derI <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
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))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if derI <= 2 and self.As[2] is None:
self.As[2] = self.checkAInBounds(-1)
if derI <= 3 and self.As[3] is None:
self.As[3] = self.checkAInBounds(-1)
if derI <= 4 and self.As[4] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A4.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a4Re = - n2Re * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
a4Im = - n2Im * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
self.As[4] = (fenics2Sparse(a4Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a4Im, parsIm, DirichletBC0, 0))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self._assembleA(mu, der, derI)

Event Timeline

c4science · Help