Page MenuHomec4science
Files F63634751

ROMApproximantLagrangeRB.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 21, 12:29

ROMApproximantLagrangeRB.py
View Options

#!/usr/bin/python
from copy import copy
import warnings
import numpy as np
import scipy as sp
import utilities
from ROMApproximantLagrange import ROMApproximantLagrange
from RROMPyTypes import Np1D, ListAny, DictAny, HFEng, HSEng
class ROMApproximantLagrangeRB(ROMApproximantLagrange):
"""
ROM RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver. Should have members:
- Vdim: domain dimension;
- 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 mu to prescribed values;
- getLSBlocks: get algebraic system blocks;
- liftDirichletData: perform lifting of Dirichlet BC as numpy
complex 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.
mus(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshot weigths (unused). Defaults to
np.ones_like(self.mus).
w(optional): Weight for norm computation. If None, set to
Re(np.mean(self.mus)). 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 2;
- 'R': rank for Galerkin projection; defaults to S.
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.
w: Weight for norm computation.
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'R': rank for Galerkin projection.
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.
R: Rank for Galerkin projection.
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.
projMat: Projection matrix for RB system assembly.
autoNode: Type of nodes, if automatically generated. Otherwise None.
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.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix wrt mu.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS wrt mu.
"""
extraApproxParameters = ["R"]
def __init__(self, HFEngine:HFEng, HSEngine:HSEng,
mus : Np1D = np.array([0]), w : float = None,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
ROMApproximantLagrange.__init__(self, HFEngine = HFEngine,
HSEngine = HSEngine, mus = mus, w = w,
approxParameters = approxParameters,
plotSnap = plotSnap,
plotSnapSpecs = plotSnapSpecs)
self.solLifting = self.HFEngine.liftDirichletData()
def resetSamples(self):
"""Reset samples."""
ROMApproximantLagrange.resetSamples(self)
self.projMat = None
def parameterList(self) -> ListAny:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximantLagrange.parameterList(self)
+ ROMApproximantLagrangeRB.extraApproxParameters)
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = utilities.purgeDict(approxParams,
ROMApproximantLagrangeRB.parameterList(self),
dictname = self.name() + ".approxParameters")
keyList = list(approxParameters.keys())
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantLagrangeRB.extraApproxParameters,
True, True)
ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy)
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.S
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "S") and self.S < self.R:
warnings.warn("Prescribed S is too small. Updating S to R.",
stacklevel = 2)
self.S = self.R
def manageSnapshots(self, u:ListAny, pos:int) -> Np1D:
"""
Post-process snapshots of solution map.
Args:
u: Numpy 1D array of FE dofs of snapshot.
pos: Derivative index.
Returns:
Numpy 1D array with adjusted snapshot dofs.
"""
return ROMApproximantLagrange.manageSnapshots(self,
u - self.solLifting, 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.projMat is not None
and ROMApproximantLagrange.checkComputedApprox(self))
def setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
self.computeSnapshots()
if self.POD:
U, _, _ = np.linalg.svd(self.RPOD, full_matrices = False)
self.projMat = self.solSnapshots.dot(U[:, : self.R])
else:
self.projMat = self.solSnapshots[:, : self.R]
self.lastApproxParameters = copy(self.approxParameters)
self.assembleReducedSystem()
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
self.setupApprox()
projMatH = self.projMat.T.conjugate()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * max(len(self.As), len(self.bs))
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
if j < len(self.bs):
self.bRBs[j] = projMatH.dot(self.bs[j]
- self.As[j].dot(self.solLifting))
else:
self.bRBs[j] = - projMatH.dot(self.As[j].dot(self.solLifting))
for j in range(len(self.As), len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
ARBmu = self.ARBs[0][:self.R,:self.R]
bRBmu = self.bRBs[0][:self.R]
for j in range(1, len(self.ARBs)):
ARBmu = ARBmu + np.power(mu, j) * self.ARBs[j][:self.R, :self.R]
for j in range(1, len(self.bRBs)):
bRBmu = bRBmu + np.power(mu, j) * self.bRBs[j][:self.R]
return self.projMat[:, :self.R].dot(np.linalg.solve(ARBmu, bRBmu))
def evalApprox(self, mu:complex) -> Np1D:
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Numpy 1D array with approximant dofs.
"""
self.setupApprox()
self.uApp = self.solLifting + self.solveReducedSystem(mu)
return self.uApp
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
try:
return sp.linalg.eigvals(A, B)
except:
warnings.warn("Generalized eig algorithm did not converge.",
stacklevel = 2)
x = np.empty(A.shape[0])
x[:] = np.NaN
return x

Event Timeline

c4science · Help