Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
88e97ff5f796feature_greedy_coarseni…

Last Change
Raw File

approximant_lagrange_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_lagrange import GenericApproximantLagrange
from rrompy.utilities.base.types import Np1D, DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangeRB']
class ApproximantLagrangeRB(GenericApproximantLagrange):
"""
ROM RB approximant computation for parametric problems.
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;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'R': rank for Galerkin projection; defaults to S.
Defaults to empty dict.
homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths (unused).
homogeneized: Whether to homogeneize Dirichlet BCs.
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 in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'R': rank for Galerkin projection.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
R: Rank for Galerkin projection.
POD: Whether to compute POD of snapshots.
samplingEngine: Sampling engine.
projMat: Projection matrix for RB system assembly.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt theta(mu).
bs: List of numpy vectors representing coefficients of linear system
RHS wrt theta(mu).
thetaAs: List of callables representing coefficients of linear system
matrix wrt mu.
thetabs: List of callables representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix wrt theta(mu).
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS wrt theta(mu).
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
self._addParametersToList(["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)
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
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:
warn("Prescribed S is too small. Updating S to R.")
self.S = self.R
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 super().checkComputedApprox())
def setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeSnapshots()
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing projection matrix.",
end = "")
if self.POD:
U, _, _ = np.linalg.svd(self.samplingEngine.RPOD,
full_matrices = False)
self.projMat = self.samplingEngine.samples.dot(U[:, : self.R])
else:
self.projMat = self.samplingEngine.samples[:, : self.R]
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
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.assembleReducedSystem()
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.assembleReducedSystem()
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.assembleReducedSystem()
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))
c4science · Help