Page MenuHomec4science
Files F65533039

ROMApproximantTaylorRB.py
No OneTemporary
Actions

File Metadata

Created
Tue, Jun 4, 11:42

ROMApproximantTaylorRB.py
View Options

#!/usr/bin/python
from copy import copy
import warnings
import numpy as np
import scipy as sp
import utilities
from ROMApproximantTaylor import ROMApproximantTaylor
class ROMApproximantTaylorRB(ROMApproximantTaylor):
"""
ROM single-point fast RB approximant computation for parametric problems
with polynomial dependence up to degree 2.
Args:
HFEngine: HF problem solver. Should have members:
- energyNormMatrix: sparse matrix (in CSC format) associated to
w-weighted H10 norm;
- problemData: list of HF problem data (excluding k);
- setProblemData: set HF problem data and k to prescribed values;
- getLSBlocks: get blocks of HF linear system as sparse matrices in
CSC format;
- liftDirichletData: perform lifting of Dirichlet BC as numpy
complex vector;
- 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;
- 'R': rank for Galerkin projection; defaults to E + 1;
- '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;
- 'R': rank for Galerkin projection;
- '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.
R: Rank for Galerkin projection.
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, i.e. 'KRYLOV'.
plotDer: Whether to plot derivatives of the Helmholtz solution map.
energyNormMatrix: Sparse matrix (in CSC format) associated 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.
solLifting: Numpy complex vector with lifting of real part of Dirichlet
boundary datum.
projMat: Numpy matrix representing projection onto RB space.
projMat: Numpy matrix representing projection onto RB space.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt k.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt k.
ARBs: List of sparse matrices (in CSC format) representing RB
coefficients of linear system matrix wrt k.
bRBs: List of numpy vectors representing RB coefficients of linear
system RHS wrt k.
"""
extraApproxParameters = ["R"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
k0 : complex = 0, w : float = None,
approxParameters : dict = {}, plotDer : bool = False):
ROMApproximantTaylor.__init__(self, HFEngine, HSEngine, k0, w,
approxParameters, plotDer)
self.solLifting = self.HFEngine.liftDirichletData()
def resetSamples(self):
"""Reset samples."""
ROMApproximantTaylor.resetSamples(self)
self.projMat = None
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximantTaylor.parameterList(self)
+ ROMApproximantTaylorRB.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,
ROMApproximantTaylorRB.parameterList(self),
dictname = "approxParameters")
keyList = list(approxParameters.keys())
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantTaylorRB.extraApproxParameters,
True, True)
ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy)
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.E + 1
@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, "E") and self.E + 1 < self.R:
warnings.warn("Prescribed E is too small. Updating E to R - 1.")
self.E = self.R - 1
def manageDerivatives(self, u:"2-tuple of Fenics function", pos:int)\
-> ("Fenics function", "Fenics function"):
"""
Store derivatives of solution map. Subtract lifting for order 0.
Args:
u: 2-tuple containing real and imaginary parts of FE dofs of
derivative.
pos: Derivative index.
Returns:
Real and imaginary parts of derivative (possibly adjusted).
"""
if pos == 0:
self.As, self.bs = self.HFEngine.getLSBlocks()
u = u - self.solLifting
return ROMApproximantTaylor.manageDerivatives(self, u, pos)
def setupApprox(self):
"""
Setup RB system. For usage of correlation matrix in POD see Section
6.3.1 in 'Reduced Basis Methods for PDEs. An introduction', A.
Quarteroni, A. Manzoni, F. Negri, Springer, 2016.
"""
need2Setup = (self.solDerivatives is None) or (self.projMat is None)
self.computeDerivatives()
if need2Setup:
if self.POD and not self.sampleType == "ARNOLDI":
A = copy(self.solDerivatives)
M = self.energyNormMatrix
V = np.zeros_like(A, dtype = np.complex)
Q = np.zeros_like(A, dtype = np.complex)
R = np.zeros((A.shape[1], A.shape[1]), dtype = np.complex)
for k in range(A.shape[1]):
Q[:, k] = (np.random.randn(Q.shape[0])
+ 1.j * np.random.randn(Q.shape[0]))
if k > 0:
for times in range(2):
Q[:, k] = Q[:, k] - Q[:, :k].dot(
(Q[:, k].conj().dot(M.dot(Q[:, :k]))).conj())
Q[:, k] = Q[:, k]/(Q[:, k].conj().dot(M.dot(Q[:, k])))**.5
R[k, k] = np.abs(A[:, k].conj().dot(M.dot(A[:, k]))) ** .5
alpha = Q[:, k].conj().dot(M.dot(A[:, k]))
if np.isclose(np.abs(alpha), 0.): s = 1.
else: s = - alpha / np.abs(alpha)
Q[:, k] = s * Q[:, k]
v = R[k, k] * Q[:, k] - A[:, k]
for times in range(2):
v = v - Q[:, :k].dot((v.conj().dot(M.dot(Q[:, :k]))
).conj())
sigma = np.abs(v.conj().dot(M.dot(v))) ** .5
if np.isclose(sigma, 0.): v = Q[:, k]
else: v = v / sigma
V[:, k] = v
J = np.arange(k + 1, A.shape[1])
vtQ = v.conj().dot(M.dot(A[:, J]))
A[:, J] = A[:, J] - 2 * np.outer(v, vtQ)
R[k, J] = Q[:, k].conj().dot(M.dot(A[:, J]))
A[:, J] = A[:, J] - np.outer(Q[:, k], R[k, J])
for k in range(A.shape[1] - 1, -1, -1):
v = V[:, k]
J = np.arange(k, A.shape[1])
vtQ = v.conj().dot(M.dot(Q[:, J]))
Q[:, J] = Q[:, J] - 2 * np.outer(v, vtQ)
self.projMatQ = Q
self.projMatR = R
if self.POD and not self.sampleType == "ARNOLDI":
U, _, _ = np.linalg.svd(self.projMatR[: self.R, : self.R])
self.projMat = self.projMatQ[:, : self.R].dot(U)
else:
self.projMat = self.solDerivatives[:, : 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

Event Timeline

c4science · Help