Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/main/ROMApproximantLagrange.py
4c3d184239eemaster

Last Change
Raw File

ROMApproximantLagrange.py
View Options

#!/usr/bin/python
import numpy as np
import utilities
from ROMApproximant import ROMApproximant
class ROMApproximantLagrange(ROMApproximant):
"""
ROM Lagrange interpolant computation for parametric problems.
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 k);
- 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.
ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
w(optional): Weight for norm computation. If None, set to
Re(np.mean(ks)). Defaults to None.
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 [].
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
k0: Default parameter.
ks: Array of snapshot parameters.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of snapshots current approximant relies upon.
S: Number of solution snapshots over which current approximant is
based upon.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
POD: Whether to compute POD of snapshots.
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.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
extraApproxParameters = ["S"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
ks : "Numpy 1D array" = np.array([0]), w : float = None,
approxParameters : dict = {}, plotSnap : list = []):
self.ks = ks
ROMApproximant.__init__(self, HFEngine, HSEngine,
np.mean(ks), w, approxParameters)
if w is None:
w = np.mean(self.ksRe)
self.w = w
self.plotSnap = plotSnap
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximant.parameterList(self)
+ ROMApproximantLagrange.extraApproxParameters)
@property
def ks(self):
"""Value of ks. Its assignment may reset snapshots."""
return self._ks
@ks.setter
def ks(self, ks):
if hasattr(self, 'ks'):
ksOld = self.ks
else:
ksOld = None
self._ks = np.array(ks)
if (ksOld is None or len(self.ks) != len(ksOld)
or not np.allclose(self.ks, ksOld, 1e-14)):
self.resetSamples()
@property
def ksRe(self):
"""Real part of ks."""
return np.real(self.ks)
@property
def ksIm(self):
"""Imaginary part of ks."""
return np.imag(self.ks)
@property
def zs(self):
"""Square of ks."""
return np.power(self.ks, 2.)
@property
def approxParameters(self):
"""Value of approximant parameters. Its assignment may change S."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = utilities.purgeDict(approxParams,
ROMApproximantLagrange.parameterList(self),
dictname = self.name() + ".approxParameters")
keyList = list(approxParameters.keys())
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantLagrange.extraApproxParameters,
True, True)
ROMApproximant.approxParameters.fset(self, approxParametersCopy)
if "S" in keyList:
self.S = approxParameters["S"]
elif hasattr(self, "S"):
self.S = self.S
else:
self.S = 2
@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. (ABSTRACT)"""
self.solSnapshots = None
self.RPOD = None
def computeSnapshots(self):
"""
Compute snapshots of solution map.
"""
if self.solSnapshots is None:
for j, k in enumerate(self.ks):
self.solveHF(k)
self.HSEngine.plot(self.uHF, name = "u({:.4f})".format(k),
what = self.plotSnap)
self.manageSnapshots(self.uHF, j)
def manageSnapshots(self, u:"numpy 1D array", 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):
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]
c4science · Help