Page MenuHomec4science
Files F91194705

fenics_engine_base.py
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 8, 20:53

fenics_engine_base.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/>.
#
from os import path, mkdir
import fenics as fen
import numpy as np
from matplotlib import pyplot as plt
from .hfengine_base import HFEngineBase
from rrompy.utilities.base.types import Np1D, strLst, FenFunc, Tuple, List
from rrompy.utilities.base import (purgeList, getNewFilename,
verbosityManager as vbMng)
from rrompy.solver.fenics import L2NormMatrix, fenplot, interp_project
from .boundary_conditions import BoundaryConditions
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['FenicsEngineBase']
class FenicsEngineBase(HFEngineBase):
"""Generic solver for parametric fenics problems."""
def __init__(self, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(verbosity = verbosity, timestamp = timestamp)
self.BCManager = BoundaryConditions("Dirichlet")
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.degree_threshold = degree_threshold
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
if not type(V).__name__ == 'FunctionSpace':
raise RROMPyException("V type not recognized.")
self.dsToBeSet = True
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
@property
def spacedim(self):
if hasattr(self, "_V"): return self.V.dim()
return super().spacedim
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
vbMng(self, "INIT", "Initializing boundary measures.", 20)
mesh = self.V.mesh()
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = mesh,
subdomain_data = boundary_markers)
self.dsToBeSet = False
vbMng(self, "DEL", "Done assembling boundary measures.", 20)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self.energyNormMatrix = L2NormMatrix(self.V)
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product.
"""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
self.energyNormDualMatrix = self.energyNormMatrix
def buildEnergyNormPartialDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
self.energyNormPartialDualMatrix = self.energyNormMatrix
def liftDirichletData(self) -> Np1D:
"""Lift Dirichlet datum."""
if not hasattr(self, "_liftedDirichletDatum"):
liftRe = interp_project(self.DirichletDatum[0], self.V)
liftIm = interp_project(self.DirichletDatum[1], self.V)
self._liftedDirichletDatum = (np.array(liftRe.vector())
+ 1.j * np.array(liftIm.vector()))
return self._liftedDirichletDatum
def reduceQuadratureDegree(self, fun:FenFunc, name:str):
"""Check whether to reduce compiler parameters to degree threshold."""
if not np.isinf(self.degree_threshold):
from ufl.algorithms.estimate_degrees import (
estimate_total_polynomial_degree as ETPD)
try:
deg = ETPD(fun)
except:
return False
if deg > self.degree_threshold:
vbMng(self, "MAIN",
("Reducing quadrature degree from {} to {} for "
"{}.").format(deg, self.degree_threshold, name), 15)
return True
return False
def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]):
"""
Iterate reduceQuadratureDegree over list and define reduce compiler
parameters.
"""
if funsNames is not None:
for fun, name in funsNames:
if self.reduceQuadratureDegree(fun, name):
return {"quadrature_degree" : self.degree_threshold}
return {}
def plot(self, u:Np1D, warping : List[callable] = None,
is_state : bool = False, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, show : bool = True, fenplotArgs : dict = {},
**figspecs) -> str:
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
warping(optional): Domain warping functions.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
fenplotArgs(optional): Optional arguments for fenplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename.
"""
if not is_state and not self.isCEye:
return super().plot(u, warping, name, save, what, saveFormat,
saveDPI, show, fenplotArgs, **figspecs)
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = (13. * len(what) / 4, 3)
subplotcode = 100 + len(what) * 10
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
uAb = fen.Function(self.V)
uAb.vector().set_local(np.abs(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fenplot(uAb, warping = warping, title = "|{0}|".format(name),
**fenplotArgs)
if self.V.mesh().geometric_dimension() > 1:
plt.colorbar(p)
if 'PHASE' in what:
uPh = fen.Function(self.V)
uPh.vector().set_local(np.angle(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fenplot(uPh, warping = warping,
title = "phase({0})".format(name), **fenplotArgs)
if self.V.mesh().geometric_dimension() > 1:
plt.colorbar(p)
if 'REAL' in what:
uRe = fen.Function(self.V)
uRe.vector().set_local(np.real(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fenplot(uRe, warping = warping, title = "Re({0})".format(name),
**fenplotArgs)
if self.V.mesh().geometric_dimension() > 1:
plt.colorbar(p)
if 'IMAG' in what:
uIm = fen.Function(self.V)
uIm.vector().set_local(np.imag(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fenplot(uIm, warping = warping, title = "Im({0})".format(name),
**fenplotArgs)
if self.V.mesh().geometric_dimension() > 1:
plt.colorbar(p)
if save is not None:
save = save.strip()
fileOut = getNewFilename("{}_fig_".format(save), saveFormat)
plt.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show:
plt.show()
plt.close()
return fileOut
def plotmesh(self, warping : List[callable] = None, name : str = "Mesh",
save : str = None, saveFormat : str = "eps",
saveDPI : int = 100, show : bool = True,
fenplotArgs : dict = {}, **figspecs) -> str:
"""
Do a nice plot of the mesh.
Args:
u: numpy complex array with function dofs.
warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
fenplotArgs(optional): Optional arguments for fenplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename.
"""
plt.figure(**figspecs)
fenplot(self.V.mesh(), warping = warping, **fenplotArgs)
if save is not None:
save = save.strip()
fileOut = getNewFilename("{}_msh_".format(save), saveFormat)
plt.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show:
plt.show()
plt.close()
return fileOut
def outParaview(self, u:Np1D, warping : List[callable] = None,
is_state : bool = False, name : str = "u",
filename : str = "out", time : float = 0.,
what : strLst = 'all', forceNewFile : bool = True,
folder : bool = False, filePW = None) -> str:
"""
Output complex-valued function with given dofs to ParaView file.
Args:
u: numpy complex array with function dofs.
warping(optional): Domain warping functions.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
time(optional): Timestamp.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
filePW(optional): Fenics File entity (for time series).
Returns:
Output filename.
"""
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if filePW is None:
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[0], self.V.mesh()))
if what == ['MESH']:
filePW << (self.V.mesh(), time)
if 'ABS' in what:
uAb = fen.Function(self.V, name = "{}_ABS".format(name))
uAb.vector().set_local(np.abs(u))
filePW << (uAb, time)
if 'PHASE' in what:
uPh = fen.Function(self.V, name = "{}_PHASE".format(name))
uPh.vector().set_local(np.angle(u))
filePW << (uPh, time)
if 'REAL' in what:
uRe = fen.Function(self.V, name = "{}_REAL".format(name))
uRe.vector().set_local(np.real(u))
filePW << (uRe, time)
if 'IMAG' in what:
uIm = fen.Function(self.V, name = "{}_IMAG".format(name))
uIm.vector().set_local(np.imag(u))
filePW << (uIm, time)
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[1], self.V.mesh()))
return filePW
def outParaviewTimeDomain(self, u:Np1D, omega:float,
warping : List[callable] = None,
is_state : bool = False,
timeFinal : float = None,
periodResolution : int = 20, name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folder : bool = False) -> str:
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
Args:
u: numpy complex array with function dofs.
omega: frequency.
warping(optional): Domain warping functions.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
Returns:
Output filename.
"""
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
omega = np.abs(omega)
t = 0.
dt = 2. * np.pi / omega / periodResolution
if timeFinal is None: timeFinal = 2. * np.pi / omega - dt
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[0], self.V.mesh()))
for j in range(int(np.ceil(timeFinal / dt)) + 1):
ut = fen.Function(self.V, name = name)
ut.vector().set_local(np.real(u) * np.cos(omega * t)
+ np.imag(u) * np.sin(omega * t))
filePW << (ut, t)
t += dt
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[1], self.V.mesh()))
return filePW

Event Timeline

c4science · Help