Page MenuHomec4science
Files F61827797

FreeFemHelmholtzEngine.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 9, 05:32

FreeFemHelmholtzEngine.py
View Options

#!/usr/bin/python
from copy import copy
import os
import warnings
import subprocess
import numpy as np
import scipy.sparse as scsp
import scipy.sparse.linalg as spla
import fenics as fen
import FreeFemConversionTools as FFCT
PI = np.pi
fenZERO = fen.Constant(0)
fenZEROC = [fenZERO, fenZERO]
class DomainError(Exception):
"""
Domain error exception.
Args:
value: Human readable string describing the exception.
Attributes:
value: Human readable string describing the exception.
"""
def __init__(self, value:str):
self.value = value
def __str__(self):
return self.value
class FreeFemHelmholtzEngine:
"""
FreeFem++-based solver for Helmholtz problems.
- \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Args:
Attributes:
"""
def __init__(self, V:str,
wavenumber:complex,
compl : bool = True,
dim : int = 2,
meshname : str = "Th",
Vname : str = "V",
diffusivity : str = "1",
refractionIndex : str = "1",
forcingTerm : str = "0",
DirichletBoundary : str = "",
NeumannBoundary : str = "",
RobinBoundary : str = "",
DirichletDatum : str = "0",
NeumannDatum : str = "0",
RobinDatum : str = "0",
RobinDatum2 : str = "0"):
# process boundaries
boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary]
for i in range(3):
boundariesList[i] = boundariesList[i].strip().upper()
DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList
if boundariesList.count("") == 3:
raise DomainError("At least one boundary must be prescribed.")
self.DirichletBoundary = DirichletBoundary
self.NeumannBoundary = NeumannBoundary
self.RobinBoundary = RobinBoundary
if compl:
self.compl = "<complex>"
else:
self.compl = ""
self.V = V
self.dim = dim
self.meshname = meshname
self.Vname = Vname
self.wavenumber2 = wavenumber ** 2
self.diffusivity = diffusivity
self.refractionIndex = refractionIndex
self.forcingTerm = forcingTerm
self.DirichletDatum = DirichletDatum
self.NeumannDatum = NeumannDatum
self.RobinDatum = RobinDatum
self.RobinDatum2 = RobinDatum2
def energyNormMatrix(self, w:float) -> "CSC sparse matrix":
"""
Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
Args:
w: Weight.
Returns:
Sparse matrix in CSR format.
"""
filename = FFCT.getNewFilename("temp", "edp")
Mfile = FFCT.getNewFilename("m", "tmp")
with open(filename, 'w') as f:
f.write(self.V + "\n")
f.write("{}{} u, v;\n".format(self.Vname, self.compl))
f.write("varf energy(u, v) = int{}d({})({})".format(self.dim,
self.meshname,
FFCT.diff(dim = self.dim,
compl = not self.compl == "")))
f.write(" + {} * int{}d({})(abs(u)^2);\n".format(w ** 2,
self.dim,
self.meshname))
f.write("matrix{0} M = energy({1}, {1});\n".format(self.compl,
self.Vname))
f.write("ofstream Mout(\"{}\");\n".format(Mfile))
f.write("Mout << M;\n")
bashCommand = "FreeFem++ -v 0 {}".format(filename)
process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
M = FFCT.ffmat2spmat(Mfile)
os.remove(filename)
os.remove(Mfile)
return M.tocsr()
def problemData(self):
"""List of HF problem data."""
dataDict = {}
dataDict["forcingTerm"] = self.forcingTerm
dataDict["DirichletDatum"] = self.DirichletDatum
dataDict["NeumannDatum"] = self.NeumannDatum
dataDict["RobinDatum2"] = self.RobinDatum2
return dataDict
def setProblemData(self, dataDict:dict, k2:complex):
"""
Set problem data.
Args:
dataDict: Dictionary of problem data.
k2: Parameter value.
"""
self.wavenumber2 = k2
self.forcingTerm = dataDict["forcingTerm"]
self.DirichletDatum = dataDict["DirichletDatum"]
self.NeumannDatum = dataDict["NeumannDatum"]
self.RobinDatum2 = dataDict["RobinDatum2"]
def getLSBlocks(self) -> ("numpy 2D array" * 3, "numpy 1D array"):
"""
Get blocks of linear system.
Returns:
Blocks of system (\sum_{j=0}^J k^j A_j = b)
"""
self.assembleA()
self.assembleb()
return [self.AL + self.AR, - self.AM], [self.b]
def liftDirichletData(self):
"""Lift Dirichlet datum."""
raise Exception("Not implemented.")
return
def setupDerivativeComputation(self, j:int, up:"numpy 1D array",
upp : "numpy 1D array" = None):
"""
Setup problem data to compute solution derivatives.
Args:
j: Derivative index.
up: Adjusted previous derivative.
upp: Adjusted pre-previous derivative.
"""
raise Exception("Not implemented.")
upRe = fen.Function(self.V)
upIm = fen.Function(self.V)
upRe.vector()[:] = np.array(np.real(up), dtype = float)
upIm.vector()[:] = np.array(np.imag(up), dtype = float)
upRe, upIm = self.n2Re * upRe - self.n2Im * upIm,\
self.n2Re * upIm + self.n2Im * upRe
self.forcingTerm = [upRe, upIm]
self.DirichletDatum = fenZEROC
self.NeumannDatum = fenZEROC
self.RobinDatum2 = fenZEROC
@property
def wavenumber2(self):
"""Value of k^2."""
return self._wavenumber2
@wavenumber2.setter
def wavenumber2(self, wavenumber2):
if hasattr(self, "A"): del self.A
if hasattr(self, "u"): del self.u
self._wavenumber2 = wavenumber2
@property
def diffusivity(self):
"""Value of a. Its assignment changes aRe and aIm."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
if hasattr(self, "A"): del self.A
if hasattr(self, "AL"): del self.AL
if hasattr(self, "u"): del self.u
self.a = diffusivity
self._diffusivity = copy(diffusivity)
@property
def refractionIndex(self):
"""Value of n. Its assignment changes nRe, nIm, n2Re and n2Im."""
return self._refractionIndex
@refractionIndex.setter
def refractionIndex(self, refractionIndex):
if hasattr(self, "A"): del self.A
if hasattr(self, "AM"): del self.AM
if hasattr(self, "u"): del self.u
self.n = refractionIndex
self._refractionIndex = copy(refractionIndex)
@property
def forcingTerm(self):
"""Value of f. Its assignment changes fRe and fIm."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
if hasattr(self, "b"): del self.b
if hasattr(self, "bF"): del self.bF
if hasattr(self, "u"): del self.u
self.f = forcingTerm
self._forcingTerm = copy(forcingTerm)
@property
def DirichletDatum(self):
"""
Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
DirichletBCIm.
"""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
if hasattr(self, "b"): del self.b
if hasattr(self, "u"): del self.u
self.u0 = DirichletDatum
self._DirichletDatum = copy(DirichletDatum)
@property
def NeumannDatum(self):
"""Value of g1. Its assignment changes g1Re and g1Im."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
if hasattr(self, "b"): del self.b
if hasattr(self, "bN"): del self.bN
if hasattr(self, "u"): del self.u
self.g1 = NeumannDatum
self._NeumannDatum = copy(NeumannDatum)
@property
def RobinDatum(self):
"""Value of h. Its assignment changes hRe and hIm."""
return self._RobinDatum
@RobinDatum.setter
def RobinDatum(self, RobinDatum):
if hasattr(self, "A"): del self.A
if hasattr(self, "AR"): del self.AR
if hasattr(self, "u"): del self.u
self.h = RobinDatum
self._RobinDatum = copy(RobinDatum)
@property
def RobinDatum2(self):
"""Value of g2. Its assignment changes g2Re and g2Im."""
return self._RobinDatum2
@RobinDatum2.setter
def RobinDatum2(self, RobinDatum2):
if hasattr(self, "b"): del self.b
if hasattr(self, "bR"): del self.bR
if hasattr(self, "u"): del self.u
self.g2 = RobinDatum2
self._RobinDatum2 = copy(RobinDatum2)
def assembleA(self, Rfact : complex = 1.):
"""Assemble matrix of linear system."""
raise Exception("Not implemented.")
if not hasattr(self, "A"):
if not hasattr(self, "AL"):
aLRe = self.aRe * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
aLIm = self.aIm * fen.dot(fen.grad(self.v1),
fen.grad(self.v2)) * fen.dx
ALRe = fen.assemble(aLRe)
ALIm = fen.assemble(aLIm)
self.DirichletBC0.apply(ALRe)
self.DirichletBC0.zero(ALIm)
ALReMat = fen.as_backend_type(ALRe).mat()
ALRer, ALRec, ALRev = ALReMat.getValuesCSR()
ALImr, ALImc, ALImv = fen.as_backend_type(
ALIm).mat().getValuesCSR()
self.AL = (scsp.csr_matrix((ALRev, ALRec, ALRer),
shape = ALReMat.size)
+ 1.j * scsp.csr_matrix((ALImv, ALImc, ALImr),
shape = ALReMat.size))
if not hasattr(self, "AM"):
aMRe = self.n2Re * fen.dot(self.v1, self.v2) * fen.dx
aMIm = self.n2Im * fen.dot(self.v1, self.v2) * fen.dx
AMRe = fen.assemble(aMRe)
AMIm = fen.assemble(aMIm)
self.DirichletBC0.zero(AMRe)
self.DirichletBC0.zero(AMIm)
AMRer, AMRec, AMRev = fen.as_backend_type(
AMRe).mat().getValuesCSR()
AMImr, AMImc, AMImv = fen.as_backend_type(
AMIm).mat().getValuesCSR()
self.AM = (scsp.csr_matrix((AMRev, AMRec, AMRer),
shape = self.AL.shape)
+ 1.j * scsp.csr_matrix((AMImv, AMImc, AMImr),
shape = self.AL.shape))
if not hasattr(self, "AR"):
aRRe = self.hRe * fen.dot(self.v1, self.v2) * self.ds(1)
aRIm = self.hIm * fen.dot(self.v1, self.v2) * self.ds(1)
ARRe = fen.assemble(aRRe)
ARIm = fen.assemble(aRIm)
self.DirichletBC0.zero(ARRe)
self.DirichletBC0.zero(ARIm)
ARRer, ARRec, ARRev = fen.as_backend_type(
ARRe).mat().getValuesCSR()
ARImr, ARImc, ARImv = fen.as_backend_type(
ARIm).mat().getValuesCSR()
self.AR = (scsp.csr_matrix((ARRev, ARRec, ARRer),
shape = self.AL.shape)
+ 1.j * scsp.csr_matrix((ARImv, ARImc, ARImr),
shape = self.AL.shape))
self.A = self.AL - self.wavenumber2 * self.AM + Rfact * self.AR
if hasattr(self, "invA"): del self.invA
def assembleb(self):
"""Assemble RHS of linear system."""
raise Exception("Not implemented.")
if not hasattr(self, "b"):
if not hasattr(self, "bF"):
LFRe = fen.dot(self.fRe, self.v2) * fen.dx
LFIm = fen.dot(self.fIm, self.v2) * fen.dx
bFRe = fen.assemble(LFRe)
bFIm = fen.assemble(LFIm)
self.DirichletBCRe.apply(bFRe)
self.DirichletBCIm.apply(bFIm)
self.bF = np.array(bFRe.array()[:] + 1.j * bFIm.array()[:],
dtype = np.complex)
if not hasattr(self, "bN"):
LNRe = fen.dot(self.g1Re, self.v2) * self.ds(0)
LNIm = fen.dot(self.g1Im, self.v2) * self.ds(0)
bNRe = fen.assemble(LNRe)
bNIm = fen.assemble(LNIm)
self.DirichletBC0.apply(bNRe)
self.DirichletBC0.apply(bNIm)
self.bN = np.array(bNRe.array()[:] + 1.j * bNIm.array()[:],
dtype = np.complex)
if not hasattr(self, "bR"):
LRRe = fen.dot(self.g2Re, self.v2) * self.ds(1)
LRIm = fen.dot(self.g2Im, self.v2) * self.ds(1)
bRRe = fen.assemble(LRRe)
bRIm = fen.assemble(LRIm)
self.DirichletBC0.apply(bRRe)
self.DirichletBC0.apply(bRIm)
self.bR = np.array(bRRe.array()[:] + 1.j * bRIm.array()[:],
dtype = np.complex)
self.b = self.bF + self.bN + self.bR
def buildLU(self):
"""Build LU decomposition of A."""
raise Exception("Not implemented.")
if not hasattr(self, "A"): self.assembleA()
if not hasattr(self, "invA"):
warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning)
self.invA = spla.splu(self.A)
warnings.simplefilter("default", scsp.SparseEfficiencyWarning)
def solve(self, LU : bool = False)\
-> ("Fenics function", "Fenics function"):
"""
Find solution of linear system.
Args:
LU(optional): Whether to use a LU solver for the system. Defaults
to False.
Returns:
Real and imaginary parts of solution.
"""
raise Exception("Not implemented.")
if not hasattr(self, "u"):
if not hasattr(self, "A"): self.assembleA()
if not hasattr(self, "b"): self.assembleb()
if LU:
self.buildLU()
self.u = self.invA.solve(self.b)
else:
self.u = spla.spsolve(self.A, self.b)
self.__solved = True
return self.u

Event Timeline

c4science · Help