Page MenuHomec4science
Files F62710254

ROMApproximant.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 15, 01:34

ROMApproximant.py
View Options

#!/usr/bin/python
from abc import abstractmethod
import numpy as np
import utilities
from copy import copy
import scipy.sparse.linalg as spla
class ROMApproximant:
"""
ROM 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;
- 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.
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.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
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.
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.
"""
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
k0 : complex = 0, w : float = None,
approxParameters : dict = {}):
self.HFEngine = HFEngine
self.HSEngine = HSEngine
self.initialHFData = self.HFEngine.problemData()
self.k0 = k0
if w is None:
w = np.real(self.k0)
self.w = w
self.approxParameters = approxParameters
self.energyNormMatrix = self.HFEngine.energyNormMatrix(self.w)
self.As, self.bs = self.HFEngine.getLSBlocks()
def name(self) -> str:
"""Approximant label."""
return self.__class__.__name__
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return ["POD"]
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = utilities.purgeDict(approxParams,
ROMApproximant.parameterList(self),
dictname = self.name() + ".approxParameters")
keyList = list(approxParameters.keys())
if "POD" in keyList:
self.POD = approxParameters["POD"]
elif hasattr(self, "POD"):
self.POD = self.POD
else:
self.POD = True
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD:bool):
if hasattr(self, "POD"): PODold = self.POD
else: PODold = -1
self._POD = POD
self._approxParameters["POD"] = self.POD
if PODold != self.POD:
self.resetSamples()
def solveHF(self, k : complex = None):
"""
Find high fidelity solution with original parameters and arbitrary
wavenumber.
Args:
k: Target wavenumber.
"""
if k is None: k = self.k0
if (not hasattr(self, "lastSolvedHF")
or not np.isclose(self.lastSolvedHF, k)):
A, b = self.buildLS(k)
self.uHF = spla.spsolve(A, b)
return self.uHF
def buildLS(self, k:complex, As : list = None, bs : list = None):
"""
Build linear system from polynomial blocks.
Args:
k: Parameter value.
As: Matrix blocks. If None, set to self.As. Defaults to None.
bs: RHS blocks. If None, set to self.bs. Defaults to None.
Returns:
Matrix and RHS of system.
"""
if As is None: As = self.As
if bs is None: bs = self.bs
A = copy(As[0])
b = copy(bs[0])
for j in range(1, len(As)):
A = A + np.power(k, j) * As[j]
for j in range(1, len(bs)):
b = b + np.power(k, j) * bs[j]
return A, b
@abstractmethod
def resetSamples(self):
"""Reset samples. (ABSTRACT)"""
pass
@abstractmethod
def setupApprox(self):
"""
Setup approximant. (ABSTRACT)
Any specialization should include something like
self.computeDerivatives()
if not self.checkComputedApprox():
...
self.lastApproxParameters = copy(self.approxParameters)
"""
pass
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (hasattr(self, "lastApproxParameters")
and self.approxParameters == self.lastApproxParameters)
@abstractmethod
def evalApprox(self, k:complex) -> ("numpy 1D array"):
"""
Evaluate approximant at arbitrary wavenumber. (ABSTRACT)
Any specialization should include something like
self.setupApprox()
self.uApp = ...
Args:
k: Target wavenumber.
Returns:
Approximant as numpy complex vector.
"""
pass
@abstractmethod
def getPoles(self, centered : bool = False) -> "numpy 1D array":
"""
Obtain approximant poles.
Args:
centered(optional): Whether to return pole positions relative to
approximation center. Defaults to False.
Returns:
Numpy complex vector of poles.
"""
pass
def HFNorm(self, k:complex, normType : "number or str" = None) -> float:
"""
Compute norm of HF solution at arbitrary wavenumber.
Args:
k: Target wavenumber.
normType(optional): Target norm identifier. If number, target norm
is weighted H^1 norm with given weight. If string, must be
recognizable by Fenics norm command. If None, set to w.
Defaults to None.
Returns:
Target norm of HFsolution.
"""
self.solveHF(k)
if normType is None: normType = self.w
return self.HSEngine.norm(self.uHF, normType)
def approxNorm(self, k:complex, normType : "number or str" = None)-> float:
"""
Compute norm of (translated) approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
normType(optional): Target norm identifier. If number, target norm
is weighted H^1 norm with given weight. If string, must be
recognizable by Fenics norm command. If None, set to w.
Defaults to None.
Returns:
Target norm of approximant.
"""
self.evalApprox(k)
if normType is None: normType = self.w
return self.HSEngine.norm(self.uApp, normType)
def approxError(self, k:complex, normType : "number or str" = None)\
-> float:
"""
Compute norm of approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
normType(optional): Target norm identifier. If number, target norm
is weighted H^1 norm with given weight. If string, must be
recognizable by Fenics norm command. If None, set to w.
Defaults to None.
Returns:
Target norm of (approximant - HFsolution).
"""
self.evalApprox(k)
self.solveHF(k)
self.evalApprox(k)
if normType is None: normType = self.w
return self.HSEngine.norm(self.uApp - self.uHF, normType)
def getHF(self, k:complex) -> "numpy 1D array":
"""
Get HF solution at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
HFsolution as numpy complex vector.
"""
self.solveHF(k)
return self.uHF
def getApp(self, k:complex) -> "numpy 1D array":
"""
Get approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Approximant as numpy complex vector.
"""
self.evalApprox(k)
return self.uApp
def plotHF(self, k:complex, name : str = "u", **figspecs):
"""
Do some nice plots of the HF solution at arbitrary wavenumber.
Args:
k: Target wavenumber.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
self.solveHF(k)
self.HSEngine.plot(self.uHF, name, **figspecs)
def plotApp(self, k:complex, name : str = "u", **figspecs):
"""
Do some nice plots of approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
self.evalApprox(k)
self.HSEngine.plot(self.uApp, name, **figspecs)
def plotErr(self, k:complex, name : str = "u", **figspecs):
"""
Do some nice plots of approximation error at arbitrary wavenumber.
Args:
k: Target wavenumber.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
self.evalApprox(k)
self.solveHF(k)
self.HSEngine.plot(self.uApp - self.uHF, name, **figspecs)

Event Timeline

c4science · Help