Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/main/ROMApproximantLagrangeRB.py
4c3d184239eefeature_greedy_coarseni…

Last Change
Raw File

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
PI = np.pi
class ROMApproximantLagrangeRB(ROMApproximantLagrange):
"""
ROM RB 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;
- 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.
ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshots weights. Defaults to uniform = 1.
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 2;
- 'R': rank for Galerkin projection; defaults to S;
- 'nodesType': sampling node type; available values are:
- 'CHEBYSHEV': Chebyshev nodes;
- 'GAUSSLEGENDRE': Gauss-Legendre nodes;
- 'MANUAL': manual selection;
defaults to 'MANUAL'.
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.
ws: Array of snapshots weights.
w: Weight for norm computation.
approxRadius: Dummy radius of approximant (i.e. distance from k0 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;
- 'nodesType': sampling node type.
S: Number of solution snapshots over which current approximant is
based upon.
R: Rank for Galerkin projection.
nodesType: Sampling node type.
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 = ["R", "nodesType"]
nodesTypeParameters = ["CHEBYSHEV", "GAUSSLEGENDRE", "MANUAL"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
ks : "Numpy 1D array" = np.array([0]),
ws : "Numpy 1D array" = None, w : float = None,
approxParameters : dict = {}, plotSnap : list = []):
ROMApproximantLagrange.__init__(self, HFEngine, HSEngine, ks, w,
approxParameters, plotSnap)
if ws == None:
ws = np.ones(np.size(self.ks))
self.ws = ws
self.solLifting = self.HFEngine.liftDirichletData()
@property
def k0(self):
"""Dummy center of approximant (i.e. k0)."""
self.k0 = np.mean(self.ks)
return self._k0
@k0.setter
def k0(self, k0:bool):
self._k0 = k0
def resetSamples(self):
"""Reset samples."""
ROMApproximantLagrange.resetSamples(self)
self.projMat = None
def parameterList(self) -> list:
"""
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
if "nodesType" in keyList:
self.nodesType = approxParameters["nodesType"]
elif hasattr(self, "nodesType"):
self.nodesType = self.nodesType
else:
self.nodesType = "MANUAL"
@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.")
self.S = self.R
@property
def nodesType(self):
"""Value of nodesType."""
return self._nodesType
@nodesType.setter
def nodesType(self, nodesType):
if hasattr(self, "nodesType"): nodesTypeold = self.nodesType
else: nodesTypeold = -1
try:
nodesType = nodesType.upper().strip().replace(" ","")
if nodesType not in self.nodesTypeParameters: raise
except:
warnings.warn(("Prescribed nodesType not recognized. Overriding to"
" 'MANUAL'."))
nodesType = "MANUAL"
self._nodesType = nodesType
self._approxParameters["nodesType"] = self.nodesType
if nodesTypeold != self.nodesType:
self.resetSamples()
def manageSnapshots(self, u:"2-tuple of Fenics function", pos:int):
"""
Post-process snapshots of solution map.
Any specialization should include something like
self.solSnapshots[:, pos] = (np.array(u[0].vector()[:])
+ 1.j * np.array(u[1].vector()[:]))
Args:
u: 2-tuple containing real and imaginary parts of FE dofs of
snapshot.
pos: Derivative index.
"""
return ROMApproximantLagrange.manageSnapshots(self, u - self.solLifting, pos)
def setupApprox(self):
"""
Compute RB projection matrix.
See ``Householder triangularization of a quasimatrix'', L.Trefethen,
2008 for QR algorithm.
"""
if self.solSnapshots is None:
if self.nodesType == "MANUAL" and len(self.ks) != self.S:
warnings.warn(("Number of prescribed nodes different from S. "
"Overriding S to len(ks)"))
self.S = len(self.ks)
self.approxRadius = np.max(np.abs(self.k0 - self.ks))
if self.nodesType != "MANUAL":
if self.nodesType == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(self.S)
self.ws = weights * 2. / PI
elif self.nodesType == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(self.S)
self.ws = weights
self.ks = nodes * self.approxRadius + self.k0
self.ws = self.ws[:, None]
self.computeSnapshots()
if not self.checkComputedApprox():
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.assembleReducedSystem()
self.lastApproxParameters = copy(self.approxParameters)
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
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, k:complex) -> "Numpy 1D array":
"""
Solve RB linear system.
Args:
k: Target wavenumber.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
ARBk = self.ARBs[0][:self.R,:self.R]
bRBk = self.bRBs[0][:self.R]
for j in range(1, len(self.ARBs)):
ARBk = ARBk + np.power(k, j) * self.ARBs[j][:self.R, :self.R]
for j in range(1, len(self.bRBs)):
bRBk = bRBk + np.power(k, j) * self.bRBs[j][:self.R]
return self.projMat[:, :self.R].dot(np.linalg.solve(ARBk, bRBk))
def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Real and imaginary parts of approximant.
"""
self.setupApprox()
self.uApp = self.solLifting + self.solveReducedSystem(k)
return self.uApp
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.
"""
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) - self.k0 * (not centered)
except:
warnings.warn("Generalized eig algorithm did not converge.")
x = np.empty(A.shape[0])
x[:] = np.NaN
return x
c4science · Help