Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/main/ROMApproximantTaylor.py
7e2d435a8d1bmaster

Last Change
Raw File

ROMApproximantTaylor.py
View Options

#!/usr/bin/python
import warnings
import numpy as np
import utilities
from ROMApproximant import ROMApproximant
class ROMApproximantTaylor(ROMApproximant):
"""
ROM single-point approximant 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;
- setupDerivativeComputation: setup HF problem data to compute j-th
solution derivative at k0;
- solve: get HF solution as complex numpy vector.
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.
k0(optional): Default parameter. Defaults to 0.
w(optional): Weight for norm computation. If None, set to Re(k0).
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;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): Whether to plot derivatives of the Helmholtz
solution map. Defaults to False.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
k0: Default parameter.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type.
plotDer: Whether to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
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 = ["E", "Emax", "sampleType"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
k0 : complex = 0, w : float = None,
approxParameters : dict = {}, plotDer : bool = False):
self.plotDer = plotDer
ROMApproximant.__init__(self, HFEngine, HSEngine,
k0, w, approxParameters)
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximant.parameterList(self)
+ ROMApproximantTaylor.extraApproxParameters)
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change E and Emax.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = utilities.purgeDict(approxParams,
ROMApproximantTaylor.parameterList(self),
dictname = "approxParameters")
keyList = list(approxParameters.keys())
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantTaylor.extraApproxParameters,
True, True)
ROMApproximant.approxParameters.fset(self, approxParametersCopy)
if "E" in keyList:
self._E = approxParameters["E"]
self._approxParameters["E"] = self.E
if "Emax" in keyList:
self.Emax = approxParameters["Emax"]
else:
if not hasattr(self, "Emax"):
self.Emax = self.E
else:
self.Emax = self.Emax
else:
if "Emax" in keyList:
self._E = approxParameters["Emax"]
self._approxParameters["E"] = self.E
self.Emax = self.E
else:
if not (hasattr(self, "Emax") and hasattr(self, "E")):
raise Exception("At least one of E and Emax must be set.")
if "sampleType" in keyList:
self.sampleType = approxParameters["sampleType"]
elif hasattr(self, "sampleType"):
self.sampleType = self.sampleType
else:
self.sampleType = "KRYLOV"
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warnings.warn("Prescribed E is too large. Updating Emax to E.")
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
self._Emax = Emax
if hasattr(self, "E") and self.Emax < self.E:
warnings.warn("Prescribed Emax is too small. Updating Emax to E.")
self.Emax = self.E
else:
self._approxParameters["Emax"] = self.Emax
if EmaxOld >= self.Emax and self.solDerivatives is not None:
self.solDerivatives = self.solDerivatives[:, : self.Emax + 1]
if hasattr(self, "HArnoldi"):
self.HArnoldi = self.HArnoldi[: self.Emax + 1,
: self.Emax + 1]
self.RArnoldi = self.RArnoldi[: self.Emax + 1,
: self.Emax + 1]
else:
self.resetSamples()
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType
else: sampleTypeOld = -1
try:
sampleType = sampleType.upper().strip().replace(" ","")
if sampleType not in ["ARNOLDI", "KRYLOV"]: raise
if sampleType == "ARNOLDI" and not self.POD:
warnings.warn(("Prescribed sampleType not compatible with POD "
"value. Overriding to 'KRYLOV'."))
sampleType = "KRYLOV"
self._sampleType = sampleType
except:
warnings.warn(("Prescribed sampleType not recognized. Overriding "
"to 'KRYLOV'."))
self.sampleType = "KRYLOV"
self._approxParameters["sampleType"] = self.sampleType
if sampleTypeOld != self.sampleType:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
self.solDerivatives = None
if hasattr(self, "HArnoldi"):
del self.HArnoldi
if hasattr(self, "RArnoldi"):
del self.RArnoldi
def computeDerivatives(self):
"""
Compute derivatives of solution map starting from order 0.
"""
if self.solDerivatives is None:
up = np.zeros(self.HSEngine.V.dim())
upp = up
for j in range(self.Emax + 1):
if j == 0:
self.HFEngine.setProblemData(self.initialHFData, self.k0)
else:
self.HFEngine.setupDerivativeComputation(j, up, upp)
uj = self.HFEngine.solve(True)
if self.plotDer:
self.HSEngine.plot(uj, name = "u_{0}".format(j))
upp = up
up = self.manageDerivatives(uj, j)
def manageDerivatives(self, u:"numpy 1D array",
pos:int) -> ("numpy 1D array"):
"""
Store derivatives of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
Returns:
Adjusted solution derivative as numpy complex vector.
"""
if pos == 0:
self.solDerivatives = np.empty((u.shape[0], self.Emax + 1),
dtype = np.complex)
self.solDerivatives[:, pos] = u
if self.sampleType == "ARNOLDI":
if pos == 0:
self.HArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
dtype = np.complex)
self.RArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
dtype = np.complex)
self.HArnoldi[: pos, pos] = self.solDerivatives[:, pos].conj().dot(
self.energyNormMatrix.dot(self.solDerivatives[:, : pos])
).conj()
self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
- self.solDerivatives[:, : pos].dot(
self.HArnoldi[: pos, pos]))
self.HArnoldi[pos, pos] = (self.solDerivatives[:, pos].conj().dot(
self.energyNormMatrix.dot(self.solDerivatives[:, pos])))**.5
self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
/ self.HArnoldi[pos, pos])
if pos == 0:
self.RArnoldi[pos, pos] = self.HArnoldi[pos, pos]
else:
self.RArnoldi[:pos+1, pos] = self.HArnoldi[: pos+1, 1 : pos+1]\
.dot(self.RArnoldi[: pos, pos - 1])
return self.solDerivatives[:, pos]
c4science · Help