Page MenuHomec4science
Files F90026510

ROMApproximantLagrangeVF.py
No OneTemporary
Actions

File Metadata

Created
Mon, Oct 28, 14:49

ROMApproximantLagrangeVF.py
View Options

#!/usr/bin/python
from copy import copy
import warnings
import numpy as np
import utilities
from ROMApproximantLagrange import ROMApproximantLagrange
from RROMPyTypes import Np1D, ListAny
class ROMApproximantLagrangeVF(ROMApproximantLagrange):
"""
ROM Lagrange VF interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver. Should have members:
- Vdim: domain dimension;
- energyNormMatrix: assemble sparse matrix (in CSC format)
associated to weighted H10 norm;
- problemData: list of HF problem data (excluding mu);
- setProblemData: set HF problem data and mu to prescribed values;
- getLSBlocks: get algebraic system blocks.
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.
mus(optional): Boundaries of snapshot positions. Defaults to
np.array([0, 1]).
ws(optional): Array of snapshot weigths. Defaults to
np.ones_like(self.mus).
w(optional): Weight for norm computation. If None, set to
Re(np.mean(self.mus)). 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 VF interpolant numerator; defaults to 0;
- 'N': degree of VF interpolant denominator; defaults to 0.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Boundaries of snapshot positions.
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 VF interpolant numerator;
- 'N': degree of VF interpolant denominator;
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
autoNodeTypes: List of possible values for autoNode.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
S: Number of solution snapshots over which current approximant is
based upon.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
autoNode: Type of nodes, if automatically generated. Otherwise None.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
"""
extraApproxParameters = ["M", "N"]
def parameterList(self) -> ListAny:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return (ROMApproximantLagrange.parameterList(self)
+ ROMApproximantLagrangeVF.extraApproxParameters)
@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,
ROMApproximantLagrangeVF.parameterList(self),
dictname = self.name() + ".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,
ROMApproximantLagrangeVF.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
@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.",
stacklevel = 2)
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.",
stacklevel = 2)
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]), stacklevel = 2)
self.Emax = vals[idxmax] + 1
else:
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
def approxRadius(self) -> float:
"""Sample radius."""
return np.max(self.mu0 - self.mus)
def setupApprox(self):
"""Compute VF interpolant. SVD-based robust eigenvalue management."""
if not self.checkComputedApprox():
S0 = self.S
M1 = self.M + 1
N1 = self.N + 1
if self.autoNode not in self.autoNodeTypes[1 :]:
raise Exception("Manual sample selection not yet implemented.")
self.computeSnapshots()
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.radiusVF(self.mus), c)
elif self.polyBasis == "GAUSSLEGENDRE":
Phi[:, j] = (j + .5) ** .5 * np.polynomial.legendre.legval(
self.radiusVF(self.mus), 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): #twice is enough!
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)
phi = eV[-1, :].conj()
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 radiusVF(self, mu:complex) -> float:
"""
Compute translated radius to be plugged into VF approximant.
Args:
mu: Target parameter.
Returns:
Translated radius to be plugged into VF approximant.
"""
return (mu - self.mu0) / self.approxRadius()
def evalApprox(self, mu:complex) -> Np1D:
"""
Evaluate VF approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Real and imaginary parts of approximant.
"""
self.setupApprox()
powerlist = np.power(self.radiusVF(mu), 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) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
return np.roots(self.Q[::-1]) * self.approxRadius() + self.mu0

Event Timeline

c4science · Help