Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/main/ROMApproximantLagrangePade.py
7e2d435a8d1bmaster

Last Change
Raw File

ROMApproximantLagrangePade.py
View Options

#!/usr/bin/python
from __future__ import print_function
from copy import copy
import warnings
import numpy as np
import utilities
from ROMApproximantLagrange import ROMApproximantLagrange
PI = np.pi
class ROMApproximantLagrangePade(ROMApproximantLagrange):
"""
ROM Lagrange Pade' interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver. Should have members:
- energyNormMatrix: assemble sparse matrix (in CSC format)
associated to weighted H10 norm;
- problemData: list of HF problem data (excluding k);
- setProblemData: set HF problem data and k to prescribed values;
- 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.
ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshots weights. Defaults to uniform = 1.
w(optional): Weight for norm computation. If None, set to
Re(np.mean(ks)). 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;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0.
- 'polyBasis': label of polynomial basis for LS problem; available
values are:
- 'CHEBYSHEV': Chebyshev nodes and weights;
- 'GAUSSLEGENDRE': Gauss-Legendre nodes and weights;
defaults to 'CHEBYSHEV'.
Defaults to empty dict.
plotSnapshots(optional): Whether to plot snapshots of the Helmholtz
solution map. Defaults to False.
verboseRobust(optional): Verbosity level for robustness-related
parameter modifications. Defaults to 1.
cleanupParameters(optional): Parameter values for pole cleanup.
Available fields are:
- 'boolCondition'(bool function handle): if evaluation on pole
returns False, pole is removed; defaults to always True;
- 'residueCheck'(bool): whether to apply residue check; defaults to
False;
- 'residueNPoints'(int): number of sample points to be used in
residue estimation; defaults to 2;
- 'residueRadius'(float): sample radius to be used in residue
estimation by Cauchy formula; defaults to 1e-5;
- 'residueTol'(float): target tolerance to be used in residue
estimation; defaults to 1e-4.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
k0: Default parameter.
ks: Array of snapshot parameters.
ws: Array of snapshots weights.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- 'polyBasis': label of polynomial basis for LS problem.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
S: Number of solution snapshots over which current approximant is
based upon.
polyBasis: Polynomial basis for LS problem.
plotSnapshots: Whether to plot snapshots of the Helmholtz solution map.
POD: Whether to compute POD of snapshots.
verboseRobust: Verbosity level for robustness-related parameter
modifications.
cleanupParameters; Parameter values for pole cleanup. Available fields
are:
- 'boolCondition': if evaluation on pole returns False, pole is
removed;
- 'residueCheck': whether to apply residue check;
- 'residueNPoints': number of sample points to be used in residue
estimation;
- 'residueRadius': sample radius to be used in residue estimation
by Cauchy formula;
- 'residueTol': target tolerance to be used in residue estimation.
energyNormMatrix: Sparse matrix (in CSC format) assoociated 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.
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.
"""
extraApproxParameters = ["M", "N", "polyBasis"]
polyBasisParameters = ["CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
ks : "Numpy 1D array" = np.array([0]),
ws : "Numpy 1D array" = None, w : float = None,
approxParameters : dict = {}, plotSnapshots : bool = False,
verboseRobust : int = 1, cleanupParameters : dict = {}):
ROMApproximantLagrange.__init__(self, HFEngine, HSEngine, ks, w,
approxParameters, plotSnapshots)
if ws == None:
ws = np.ones(np.size(self.ks))
self.ws = ws
self.verboseRobust = verboseRobust
self.cleanupParameters = cleanupParameters
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximantLagrange.parameterList(self)
+ ROMApproximantLagrangePade.extraApproxParameters)
@property
def k0(self):
"""Dummy center of approximant (i.e. k0)."""
self.k0 = np.mean(self.ks)
return self._k0
@k0.setter
def k0(self, k0:bool):
self._k0 = k0
@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,
ROMApproximantLagrangePade.parameterList(self),
dictname = "approxParameters")
keyList = list(approxParameters.keys())
if "M" in keyList:
self.M = 0 #to avoid warnings if M and S are changed simultaneously
if "N" in keyList:
self.N = 0 #to avoid warnings if N and S are changed simultaneously
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantLagrangePade.extraApproxParameters,
True, True)
ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy)
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 0
if "polyBasis" in keyList:
self.polyBasis = approxParameters["polyBasis"]
elif hasattr(self, "polyBasis"):
self.polyBasis = self.polyBasis
else:
self.polyBasis = "CHEBYSHEV"
@property
def M(self):
"""Value of M. Its assignment may change S."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if hasattr(self, "S") and self.S < self.M + 1:
warnings.warn("Prescribed S is too small. Updating S to M + 1.")
self.S = self.M + 1
@property
def N(self):
"""Value of N. Its assignment may change S."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
if hasattr(self, "S") and self.S < self.N + 1:
warnings.warn("Prescribed S is too small. Updating S to N + 1.")
self.S = self.N + 1
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
vals, label = [0] * 2, {0:"M", 1:"N"}
if hasattr(self, "M"): vals[0] = self.M
if hasattr(self, "N"): vals[1] = self.N
idxmax = np.argmax(vals)
if vals[idxmax] + 1 > S:
warnings.warn("Prescribed S is too small. Updating S to {} + 1."\
.format(label[idxmax]))
self.Emax = vals[idxmax] + 1
else:
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
@property
def polyBasis(self):
"""Value of polyBasis."""
return self._polyBasis
@polyBasis.setter
def polyBasis(self, polyBasis):
if hasattr(self, "polyBasis"): polyBasisold = self.polyBasis
else: polyBasisold = -1
try:
polyBasis = polyBasis.upper().strip().replace(" ","")
if polyBasis not in self.polyBasisParameters: raise
except:
warnings.warn(("Prescribed polyBasis not recognized. Overriding to"
" 'CHEBYSHEV'."))
polyBasis = "CHEBYSHEV"
self._polyBasis = polyBasis
self._approxParameters["polyBasis"] = self.polyBasis
if polyBasisold != self.polyBasis:
self.resetSamples()
@property
def cleanupParameters(self):
"""Value of cleanupParameters."""
return self._cleanupParameters
@cleanupParameters.setter
def cleanupParameters(self, cleanupParameters):
allowedCleanupKeys = ['boolCondition','residueCheck','residueNPoints',
'residueRadius','residueTol']
cleanupKeys = cleanupParameters.keys()
if 'boolCondition' not in cleanupKeys:
cleanupParameters['boolCondition'] = lambda x : True
if 'residueCheck' not in cleanupKeys:
cleanupParameters['residueCheck'] = False
if 'residueNPoints' not in cleanupKeys:
cleanupParameters['residueNPoints'] = 2
if 'residueRadius' not in cleanupKeys:
cleanupParameters['residueRadius'] = 1e-5
if 'residueTol' not in cleanupKeys:
cleanupParameters['residueTol'] = 1e-4
cleanupParameters['boolCondition'] = np.vectorize(
cleanupParameters['boolCondition'])
self._cleanupParameters = {key : cleanupParameters[key]
for key in allowedCleanupKeys}
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
S0 = self.S
M1 = self.M + 1
N1 = self.N + 1
if self.solSnapshots is None:
self.approxRadius = np.max(np.abs(self.k0 - self.ks))
if self.polyBasis == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(S0)
self.ws = weights * 2. / PI
elif self.polyBasis == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(S0)
self.ws = weights
self.ks = nodes * self.approxRadius + self.k0
self.ws = self.ws[:, None]
self.computeSnapshots()
if not self.checkComputedApprox():
Phi = np.zeros((S0, max(M1, N1)), dtype = np.complex)
Phi[:, 0] = np.ones((S0,)) / 2**.5
for j in range(1, max(M1, N1)):
c = np.zeros((j + 1,))
c[-1] = 1.
if self.polyBasis == "CHEBYSHEV":
Phi[:, j] = np.polynomial.chebyshev.chebval(
self.radiusPade(self.ks), c)
elif self.polyBasis == "GAUSSLEGENDRE":
Phi[:, j] = (j + .5) ** .5 * np.polynomial.legendre.legval(
self.radiusPade(self.ks), c)
if self.POD:
II = np.array(np.arange(0, S0**3, S0)
+ np.kron(np.arange(S0), np.ones(S0)),
dtype = np.int)
Rtall = np.zeros(S0**3, dtype = np.complex)
Rtall[II] = np.reshape(self.RPOD.T, (S0**2,))
Rtall = np.reshape(Rtall, (S0, S0**2)).T
B = Rtall.dot(Phi[:, : N1])
Z = copy(B)
Z = np.reshape(Z.T, (S0 * (N1), S0)).T
for j in range(2):
Z = Z - Phi[:, : M1].dot(Phi[:, : M1].conj().T.dot(
np.multiply(self.ws, Z)))
Z = np.reshape(Z.T, (N1, S0**2)).T
else:
ker = self.solSnapshots.conj().T.dot(self.energyNormMat.dot(
self.solSnapshots))
WPhi = np.reshape(np.multiply(self.ws, Phi[:, : M1]).T,
(S0 * M1)).conj()[:, None]
Y = np.multiply(WPhi, np.kron(np.ones((M1, 1)), Phi[:, : N1]))
Ylarge = np.reshape(Y.T, (M1 * N1, S0)).T
B = np.reshape(ker.dot(Ylarge).T, (N1, M1 * S0)).T
D = np.multiply(self.ws, np.diag(ker)[:, None])
D = Phi[:, : N1].conj().T.dot(np.multiply(D, Phi[:, : N1]))
Z = B.conj().T.dot(Y) - D
_, ev, eV = np.linalg.svd(Z, full_matrices=False)
eV = eV.conj().T
phi = eV[:, np.argmin(ev)]
if self.POD:
c = Phi[:, : M1].conj().T.dot(np.multiply(self.ws,
np.reshape(B.dot(phi).T, (S0, S0)).T))
else:
c = np.reshape(Y.dot(phi).T, (M1, S0))
polybase = np.zeros((max(M1, N1), max(M1, N1)))
polybase[0, 0] = 1
polybase[1, 1] = 1
for j in range(2, max(M1, N1)):
if self.polyBasis == "CHEBYSHEV":
polybase[1 : j + 1, j] = 2 * polybase[: j, j - 1]
polybase[: j + 1, j] = (polybase[: j + 1, j]
- polybase[: j + 1, j - 2])
elif self.polyBasis == "GAUSSLEGENDRE":
polybase[1 : j + 1, j] = ((2 * j - 1) / j
* polybase[: j, j - 1])
polybase[: j + 1, j] = (polybase[: j + 1, j]
- ((j - 1) / j * polybase[: j + 1, j - 2]))
if self.polyBasis == "CHEBYSHEV":
polybase[0, 0] = .5 ** .5
elif self.polyBasis == "GAUSSLEGENDRE":
polybase = np.multiply(np.power(
np.arange(.5, max(M1, N1)), .5), polybase)
self.P = self.solSnapshots.dot(polybase[: M1, : M1].dot(c).T)
self.Q = polybase[: N1, : N1].dot(phi)
self.approxParameters = {"N" : self.N, "M" : self.M, "S" : S0}
self.lastApproxParameters = copy(self.approxParameters)
def _cleanup(self):
"""Cleanup Pade' denominator by removing unwanted poles."""
if self.N == 0: return
poles = np.roots(self.Q[::-1]) + self.k0
NpolesOld = len(poles)
poles = poles[self.cleanupParameters['boolCondition'](poles)]
if self.cleanupParameters['residueCheck']:
resR = self.cleanupParameters['residueRadius']
resTol = self.cleanupParameters['residueTol']
residues = np.zeros_like(poles)
for l, pole in enumerate(poles):
resV = np.zeros(self.solDerivatives.shape[0],
dtype = np.complex)
NPoints = self.cleanupParameters['residueNPoints']
for theta in 2 * PI * np.arange(NPoints) / NPoints:
deltapole = resR * np.exp(1.j * theta)
ksample = pole + deltapole
self.solveHF(ksample)
resV = deltapole / NPoints * (resV + self.uHF)
residues[l] = np.abs(resV.dot(
self.energyNormMatrix.dot(resV).conj())) ** .5
poles = poles[residues >= resTol]
NpolesNew = len(poles)
if NpolesOld > NpolesNew:
if self.verboseRobust >= 1:
sSing = "s" * (NpolesOld - NpolesNew > 1)
print(("Identified {0} pole{1} to be removed.\n"
"Reducing N from {2} to {3}.").format(
NpolesOld - NpolesNew, sSing, NpolesOld, NpolesNew))
self.N = NpolesNew
newQ = np.polyfit(np.append(poles - self.k0, 0.),
np.append(np.zeros(NpolesNew), 1.), NpolesNew)
self.Q = newQ[::-1] / np.linalg.norm(newQ)
def radiusPade(self, k:complex):
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
k: Target wavenumber.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
return (k - self.k0) / self.approxRadius
def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Real and imaginary parts of approximant.
"""
self.setupApprox()
powerlist = np.power(self.radiusPade(k), range(max(self.M, self.N) + 1))
self.uApp = (self.P.dot(powerlist[:self.M + 1])
/ self.Q.dot(powerlist[:self.N + 1]))
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()
return np.roots(self.Q[::-1]) * self.approxRadius + self.k0 * centered
c4science · Help