Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/hfengines/scipy/generic_problem_engine.py
0c9a0bbd9f4efeature_greedy_coarseni…

Last Change
Raw File

generic_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/>.
#
from abc import abstractmethod
import fenics as fen
import numpy as np
from scipy.sparse import csr_matrix
import scipy.sparse.linalg as scspla
from matplotlib import pyplot as plt
from rrompy.hfengines.base import ProblemEngineBase
from rrompy.utilities.base.types import Np1D, Np2D, strLst, Tuple, List
from rrompy.utilities.base import purgeList, getNewFilename
__all__ = ['GenericProblemEngine']
class GenericProblemEngine(ProblemEngineBase):
"""
Generic solver for parametric problems.
"""
def __init__(self):
super().__init__()
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def innerProduct(self, u:Np1D, v:Np1D) -> float:
"""Hilbert space scalar product."""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
return v.conj().T.dot(self.energyNormMatrix.dot(u))
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
normMat = fen.assemble(fen.dot(self.u, self.v) * fen.dx)
normMatFen = fen.as_backend_type(normMat).mat()
self.energyNormMatrix = csr_matrix(normMatFen.getValuesCSR()[::-1],
shape = normMat.size)
def solve(self, mu:complex, RHS : Np1D = None) -> Np1D:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
"""
if RHS is None: RHS = self.b(mu)
return scspla.spsolve(self.A(mu), RHS)
def residual(self, u:Np1D, mu:complex) -> Np1D:
"""
Find residual of linear system for given approximate solution.
Args:
u: numpy complex array with function dofs. If None, set to 0.
mu: parameter value.
"""
RHS = self.b(mu)
if u is None: return RHS
return RHS - self.A(mu).dot(u)
def plot(self, u:Np1D, name : str = "u", save : bool = False,
what : strLst = 'all', **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'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Whether to save plot(s). Defaults to False.
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()[:] = np.array(np.abs(u), dtype = float)
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()[:] = np.array(np.angle(u), dtype = float)
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()[:] = np.array(np.real(u), dtype = float)
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()[:] = np.array(np.imag(u), dtype = float)
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uIm, title = "Im({0})".format(name))
plt.colorbar(p)
if save:
plt.savefig(getNewFilename("fig", "eps"), format='eps', dpi=1000)
plt.show()
plt.close()
def plotmesh(self, name : str = "Mesh", save : bool = False, **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): Whether to save plot(s). Defaults to False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
plt.figure(**figspecs)
fen.plot(self.V.mesh())
if save:
plt.savefig(getNewFilename("msh", "eps"), format='eps', dpi=1000)
plt.show()
plt.close()
@abstractmethod
def A(self, mu:complex, der : int = 0) -> Np2D:
"""Assemble (derivative of) operator of linear system."""
pass
@abstractmethod
def b(self, mu:complex, der : int = 0) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
pass
@abstractmethod
def affineBlocksA(self, mu : complex = 0.) -> Tuple[List[Np2D], callable]:
"""Assemble affine blocks of operator of linear system."""
pass
@abstractmethod
def affineBlocksb(self, mu : complex = 0.) -> Tuple[List[Np1D], callable]:
"""Assemble affine blocks of RHS of linear system."""
pass
c4science · Help