Page MenuHomec4science
Files F61569432

approximant_taylor_rb.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 7, 12:36

approximant_taylor_rb.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
from copy import copy
import numpy as np
import scipy as sp
from .generic_approximant_taylor import GenericApproximantTaylor
from rrompy.sampling.base.pod_engine import PODEngine
from rrompy.utilities.base.types import Np1D, DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantTaylorRB']
class ApproximantTaylorRB(GenericApproximantTaylor):
"""
ROM single-point fast RB approximant computation for parametric problems
with polynomial dependence up to degree 2.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
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.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- '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.
POD: Whether to compute QR factorization of derivatives.
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'.
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.
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 mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing RB
coefficients of linear system matrix wrt mu.
bRBs: List of numpy vectors representing RB coefficients of linear
system RHS wrt mu.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity)
if self.verbosity >= 10:
verbosityDepth("INIT", "Computing affine blocks of system.")
self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0,
self.homogeneized)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done computing affine blocks.")
self._postInit()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@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):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["R"],
True, True, baselevel = 1)
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
else:
self.R = self.E + 1
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
GenericApproximantTaylor.POD.fset(self, POD)
if (hasattr(self, "sampleType") and self.sampleType == "ARNOLDI"
and not self.POD):
warn(("Arnoldi sampling implicitly forces POD-type derivative "
"management."))
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
GenericApproximantTaylor.sampleType.fset(self, sampleType)
if (hasattr(self, "POD") and not self.POD
and self.sampleType == "ARNOLDI"):
warn(("Arnoldi sampling implicitly forces POD-type derivative "
"management."))
@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:
warn("Prescribed E is too small. Updating E to R - 1.")
self.E = self.R - 1
def setupApprox(self):
"""Setup RB system."""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeDerivatives()
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing projection matrix.",
end = "")
if self.POD and not self.sampleType == "ARNOLDI":
self.PODEngine = PODEngine(self.HFEngine)
self.projMatQ, self.projMatR = self.PODEngine.QRHouseholder(
self.samplingEngine.samples)
if self.POD:
if self.sampleType == "ARNOLDI":
self.projMatR = self.samplingEngine.RArnoldi
self.projMatQ = self.samplingEngine.samples
U, _, _ = np.linalg.svd(self.projMatR[: self.E + 1,
: self.E + 1])
self.projMat = self.projMatQ[:, : self.E + 1].dot(U[:,
: self.R])
else:
self.projMat = self.samplingEngine.samples[:, : self.R + 1]
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
self.assembleReducedSystem()
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
if not self.checkComputedApprox():
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT", "Projecting affine terms of HF model.",
end = "")
projMatH = self.projMat.T.conj()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Assembling reduced model for mu = {}.".format(mu),
end = "")
ARBmu = self.thetaAs(mu, 0) * self.ARBs[0][:self.R,:self.R]
bRBmu = self.thetabs(mu, 0) * self.bRBs[0][:self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.ARBs)):
ARBmu += self.thetaAs(mu, j) * self.ARBs[j][:self.R, :self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.bRBs)):
bRBmu += self.thetabs(mu, j) * self.bRBs[j][:self.R]
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
if self.verbosity >= 5:
verbosityDepth("INIT",
"Solving reduced model for mu = {}.".format(mu),
end = "")
uRB = np.linalg.solve(ARBmu, bRBmu)
if self.verbosity >= 5:
verbosityDepth("DEL", " Done.", inline = True)
return uRB
def evalApproxReduced(self, mu:complex):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Computing RB solution at mu = {}.".format(mu))
self.uAppReduced = self.solveReducedSystem(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done computing RB solution.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.projMat[:, : self.R].dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
warn(("Impossible to compute poles in general affine parameter "
"dependence. Results subject to interpretation/rescaling, or "
"possibly completely wrong."))
self.setupApprox()
if len(self.ARBs) < 2:
return
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.
return self.HFEngine.rescalingInv(sp.linalg.eigvals(A, B)
+ self.HFEngine.rescaling(self.mu0))

Event Timeline

c4science · Help