Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/main/ROMApproximantTaylorPade.py
4c3d184239eefeature_greedy_coarseni…

Last Change
Raw File

ROMApproximantTaylorPade.py
View Options

#!/usr/bin/python
from __future__ import print_function
from copy import copy
import warnings
import numpy as np
import utilities
from ROMApproximantTaylor import ROMApproximantTaylor
PI = np.pi
class ROMApproximantTaylorPade(ROMApproximantTaylor):
"""
ROM single-point fast Pade' approximant 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;
- setupDerivativeComputation: setup HF problem data to compute j-th
solution derivative at k0;
- 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.
k0(optional): Default parameter. Defaults to 0.
w(optional): Weight for norm computation. If None, set to Re(k0).
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;
- 'rho': weight for computation of original Pade' approximant;
defaults to inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- '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;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0;
- '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.
plotDer(optional): Whether to plot derivatives of the Helmholtz
solution map. Defaults to False.
equilibration(optional): Whether to apply equilibration to Gram matrix.
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.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
k0: Default parameter.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
defaults to inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'robustTol': tolerance for robust Pade' denominator management;
- 'sampleType': label of sampling type.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
robustTol: tolerance for robust Pade' denominator management.
sampleType: Label of sampling type.
plotDer: Whether to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
equilibration: Whether to apply equilibration to Gram matrix.
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.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
equilPowers: Numpy 1D (column) array containing equilibration powers.
Does not exist if equilibration is disabled.
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", "robustTol", "rho"]
def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
k0 : complex = 0, w : float = None,
approxParameters : dict = {}, plotDer : list = [],
equilibration : bool = False, verboseRobust : int = 1,
cleanupParameters : dict = {}):
self.equilibration = equilibration
self.verboseRobust = verboseRobust
self.cleanupParameters = cleanupParameters
ROMApproximantTaylor.__init__(self, HFEngine, HSEngine, k0, w,
approxParameters, plotDer)
def parameterList(self) -> list:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximantTaylor.parameterList(self)
+ ROMApproximantTaylorPade.extraApproxParameters)
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N, E,
Emax and robustTol.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = utilities.purgeDict(approxParams,
ROMApproximantTaylorPade.parameterList(self),
dictname = self.name() + ".approxParameters")
keyList = list(approxParameters.keys())
if "M" in keyList:
self.M = 0 #to avoid warnings if M and E are changed simultaneously
if "N" in keyList:
self.N = 0 #to avoid warnings if N and E are changed simultaneously
approxParametersCopy = utilities.purgeDict(approxParameters,
ROMApproximantTaylorPade.extraApproxParameters,
True, True)
ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy)
if "rho" in keyList:
self.rho = approxParameters["rho"]
elif hasattr(self, "rho"):
self.rho = self.rho
else:
self.rho = np.inf
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 "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
@property
def rho(self):
"""Value of rho."""
return self._rho
@rho.setter
def rho(self, rho):
self._rho = np.abs(rho)
self._approxParameters["rho"] = self.rho
@property
def M(self):
"""Value of M. Its assignment may change Emax and E."""
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, "Emax") and self.Emax < self.M:
warnings.warn("Prescribed Emax is too small. Updating Emax to M.")
self.Emax = self.M
if hasattr(self, "E") and self.E < self.M:
warnings.warn("Prescribed E is too small. Updating E to M.")
self.E = self.M
@property
def N(self):
"""Value of N. Its assignment may change Emax."""
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, "Emax") and self.Emax < self.N:
warnings.warn("Prescribed Emax is too small. Updating Emax to N.")
self.Emax = self.N
if hasattr(self, "E") and self.E < self.N:
warnings.warn("Prescribed E is too small. Updating E to N.")
self.E = self.N
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
warnings.warn(("Overriding prescribed negative robustness "
"tolerance to 0."))
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
if hasattr(self, "M") and self.M > self.E:
warnings.warn("Prescribed E is too small. Updating E to M.")
self.E = self.M
if hasattr(self, "N") and self.N > self.E:
warnings.warn("Prescribed E is too small. Updating E to N.")
self.E = self.N
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warnings.warn("Prescribed Emax is too small. Updating Emax to E.")
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
vals, label = [0] * 3, {0:"E", 1:"M", 2:"N"}
if hasattr(self, "E"): vals[0] = self.E
if hasattr(self, "M"): vals[1] = self.M
if hasattr(self, "N"): vals[2] = self.N
idxmax = np.argmax(vals)
if vals[idxmax] > Emax:
warnings.warn("Prescribed Emax is too small. Updating Emax to {}."\
.format(label[idxmax]))
self.Emax = vals[idxmax]
else:
self._Emax = Emax
self._approxParameters["Emax"] = Emax
if (EmaxOld > self.Emax and self.solDerivatives is not None
and self.HArnoldi is not None and self.RArnoldi is not None):
self.solDerivatives = self.solDerivatives[:, : self.Emax + 1]
if self.sampleType == "ARNOLDI":
self.HArnoldi = self.HArnoldi[: self.Emax + 1,
: self.Emax + 1]
self.RArnoldi = self.RArnoldi[: self.Emax + 1,
: self.Emax + 1]
else:
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' approximant.
See Section 4 in Fast Pade' paper.
SVD-based robust eigenvalue management.
"""
self.computeDerivatives()
if not self.checkComputedApprox():
while True:
if self.POD:
if self.sampleType == "ARNOLDI":
ev, eV = self.findeveVGArnoldi()
else:
ev, eV = self.findeveVGQR()
else:
self.buildG()
ev, eV = np.linalg.eigh(self.G)
ts = self.robustTol * np.linalg.norm(ev)
Nnew = np.sum(np.abs(ev) >= ts)
diff = self.N - Nnew
if diff <= 0: break
Enew = self.E - diff
Mnew = min(self.M, Enew)
if Mnew == self.M:
strM = ""
else:
strM = ", M from {0} to {1},".format(self.M, Mnew)
if self.verboseRobust >= 1:
print(("Smallest {0} eigenvalues below tolerance.\n"
"Reducing N from {1} to {2}{5} and E from {3} to "
"{4}.").format(diff + 1, self.N, Nnew,
self.E, Enew, strM))
newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
self.approxParameters = newParameters
self.Q = eV[::-1, 0]
if self.equilibration:
self.Q = np.multiply(self.equilPowers.flatten()[::-1], self.Q)
self._cleanup()
QToeplitz = np.zeros((self.Emax + 1, self.M + 1),
dtype = np.complex)
for i in range(len(self.Q)):
rng = np.arange(self.M + 1 - i)
QToeplitz[rng, rng + i] = self.Q[i]
if self.sampleType == "ARNOLDI":
QToeplitz = self.RArnoldi.dot(QToeplitz)
self.P = self.solDerivatives.dot(QToeplitz)
self.lastApproxParameters = copy(self.approxParameters)
def buildG(self):
"""Assemble Pade' denominator matrix."""
if self.N == 0:
self.G = np.array([[1]], dtype = np.complex)
return
self.computeDerivatives()
if np.isinf(self.rho):
lgSize = self.N + 1
else:
lgSize = self.E + self.N - self.M
Gshift = self.E - lgSize + 1
largeG = np.zeros((lgSize, lgSize), dtype=np.complex)
largeG = self.solDerivatives[:, Gshift:self.E+1].conj().T.dot(
self.energyNormMatrix.dot(self.solDerivatives[:, Gshift:self.E+1]))
if np.isinf(self.rho):
self.G = largeG
else:
self.G = np.zeros((self.N + 1, self.N + 1), dtype=np.complex)
for k in range(self.E - self.M):
self.G = (self.G + self.rho ** (2 * k)
* largeG[k : k+self.N+1, k : k+self.N+1])
if self.equilibration:
Gd = np.diag(self.G)
gamma = np.average(np.abs(np.divide(Gd[1:], Gd[:-1])),
weights = np.power(1.2, np.arange(self.N)))**.5
self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
self.G = np.multiply(self.equilPowers,
np.multiply(self.equilPowers, self.G).T).T
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 findeveVGQR(self):
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor. See ``Householder triangularization of a
quasimatrix'', L.Trefethen, 2008 for QR algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self.computeDerivatives()
A = copy(self.solDerivatives[:, (self.E - self.N) : (self.E + 1)])
M = self.energyNormMatrix
E = np.zeros_like(A, dtype = np.complex)
R = np.zeros((self.N + 1, self.N + 1), dtype = np.complex)
for k in range(self.N + 1):
E[:, k] = (np.random.randn(E.shape[0])
+ 1.j * np.random.randn(E.shape[0]))
if k > 0:
for times in range(2):
E[:, k] = E[:, k] - E[:, :k].dot(
(E[:, k].conj().dot(M.dot(E[:, :k]))).conj())
E[:, k] = E[:, k] / (E[:, k].conj().dot(M.dot(E[:, k]))) ** .5
R[k, k] = np.abs(A[:, k].conj().dot(M.dot(A[:, k]))) ** .5
alpha = E[:, k].conj().dot(M.dot(A[:, k]))
if np.isclose(np.abs(alpha), 0.): s = 1.
else: s = - alpha / np.abs(alpha)
E[:, k] = s * E[:, k]
v = R[k, k] * E[:, k] - A[:, k]
for times in range(2):
v = v - E[:, :k].dot((v.conj().dot(M.dot(E[:, :k]))).conj())
sigma = np.abs(v.conj().dot(M.dot(v))) ** .5
if np.isclose(sigma, 0.): v = E[:, k]
else: v = v / sigma
J = np.arange(k + 1, self.N + 1)
vtQ = v.conj().dot(M.dot(A[:, J]))
A[:, J] = A[:, J] - 2 * np.outer(v, vtQ)
R[k, J] = E[:, k].conj().dot(M.dot(A[:, J]))
A[:, J] = A[:, J] - np.outer(E[:, k], R[k, J])
if self.equilibration:
Rd = np.diag(R)
gamma = np.average(np.abs(np.divide(Rd[1:], Rd[:-1])),
weights = np.power(1.2, np.arange(self.N)))**.5
self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
R = np.multiply(self.equilPowers, R.T).T
_, s, V = np.linalg.svd(R, full_matrices = False)
return s[::-1], V.conj().T[:, ::-1]
def findeveVGArnoldi(self):
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor of solution derivatives orthogonalized by
Arnoldi algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self.computeDerivatives()
if self.sampleType != "ARNOLDI":
raise Exception(("Eigensolver different from 'ARNOLDI'."
"Arnoldi eigenproblem solver cannot be called."))
R = self.RArnoldi[: self.E + 1, self.E - self.N : self.E + 1]
if self.equilibration:
Rd = np.diag(R)
gamma = np.average(np.abs(np.divide(Rd[1:], Rd[:-1])),
weights = np.power(1.2, np.arange(self.N)))**.5
self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
R = np.multiply(self.equilPowers, R.T).T
_, s, V = np.linalg.svd(R, full_matrices = False)
return s[::-1], V.conj().T[:, ::-1]
def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Approximant as numpy complex vector.
"""
self.setupApprox()
powerlist = np.power(k - self.k0, 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.k0 * centered
c4science · Help