Page MenuHomec4science
Files F60268763

problem_engine_base.py
No OneTemporary
Actions

File Metadata

Created
Sun, Apr 28, 19:12

problem_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 abc import abstractmethod
from os import path, mkdir
import fenics as fen
import numpy as np
from scipy.sparse import csr_matrix
from matplotlib import pyplot as plt
from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, strLst, FenFunc,
Tuple, List)
from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
from .boundary_conditions import BoundaryConditions
from .matrix_engine_base import MatrixEngineBase
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['ProblemEngineBase']
class ProblemEngineBase(MatrixEngineBase):
"""
Generic solver for parametric problems.
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
"""
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.bsmu = np.nan
self.liftDirichletDatamu = np.nan
self.mu0BC = np.nan
self.degree_threshold = degree_threshold
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
self.resetAs()
self.resetbs()
if not type(V).__name__ == 'FunctionSpace':
raise RROMPyException("V type not recognized.")
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
def spacedim(self):
return self.V.dim()
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D:
"""Hilbert space scalar product."""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
return super().innerProduct(u, v, onlyDiag)
def buildEnergyNormForm(self): # L2
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling energy matrix.",
timestamp = self.timestamp)
normMatFen = fen.assemble(fen.dot(self.u, self.v) * fen.dx)
normMat = fen.as_backend_type(normMatFen).mat()
self.energyNormMatrix = csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling energy matrix.",
timestamp = self.timestamp)
def setDirichletDatum(self, mu:complex):
"""Set Dirichlet datum if parametric."""
if hasattr(self, "liftedDirichletDatum"):
self.liftDirichletDatamu = mu
def liftDirichletData(self, mu:complex) -> Np1D:
"""Lift Dirichlet datum."""
self.setDirichletDatum(mu)
if not np.isclose(self.liftDirichletDatamu, mu):
try:
liftRe = fen.interpolate(self.DirichletDatum[0], self.V)
except:
liftRe = fen.project(self.DirichletDatum[0], self.V)
try:
liftIm = fen.interpolate(self.DirichletDatum[1], self.V)
except:
liftIm = fen.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:
if self.verbosity >= 15:
verbosityDepth("MAIN", ("Reducing quadrature degree from "
"{} to {} for {}.").format(
deg,
self.degree_threshold,
name),
timestamp = self.timestamp)
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 {}
@abstractmethod
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
if self.As[der] is None:
self.As[der] = 0.
return self.As[der]
@abstractmethod
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if self.bs[homogeneized][der] is None:
self.bs[homogeneized][der] = 0.
return self.bs[homogeneized][der]
def plot(self, u:Np1D, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, show : bool = True, **figspecs):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
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.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
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 = fen.plot(uAb, title = "|{0}|".format(name))
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 = fen.plot(uPh, title = "phase({0})".format(name))
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 = fen.plot(uRe, title = "Re({0})".format(name))
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 = fen.plot(uIm, title = "Im({0})".format(name))
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
def plotmesh(self, name : str = "Mesh", save : str = None,
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, **figspecs):
"""
Do a nice plot of the mesh.
Args:
u: numpy complex array with function dofs.
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.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
plt.figure(**figspecs)
fen.plot(self.V.mesh())
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
def outParaview(self, u:Np1D, name : str = "u", filename : str = "out",
time : float = 0., what : strLst = 'all',
forceNewFile : bool = True, folder : bool = False,
filePW = None):
"""
Output complex-valued function with given dofs to ParaView file.
Args:
u: numpy complex array with function dofs.
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).
"""
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 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)
return filePW
def outParaviewTimeDomain(self, u:Np1D, omega:float,
timeFinal : float = None,
periodResolution : int = 20, name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folder : bool = False):
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
Args:
u: numpy complex array with function dofs.
omega: frequency.
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.
"""
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
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
return filePW

Event Timeline

c4science · Help