Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py
51cd9a0a331afeature_greedy_coarseni…

Last Change
Raw File

approximant_lagrange_greedy_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/>.
#
import numpy as np
from copy import copy
from .generic_approximant_lagrange_greedy import (
GenericApproximantLagrangeGreedy)
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB
from rrompy.utilities.base.types import DictAny, HFEng, List
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangeRBGreedy']
class ApproximantLagrangeRBGreedy(GenericApproximantLagrangeGreedy,
ApproximantLagrangeRB):
"""
ROM greedy 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;
- 'muBounds': list of bounds for parameter values; defaults to
[[0], [1]];
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement; defaults to 0.2;
- 'nTrainingPoints': number of training points; defaults to
maxIter / refinementRatio;
- 'nTestPoints': number of starting test points; defaults to 1;
- 'trainingSetGenerator': training sample points generator;
defaults to uniform sampler within muBounds.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
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;
- 'muBounds': list of bounds for parameter values;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement;
- 'nTrainingPoints': number of training points;
- 'nTestPoints': number of starting test points;
- 'trainingSetGenerator': training sample points generator;
- 'robustTol': tolerance for robust Pade' denominator management.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
POD: whether to compute POD of snapshots.
muBounds: list of bounds for parameter values.
greedyTol: uniform error tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTrainingPoints: number of training points.
nTestPoints: number of starting test points.
trainingSetGenerator: training sample points generator.
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()
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()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
self._S = S
@property
def R(self):
"""Value of R."""
return self._S
@R.setter
def R(self, R):
warn(("R is used just to simplify inheritance, and its value cannot "
"be changed from that of S."))
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
self.resbb = None
self.resAb = None
self.resAA = None
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 errorEstimator(self, mus:List[np.complex],
homogeneized : bool = None) -> List[np.complex]:
"""
Standard residual-based error estimator. Unreliable for unstable
problems (inf-sup constant is missing).
"""
if homogeneized is None:
homogeneized = self.homogeneized
nmus = len(mus)
err = [0.] * nmus
self.assembleReducedSystem()
self.assembleReducedResidualBlocks(homogeneized)
assert homogeneized == self.homogeneized
nAs = len(self.As)
nbs = len(self.bs)
for j in range(nmus):
mu = mus[j]
uAppReduced = self.getAppReduced(mu)
prodbb = 0.
prodAb = 0.
prodAA = 0.
for i1 in range(nbs):
rhobi1 = self.thetabs(mu, i1)
for i2 in range(nbs):
rhobi2 = self.thetabs(mu, i2).conj()
prodbb += (rhobi1 * rhobi2
* self.resbb[homogeneized][i1][i2])
for i1 in range(nAs):
rhoAi1 = self.thetaAs(mu, i1)
for i2 in range(nbs):
rhobi2 = self.thetabs(mu, i2).conj()
prodAb += (rhoAi1 * rhobi2
* self.resAb[homogeneized][i1][i2])
for i1 in range(nAs):
rhoAi1 = self.thetaAs(mu, i1)
for i2 in range(nAs):
rhoAi2 = self.thetaAs(mu, i2).conj()
prodAA += (rhoAi1 * rhoAi2
* self.resAA[homogeneized][i1][i2])
err[j] = np.abs(((uAppReduced.T.conj().dot(prodAA.dot(uAppReduced))
- 2. * np.real(prodAb.dot(uAppReduced))) / prodbb
+ 1.)[0]) ** .5
return err
def setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.greedy()
self.S = len(self.mus)
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing projection matrix.",
end = "")
self.projMat = self.samplingEngine.samples
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 assembleReducedResidualBlocks(self, homogeneized:bool):
"""Build affine blocks of RB linear system through projections."""
self.initMassMatrixInverse()
self.assembleReducedSystem()
computeResbb = (self.resbb is None
or homogeneized not in self.resbb.keys())
computeResAb = (self.resAb is None
or homogeneized not in self.resAb.keys()
or self.resAb[homogeneized][0][0].shape[0] != self.S)
computeResAA = (self.resAA is None
or homogeneized not in self.resAA.keys()
or self.resAA[homogeneized][0][0].shape[0] != self.S)
if computeResbb or computeResAb or computeResAA:
if self.verbosity >= 7:
verbosityDepth("INIT", "Projecting affine terms of residual.",
end = "")
nAs = len(self.As)
nbs = len(self.bs)
if computeResbb:
if self.resbb is None: self.resbb = {}
self.resbb[homogeneized] = []
for i in range(nbs):
resbbi = []
Mbi = self.massInv.solve(self.bs[i])
for j in range(i):
Mbj = self.massInv.solve(self.bs[j])
resbbi += [self.HFEngine.innerProduct(Mbi, Mbj)]
resbbi += [self.HFEngine.innerProduct(Mbi, Mbi)]
self.resbb[homogeneized] += [resbbi]
for i in range(nbs):
for j in range(i + 1, nbs):
self.resbb[homogeneized][i] += [
self.resbb[homogeneized][j][i].conj()]
if computeResAb:
if self.resAb is None: self.resAb = {}
if homogeneized not in self.resAb.keys():
self.resAb[homogeneized] = []
for i in range(nAs):
resAbi = []
MAi = self.massInv.solve(self.As[i].dot(self.projMat))
for j in range(nbs):
Mbj = self.massInv.solve(self.bs[j])
resAbi += [self.HFEngine.innerProduct(MAi, Mbj)]
self.resAb[homogeneized] += [resAbi]
else:
Sold = self.resAb[homogeneized][0][0].shape[0]
for i in range(nAs):
MAi = self.massInv.solve(self.As[i].dot(
self.projMat[:, Sold :]))
for j in range(nbs):
Mbj = self.massInv.solve(self.bs[j])
self.resAb[homogeneized][i][j] = np.append(
self.resAb[homogeneized][i][j],
self.HFEngine.innerProduct(MAi, Mbj))
if computeResAA:
if self.resAA is None: self.resAA = {}
if homogeneized not in self.resAA.keys():
self.resAA[homogeneized] = []
for i in range(nAs):
resAAi = []
MAi = self.massInv.solve(self.As[i].dot(self.projMat))
for j in range(i):
MAj = self.massInv.solve(self.As[j].dot(
self.projMat))
resAAi += [self.HFEngine.innerProduct(MAi, MAj)]
resAAi += [self.HFEngine.innerProduct(MAi, MAi)]
self.resAA[homogeneized] += [resAAi]
for i in range(nAs):
for j in range(i + 1, nAs):
self.resAA[homogeneized][i] += [
self.resAA[homogeneized][j][i].conj()]
else:
Sold = self.resAA[homogeneized][0][0].shape[0]
for i in range(nAs):
MAi = self.massInv.solve(self.As[i].dot(self.projMat))
for j in range(i):
MAj = self.massInv.solve(self.As[j].dot(
self.projMat))
resAAcross1 = self.HFEngine.innerProduct(
MAi[:, Sold :], MAj[:, : Sold])
resAAcross2 = self.HFEngine.innerProduct(
MAi[:, : Sold], MAj[:, Sold :])
resAAdiag = self.HFEngine.innerProduct(
MAi[:, Sold :], MAj[:, Sold :])
self.resAA[homogeneized][i][j] = np.block(
[[self.resAA[homogeneized][i][j], resAAcross1],
[ resAAcross2, resAAdiag]])
resAAcross = self.HFEngine.innerProduct(MAi[:, Sold :],
MAi[:, : Sold])
resAAdiag = self.HFEngine.innerProduct(MAi[:, Sold :],
MAi[:, Sold :])
self.resAA[homogeneized][i][i] = np.block(
[[self.resAA[homogeneized][i][i], resAAcross],
[ resAAcross.T.conj(), resAAdiag]])
for i in range(nAs):
for j in range(i + 1, nAs):
self.resAA[homogeneized][i][j] = (
self.resAA[homogeneized][j][i].conj())
if self.verbosity >= 7:
verbosityDepth("DEL", ("Done setting up affine "
"decomposition of residual."))
c4science · Help