Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
c6af74908c2cfeature_greedy_coarseni…

Last Change
Raw File

generic_approximant_lagrange.py
View Options

#!/usr/bin/python
# 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 warnings
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng
__all__ = ['GenericApproximantLagrange']
class GenericApproximantLagrange(GenericApproximant):
"""
ROM Lagrange interpolant computation for parametric problems (ABSTRACT).
Args:
HFEngine: HF problem solver. Should have members:
- energyNormMatrix: assemble sparse matrix (in CSC format)
associated to weighted H10 norm;
- problemData: list of HF problem data (excluding mu);
- setProblemData: set HF problem data and k to prescribed values;
- getLSBlocks: get algebraic system blocks.
HSEngine: Hilbert space general purpose engine. Should have members:
- norm: compute Hilbert norm of expression represented as complex
numpy vector;
- plot: plot Hilbert expression represented as complex numpy
vector.
mus(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshot weigths (possibly unused). Defaults to
np.ones_like(self.mus).
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths (possibly unused).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of snapshots current approximant relies upon.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
autoNodeTypes: List of possible values for autoNode.
S: Number of solution snapshots over which current approximant is
based upon.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
autoNode: Type of nodes, if automatically generated. Otherwise None.
Possibly unused.
energyNormMatrix: Sparse matrix (in CSC format) assoociated to
w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
"""
autoNodeTypes = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, HFEngine:HFEng, HSEngine:HSEng,
mus : Np1D = np.array([0]), ws : Np1D = None,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
self.mus = mus
if ws is None:
ws = np.ones_like(self.mus)
self.ws = ws
super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine,
mu0 = self.mu0, approxParameters = approxParameters)
self.parameterList += ["S"]
self.plotSnap = plotSnap
self.plotSnapSpecs = plotSnapSpecs
self.resetSamples()
@property
def mu0(self):
"""Dummy center of approximant (i.e. mu0)."""
self.mu0 = np.mean(self.mus)
return self._mu0
@mu0.setter
def mu0(self, mu0:bool):
self._mu0 = mu0
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
musOld = self.mus if hasattr(self, 'mus') else None
self._mus = np.array(mus)
_, musCounts = np.unique(self._mus, return_counts = True)
if len(np.where(musCounts > 1)[0]) > 0:
raise Exception("Repeated sample points not allowed.")
if (musOld is None or len(self.mus) != len(musOld)
or not np.allclose(self.mus, musOld, 1e-14)):
self.resetSamples()
self.autoNode = None
@property
def ws(self):
"""Value of ws."""
return self._mus
@ws.setter
def ws(self, ws):
if hasattr(self, 'ws'):
wsOld = self.ws
else:
wsOld = None
ws = np.array(ws)
self._ws = np.abs(ws)
if not np.array_equal(self._ws, ws):
warnings.warn(("Negative weights not allowed. Overriding to "
"absolute values."), stacklevel = 2)
if (wsOld is None or len(self.ws) != len(wsOld)
or not np.allclose(self.ws, wsOld, 1e-14)):
self.autoNode = None
def checkNodesWeightsConsistency(self) -> bool:
"""Check if mus and ws have consistent lengths."""
return self.mus.shape == self.ws.shape
@property
def approxParameters(self):
"""Value of approximant parameters. Its assignment may change S."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = super().approxParameters.fset(self, approxParams)
keyList = list(approxParameters.keys())
if "S" in keyList:
self.S = approxParameters["S"]
elif hasattr(self, "S"):
self.S = self.S
else:
self.S = 2
return approxParameters
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
self.solSnapshots = None
self.RPOD = None
def setNodesWeights(self, xs : Np1D = np.array([0, 1]),
nodeTypes : str = "UNIFORM"):
"""
Assign sampling nodes and weights based on a quadrature formula.
Args:
xs: extremes of the sampling interval;
nodeTypes: type of nodes/weights; allowed values are in
self.autoNodeTypes.
"""
if len(xs) != 2:
raise Exception(("Value of xs not recognized. Should be numpy 1D "
"vector or list with length 2."))
try:
nodeTypes = nodeTypes.upper().strip().replace(" ","")
if nodeTypes not in self.autoNodeTypes: raise
except:
warnings.warn(("Prescribed nodeTypes not recognized. Overriding "
"to 'UNIFORM'."), stacklevel = 2)
nodeTypes = "UNIFORM"
if nodeTypes == "UNIFORM":
nodes = np.linspace(1., -1., self.S)
weights = np.ones((self.S, 1))
elif nodeTypes == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(self.S)
weights = weights[:, None] * 2. / np.pi
else: #if nodeTypes == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(self.S)
weights = weights[:, None]
mu0 = np.mean(xs)
r = (xs[0] - xs[1]) / 2.
self.mus = r * nodes + mu0
self.ws = np.abs(r) * weights
self.autoNode = nodeTypes
def computeSnapshots(self):
"""
Compute snapshots of solution map.
"""
if self.solSnapshots is None:
if len(self.mus) != self.S:
warnings.warn(("Number of prescribed nodes different from "
"S. Overriding S to len(mus)."), stacklevel = 2)
self.S = len(self.mus)
for j, mu in enumerate(self.mus):
self.solveHF(mu)
self.HSEngine.plot(self.uHF, name = "u({:.4f})".format(mu),
what = self.plotSnap, **self.plotSnapSpecs)
self.manageSnapshots(self.uHF, j)
def manageSnapshots(self, u:Np1D, pos:int):
"""
Store snapshots of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
"""
if pos == 0:
self.solSnapshots = np.empty((u.shape[0], self.S),
dtype = np.complex)
self.solSnapshots[:, pos] = u
if self.POD:
if pos == 0:
self.RPOD = np.eye(self.S, dtype = np.complex)
beta = 1
for j in range(2): #twice is enough!
nu = self.solSnapshots[:, pos].conj().dot(
self.energyNormMatrix.dot(self.solSnapshots[:, : pos])
).conj()
self.RPOD[: pos, pos] = self.RPOD[: pos, pos] + beta * nu
eta = (self.solSnapshots[:, pos]
- self.solSnapshots[:, : pos].dot(nu))
beta = eta.conj().dot(self.energyNormMatrix.dot(eta))**.5
self.solSnapshots[:, pos] = eta / beta
self.RPOD[pos, pos] = beta * self.RPOD[pos, pos]
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (self.solSnapshots is not None
and super().checkComputedApprox(self))
c4science · Help