Page MenuHomec4science

approximant_lagrange_pade_orthogonal.py
No OneTemporary

File Metadata

Created
Thu, Mar 28, 09:20

approximant_lagrange_pade_orthogonal.py

# 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
from rrompy.reduction_methods.base import checkRobustTolerance
from .approximant_lagrange_pade import ApproximantLagrangePade
from rrompy.utilities.base.types import Np1D, List
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangePadeOrthogonal']
class ApproximantLagrangePadeOrthogonal(ApproximantLagrangePade):
"""
ROM Lagrange Pade' interpolant computation for parametric problems with
rational-type interpolation basis functions.
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];
- 'E': coefficient of interpolant to be minimized; defaults to
min(S, M + 1);
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
defaults to None;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
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;
- 'E': coefficient of interpolant to be minimized;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust Pade' denominator management.
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.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
POD: Whether to compute POD of snapshots.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
@property
def sampler(self):
"""Value of sampler."""
return self._sampler
@sampler.setter
def sampler(self, sampler):
if 'generatePoints' not in dir(sampler):
raise Exception("Sampler type not recognized.")
if hasattr(self, '_sampler'):
samplerOld = self.sampler
self._sampler = sampler
self._approxParameters["sampler"] = self.sampler
if not 'samplerOld' in locals() or samplerOld != self.sampler:
if hasattr(self.sampler, "kind"):
if self.sampler.kind == "CHEBYSHEV":
self._val = np.polynomial.chebyshev.chebval
self._vander = np.polynomial.chebyshev.chebvander
self._fit = np.polynomial.chebyshev.chebfit
self._roots = np.polynomial.chebyshev.chebroots
self._domcoeff = lambda n: 2. ** max(0, n - 1)
else:#if self.sampler.kind in ["UNIFORM", "GAUSSLEGENDRE"]:
self._val = np.polynomial.legendre.legval
self._vander = np.polynomial.legendre.legvander
self._fit = np.polynomial.legendre.legfit
self._roots = np.polynomial.legendre.legroots
from scipy.special import binom
self._domcoeff = lambda n: binom(2 * n, n)
else:
self._val = np.polynomial.polynomial.polyval
self._vander = np.polynomial.polynomial.polyvander
self._fit = np.polynomial.polynomial.polyfit
self._roots = np.polynomial.polynomial.polyroots
self._domcoeff = lambda n: 1.
self.resetSamples()
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeRescaleParameter()
self.computeSnapshots()
if self.N > 0:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Starting computation of "
"denominator."))
TN = self._vander(self.radiusPade(self.mus), self.N)
while self.N > 0:
TN = TN[:, : self.N + 1]
if self.POD:
data = self.samplingEngine.RPOD
else:
data = self.samplingEngine.samples
RHSFull = np.empty((self.S, data.shape[0] * (self.N + 1)),
dtype = np.complex)
for j in range(self.S):
RHSFull[j, :] = np.kron(data[:, j], TN[j, :])
fitOut = self._fit(self.radiusPade(self.mus), RHSFull,
self.E, w = self.ws, full = True,
rcond = self.interpRcond)
if self.verbosity >= 5:
verbosityDepth("MAIN", ("Fitting {} samples with "
"degree {} through {}... "
"Conditioning of LS system: "
"{:.4e}.").format(
self.S, self.E, self._fit.__name__,
fitOut[1][2][0] / fitOut[1][2][-1]))
if fitOut[1][1] < self.E + 1:
Enew = fitOut[1][1] - 1
Nnew = min(self.N, Enew)
Mnew = min(self.M, Enew)
if Nnew == self.N:
strN = ""
else:
strN = "N from {} to {} and ".format(self.N, Nnew)
if Mnew == self.M:
strM = ""
else:
strM = "M from {} to {} and ".format(self.M, Mnew)
warn(("Polyfit is poorly conditioned.\nReducing {}{}E "
"from {} to {}.").format(strN, strM,
self.E, Enew))
newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
self.approxParameters = newParameters
continue
G = fitOut[0][-1, :].reshape((self.N + 1, data.shape[0]))
if self.POD:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square "
"root of gramian matrix."),
end = "")
_, ev, eV = np.linalg.svd(G, full_matrices = False)
ev = ev[::-1]
eV = eV[::-1, :].conj().T
else:
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.",
end = "")
G2 = self.HFEngine.innerProduct(G, G)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.",
inline = True)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving eigenvalue "
"problem for gramian "
"matrix."), end = "")
ev, eV = np.linalg.eigh(G2)
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
newParameters = checkRobustTolerance(ev, self.E,
self.robustTol)
if not newParameters:
break
self.approxParameters = newParameters
if self.N <= 0:
eV = np.ones((1, 1))
self.Q = eV[:, 0]
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.")
else:
self.Q = np.ones(1, dtype = np.complex)
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.")
self.lastApproxParameters = copy(self.approxParameters)
Qevaldiag = np.diag(self.getQVal(self.mus))
while self.M >= 0:
fitOut = self._fit(self.radiusPade(self.mus), Qevaldiag,
self.M, w = self.ws, full = True,
rcond = self.interpRcond)
if self.verbosity >= 5:
verbosityDepth("MAIN", ("Fitting {} samples with degree "
"{} through {}... Conditioning of "
"LS system: {:.4e}.").format(
self.S, self.M, self._fit.__name__,
fitOut[1][2][0] / fitOut[1][2][-1]))
if fitOut[1][1] == self.M + 1:
P = fitOut[0].T
break
warn(("Polyfit is poorly conditioned. Reducing M from {} to "
"{}. Exact snapshot interpolation not guaranteed.")\
.format(self.M, fitOut[1][1] - 1))
self.M = fitOut[1][1] - 1
self.P = np.atleast_2d(P)
if self.POD:
self.P = self.samplingEngine.RPOD.dot(self.P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.")
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def getPVal(self, mu:List[complex]):
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating numerator at mu = {}.".format(mu),
end = "")
try:
len(mu)
except:
mu = [mu]
p = self._val(self.radiusPade(mu), self.P.T)
if len(mu) == 1:
p = p.flatten()
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return p
def getQVal(self, mu:List[complex]):
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating denominator at mu = {}.".format(mu),
end = "")
q = self._val(self.radiusPade(mu), self.Q)
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
return self.HFEngine.rescalingInv(
self.scaleFactor * self._roots(self.Q)
+ self.HFEngine.rescaling(self.mu0))

Event Timeline