diff --git a/examples/fenics/HelmholtzPadeTaylorApproximant.py b/examples/fenics/HelmholtzPadeTaylorApproximant.py
index 0b39054..a03126b 100644
--- a/examples/fenics/HelmholtzPadeTaylorApproximant.py
+++ b/examples/fenics/HelmholtzPadeTaylorApproximant.py
@@ -1,69 +1,69 @@
import numpy as np
from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.hsengines.fenics import HSEngine as HS
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
-testNo = 1
+testNo = 2
if testNo == 1:
params = {'N':4, 'M':3, 'E':4, 'sampleType':'Arnoldi', 'POD':True}
z0 = 12+.5j
ztar = 11
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
plotter = HS(solver.V)
approx = Pade(solver, plotter, mu0 = z0, approxParameters = params)
approx.plotApp(ztar, name = 'u_Pade''')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 2:
params = {'N':7, 'M':6, 'E':7, 'sampleType':'Arnoldi', 'POD':True}
z0 = 16
ztar = 14
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
plotter = HS(solver.V)
approx = Pade(solver, plotter, mu0 = z0, approxParameters = params,
plotDer = 'ALL')
approx.plotApp(ztar, name = 'u_Pade''')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 3:
params = {'N':8, 'M':7, 'E':8, 'sampleType':'Krylov', 'POD':True}
k0 = 3
ktar = 4+.5j
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30)
plotter = HS(solver.V)
approx = Pade(solver, plotter, mu0 = k0, approxParameters = params)
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
diff --git a/examples/fenics/HelmholtzRBTaylorApproximant.py b/examples/fenics/HelmholtzRBTaylorApproximant.py
index 203e312..90393f7 100644
--- a/examples/fenics/HelmholtzRBTaylorApproximant.py
+++ b/examples/fenics/HelmholtzRBTaylorApproximant.py
@@ -1,69 +1,69 @@
import numpy as np
from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.hsengines.fenics import HSEngine as HS
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
-testNo = 1
+testNo = 2
if testNo == 1:
params = {'E':4, 'R':4, 'sampleType':'Arnoldi', 'POD':True}
z0 = 12+.5j
ztar = 11
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
plotter = HS(solver.V)
approx = RB(solver, plotter, mu0 = z0, approxParameters = params)
approx.plotApp(ztar, name = 'u_RB')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}'.format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles RB:')
print(approx.getPoles())
############
elif testNo == 2:
- params = {'E':7, 'R':7, 'sampleType':'Arnoldi', 'POD':True}
+ params = {'E':20, 'R':21, 'sampleType':'Arnoldi', 'POD':True}
z0 = 16
ztar = 14
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
plotter = HS(solver.V)
approx = RB(solver, plotter, mu0 = z0, approxParameters = params,
plotDer = 'ALL')
approx.plotApp(ztar, name = 'u_RB')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles RB:')
print(approx.getPoles())
############
elif testNo == 3:
params = {'E':8, 'R':8, 'sampleType':'Krylov', 'POD':True}
k0 = 3
ktar = 4+.5j
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30)
plotter = HS(solver.V)
approx = RB(solver, plotter, mu0 = k0, approxParameters = params)
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles RB:')
print(approx.getPoles())
diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py
index 74face3..00412af 100644
--- a/rrompy/reduction_methods/base/__init__.py
+++ b/rrompy/reduction_methods/base/__init__.py
@@ -1,27 +1,29 @@
# 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 .
#
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
from rrompy.reduction_methods.base.parameter_sweeper import ParameterSweeper
+from rrompy.reduction_methods.base.pod_engine import PODEngine
__all__ = [
'GenericApproximant',
- 'ParameterSweeper'
+ 'ParameterSweeper',
+ 'PODEngine'
]
diff --git a/rrompy/reduction_methods/base/pod_engine.py b/rrompy/reduction_methods/base/pod_engine.py
new file mode 100644
index 0000000..1ba4793
--- /dev/null
+++ b/rrompy/reduction_methods/base/pod_engine.py
@@ -0,0 +1,133 @@
+# 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 .
+#
+
+import numpy as np
+from copy import copy
+from rrompy.utilities.types import Np1D, Np2D, Tuple
+
+__all__ = ['PODEngine']
+
+class PODEngine:
+ """
+ POD engine for general matrix orthogonalization.
+
+ Args/Attributes:
+ M: Square matrix associated to scalar product.
+ """
+
+ def __init__(self, M:Np2D):
+ self.M = M
+
+ def norm(self, a:Np1D) -> float:
+ """Compute norm of a vector."""
+ return np.sqrt(np.abs(a.conj().T.dot(self.M.dot(a))))
+
+ def GS(self, a:Np1D, Q:Np2D, n : int = None) -> Tuple[Np1D, Np1D]:
+ """
+ Compute 1 Gram-Schmidt step with given projector.
+
+ Args:
+ a: vector to be projected;
+ Q: orthogonal projection matrix;
+ n: number of columns of Q to be considered.
+
+ Returns:
+ Resulting normalized vector, coefficients of a wrt the updated
+ basis.
+ """
+ if n is None:
+ n = Q.shape[1]
+ Q = Q[:, : n]
+ r = np.zeros((n + 1,), dtype = a.dtype)
+ if n > 0:
+ for j in range(2): # twice is enough!
+ nu = Q.conj().T.dot(self.M.dot(a))
+ a = a - Q.dot(nu)
+ r[:-1] = r[:-1] + nu
+ r[-1] = self.norm(a)
+ if np.isclose(np.abs(r[-1]), 0.):
+ r[-1] = 1.
+ a = a / r[-1]
+ return a, r
+
+ def QRGramSchmidt(self, A:Np2D,
+ only_R : bool = False) -> Tuple[Np1D, Np1D]:
+ """
+ Compute QR decomposition of a matrix through Gram-Schmidt method.
+
+ Args:
+ A: matrix to be decomposed;
+ only_R(optional): whether to skip reconstruction of Q; defaults to
+ False.
+
+ Returns:
+ Resulting orthogonal and upper-triangular factors.
+ """
+ N = A.shape[1]
+ Q = np.zeros_like(A, dtype = A.dtype)
+ R = np.zeros((N, N), dtype = A.dtype)
+ for k in range(N):
+ Q[:, k], R[: k + 1, k] = self.GS(A[:, k], Q, k)
+ if only_R:
+ return R
+ return Q, R
+
+ def QRHouseholder(self, A:Np2D, Q0 : Np2D = None,
+ only_R : bool = False) -> Tuple[Np1D, Np1D]:
+ """
+ Compute QR decomposition of a matrix through Householder method.
+
+ Args:
+ A: matrix to be decomposed;
+ Q0(optional): initial orthogonal guess for Q; defaults to random;
+ only_R(optional): whether to skip reconstruction of Q; defaults to
+ False.
+
+ Returns:
+ Resulting (orthogonal and )upper-triangular factor(s).
+ """
+ B = copy(A)
+ N = B.shape[1]
+ V = np.zeros_like(B, dtype = B.dtype)
+ R = np.zeros((N, N), dtype = B.dtype)
+ if Q0 is None:
+ Q = np.zeros_like(B, dtype = B.dtype) + np.random.randn(*(B.shape))
+ else:
+ Q = copy(Q0)
+ for k in range(N):
+ if Q0 is None:
+ Q[:, k], _ = self.GS(Q[:, k], Q, k)
+ a = B[:, k]
+ R[k, k] = self.norm(a)
+ alpha = Q[:, k].conj().T.dot(self.M.dot(a))
+ if np.isclose(np.abs(alpha), 0.): s = 1.
+ else: s = - alpha / np.abs(alpha)
+ Q[:, k] = s * Q[:, k]
+ V[:, k], _ = self.GS(R[k, k] * Q[:, k] - a, Q, k)
+ J = np.arange(k + 1, N)
+ vtB = V[:, k].conj().dot(self.M.dot(B[:, J]))
+ B[:, J] = B[:, J] - 2 * np.outer(V[:, k], vtB)
+ R[k, J] = Q[:, k].conj().T.dot(self.M.dot(B[:, J]))
+ B[:, J] = B[:, J] - np.outer(Q[:, k], R[k, J])
+ if only_R:
+ return R
+ for k in range(N - 1, -1, -1):
+ J = np.arange(k, N)
+ vtQ = V[:, k].conj().T.dot(self.M.dot(Q[:, J]))
+ Q[:, J] = Q[:, J] - 2 * np.outer(V[:, k], vtQ)
+ return Q, R
diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
index 6f7865d..d6c1033 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
@@ -1,250 +1,251 @@
# 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 .
#
from copy import copy
import warnings
import numpy as np
import scipy as sp
from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import GenericApproximantLagrange
from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng
from rrompy.utilities import purgeDict
__all__ = ['ApproximantLagrangeRB']
class ApproximantLagrangeRB(GenericApproximantLagrange):
"""
ROM RB 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 mu);
- setProblemData: set HF problem data and mu to prescribed values;
- getLSBlocks: get algebraic system blocks;
- liftDirichletData: perform lifting of Dirichlet BC as numpy
complex 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.
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];
- 'R': rank for Galerkin projection; defaults to S.
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: Array of snapshot parameters.
ws: Array of snapshot weigths (unused).
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
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;
- 'R': rank for Galerkin projection.
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.
R: Rank for Galerkin projection.
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.
projMat: Projection matrix for RB system assembly.
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.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
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.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix wrt mu.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS wrt mu.
"""
def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
approxParameters = approxParameters,
plotSnap = plotSnap, plotSnapSpecs = plotSnapSpecs)
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@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):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters")
approxParametersCopy = purgeDict(approxParameters, ["R"], True, True)
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.S
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "S") and self.S < self.R:
warnings.warn("Prescribed S is too small. Updating S to R.",
stacklevel = 2)
self.S = self.R
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 setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
self.computeSnapshots()
if self.POD:
U, _, _ = np.linalg.svd(self.RPOD, full_matrices = False)
self.projMat = self.solSnapshots.dot(U[:, : self.R])
else:
self.projMat = self.solSnapshots[:, : self.R]
self.lastApproxParameters = copy(self.approxParameters)
self.assembleReducedSystem()
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
self.setupApprox()
projMatH = self.projMat.T.conjugate()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
ARBmu = self.ARBs[0][:self.R,:self.R]
bRBmu = self.bRBs[0][:self.R]
for j in range(1, len(self.ARBs)):
ARBmu = ARBmu + np.power(mu, j) * self.ARBs[j][:self.R, :self.R]
for j in range(1, len(self.bRBs)):
bRBmu = bRBmu + np.power(mu, j) * self.bRBs[j][:self.R]
return self.projMat[:, :self.R].dot(np.linalg.solve(ARBmu, bRBmu))
def evalApprox(self, mu:complex) -> Np1D:
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Numpy 1D array with approximant dofs.
"""
self.setupApprox()
self.uApp = self.solveReducedSystem(mu)
return self.uApp
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
try:
return sp.linalg.eigvals(A, B)
except:
warnings.warn("Generalized eig algorithm did not converge.",
stacklevel = 2)
x = np.empty(A.shape[0])
x[:] = np.NaN
return x
+
diff --git a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
index 1ecd211..5eb3dd0 100644
--- a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
+++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
@@ -1,229 +1,223 @@
# 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 .
#
import numpy as np
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
+from rrompy.reduction_methods.base.pod_engine import PODEngine
from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng
from rrompy.utilities import purgeDict
__all__ = ['GenericApproximantLagrange']
class GenericApproximantLagrange(GenericApproximant):
"""
ROM Lagrange interpolant computation for parametric problems (ABSTRACT).
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 mu);
- setProblemData: set HF problem data and k 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.
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;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1].
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: 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 snapshots current approximant relies upon;
- 'sampler': sample point generator.
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.
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.
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.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
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.
"""
def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["S", "sampler"]
super().__init__(HFEngine = HFEngine, HSEngine = HSEngine,
mu0 = mu0, approxParameters = approxParameters)
self.plotSnap = plotSnap
self.plotSnapSpecs = plotSnapSpecs
self.resetSamples()
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
musOld = self.mus if hasattr(self, 'mus') else None
self._mus = np.array(mus)
_, musCounts = np.unique(self._mus, return_counts = True)
if len(np.where(musCounts > 1)[0]) > 0:
raise Exception("Repeated sample points not allowed.")
if (musOld is None or len(self.mus) != len(musOld)
or not np.allclose(self.mus, musOld, 1e-14)):
self.resetSamples()
self.autoNode = None
@property
def approxParameters(self):
"""Value of approximant parameters. Its assignment may change S."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters")
approxParametersCopy = purgeDict(approxParameters, ["S", "sampler"],
True, True)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "S" in keyList:
self.S = approxParameters["S"]
elif hasattr(self, "S"):
self.S = self.S
else:
self.S = 2
if "sampler" in keyList:
self.sampler = approxParameters["sampler"]
elif not hasattr(self, "S"):
from rrompy.sampling import QuadratureSampler
self.sampler = QuadratureSampler([0., 1.], "UNIFORM")
del QuadratureSampler
@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
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
@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:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
self.solSnapshots = None
self.RPOD = None
def computeSnapshots(self):
"""
Compute snapshots of solution map.
"""
if self.solSnapshots is None:
self.mus, self.ws = self.sampler.generatePoints(self.S)
self.mus = np.array([x[0] for x in self.mus])
for j, mu in enumerate(self.mus):
self.solveHF(mu)
self.HSEngine.plot(self.uHF, name = "u({:.4f})".format(mu),
what = self.plotSnap, **self.plotSnapSpecs)
self.manageSnapshots(self.uHF, j)
def manageSnapshots(self, u:Np1D, pos:int):
"""
Store snapshots of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
"""
if pos == 0:
self.solSnapshots = np.empty((u.shape[0], self.S),
dtype = np.complex)
self.solSnapshots[:, pos] = u
if self.POD:
if pos == 0:
+ self.PODEngine = PODEngine(self.energyNormMatrix)
self.RPOD = np.eye(self.S, dtype = np.complex)
- beta = 1
- for j in range(2): #twice is enough!
- nu = self.solSnapshots[:, pos].conj().dot(
- self.energyNormMatrix.dot(self.solSnapshots[:, : pos])
- ).conj()
- self.RPOD[: pos, pos] = self.RPOD[: pos, pos] + beta * nu
- eta = (self.solSnapshots[:, pos]
- - self.solSnapshots[:, : pos].dot(nu))
- beta = eta.conj().dot(self.energyNormMatrix.dot(eta))**.5
- self.solSnapshots[:, pos] = eta / beta
- self.RPOD[pos, pos] = beta * self.RPOD[pos, pos]
+ self.solDerivatives[:, pos], self.RPOD[: pos + 1, pos] = (
+ self.PODEngine.GS(self.solDerivatives[:, pos],
+ self.solDerivatives, pos))
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.solSnapshots is not None and super().checkComputedApprox()
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
index 92d2d6c..831f930 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
@@ -1,458 +1,433 @@
# 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 .
#
from copy import copy
import warnings
import numpy as np
from rrompy.reduction_methods.taylor.generic_approximant_taylor import GenericApproximantTaylor
+from rrompy.reduction_methods.base.pod_engine import PODEngine
from rrompy.utilities.types import Np1D, Np2D, ListAny, Tuple, DictAny
from rrompy.utilities.types import HFEng, HSEng
from rrompy.utilities import purgeDict
__all__ = ['ApproximantTaylorPade']
class ApproximantTaylorPade(GenericApproximantTaylor):
"""
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 mu);
- setProblemData: set HF problem data and mu to prescribed values;
- setupDerivativeComputation: setup HF problem data to compute j-th
solution derivative at mu0;
- 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.
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;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.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): What to plot of derivatives of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives 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.
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;
- 'rho': weight for computation of original Pade' 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.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
RArnoldi: Right factor of derivatives QR factorization. If no Arnoldi,
set to None.
HArnoldi: Hessenberg factor of derivatives Arnoldi algorithm. If no
Arnoldi, set to None.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
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: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
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.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
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.
"""
def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
approxParameters : DictAny = {}, plotDer : ListAny = [],
plotDerSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["M", "N", "robustTol", "rho"]
super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
approxParameters = approxParameters,
plotDer = plotDer, plotDerSpecs = plotDerSpecs)
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters")
approxParametersCopy = purgeDict(approxParameters,
["M", "N", "robustTol", "rho"],
True, True)
keyList = list(approxParameters.keys())
if "rho" in keyList:
self._rho = approxParameters["rho"]
elif hasattr(self, "rho"):
self._rho = self.rho
else:
self._rho = np.inf
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
self.rho = self._rho
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
self._ignoreParWarnings = True
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
del self._ignoreParWarnings
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 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 not hasattr(self, "_ignoreParWarnings"):
self.checkMNEEmax()
@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
self.checkMNEEmax()
def checkMNEEmax(self):
"""Check consistency of M, N, E, and Emax."""
M = self.M if hasattr(self, "_M") else 0
N = self.N if hasattr(self, "_N") else 0
E = self.E if hasattr(self, "_E") else M + N
Emax = self.Emax if hasattr(self, "_Emax") else M + N
if self.rho == np.inf:
if Emax < max(N, M):
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(M, N)."), stacklevel = 3)
self.Emax = max(N, M)
if E < max(N, M):
warnings.warn(("Prescribed E is too small. Updating E to "
"max(M, N)."), stacklevel = 3)
self.E = max(N, M)
else:
if Emax < N + M:
warnings.warn(("Prescribed Emax is too small. Updating "
"Emax to M + N."), stacklevel = 3)
self.Emax = self.N + M
if E < N + M:
warnings.warn(("Prescribed E is too small. Updating E to "
"M + N."), stacklevel = 3)
self.E = self.N + M
@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."), stacklevel = 2)
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
self.checkMNEEmax()
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."), stacklevel = 2)
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
if hasattr(self, "_N"): N = self.N
else: N = 0
if hasattr(self, "_M"): M = self.M
else: M = 0
if hasattr(self, "_E"): E = self.E
else: E = 0
if self.rho == np.inf:
if max(N, M, E) > Emax:
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(N, M, E)."), stacklevel = 2)
Emax = max(N, M, E)
else:
if max(N + M, E) > Emax:
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(N + M, E)."), stacklevel = 2)
Emax = max(N + M, E)
self._Emax = Emax
self._approxParameters["Emax"] = Emax
if (EmaxOld > self.Emax and self.solDerivatives is not None
and hasattr(self, 'HArnoldi') and hasattr(self, 'RArnoldi')
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()
def setupApprox(self):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
self.computeDerivatives()
while True:
if self.POD:
ev, eV = self.findeveVGQR()
else:
ev, eV = self.findeveVGExplicit()
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)
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]
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."""
self.computeDerivatives()
if self.rho == np.inf:
if self.sampleType == "KRYLOV":
DerE = self.solDerivatives[:, self.E - self.N : self.E + 1]
self.G = DerE.conj().T.dot(self.energyNormMatrix.dot(DerE))
else:
RArnE = self.RArnoldi[: self.E + 1,
self.E - self.N : self.E + 1]
self.G = RArnE.conj().T.dot(RArnE)
else:
if self.sampleType == "KRYLOV":
DerE = self.solDerivatives[:, self.M - self.N + 1 : self.E + 1]
Gbig = DerE.conj().T.dot(self.energyNormMatrix.dot(DerE))
else:
RArnE = self.RArnoldi[: self.E + 1,
self.M - self.N + 1 : self.E + 1]
Gbig = RArnE.conj().T.dot(RArnE)
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) * Gbig[k : k + self.N + 1,
k : k + self.N + 1])
def findeveVGExplicit(self) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
self.buildG()
ev, eV = np.linalg.eigh(self.G)
return ev, eV
def findeveVGQR(self) -> Tuple[Np1D, Np2D]:
"""
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()
if self.sampleType == "KRYLOV":
if self.rho == np.inf:
A = copy(self.solDerivatives[:, self.E - self.N : self.E + 1])
else:
A = copy(self.solDerivatives[:,
self.M - self.N + 1 : self.E + 1])
- M = self.energyNormMatrix
- E = np.zeros_like(A, dtype = np.complex)
- R = np.zeros((A.shape[1], A.shape[1]), dtype = np.complex)
- for k in range(A.shape[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): #twice is enough!
- Ekproj = E[:, k].conj().dot(M.dot(E[:, :k])).conj()
- E[:, k] = E[:, k] - E[:, :k].dot(Ekproj)
- 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): #twice is enough!
- vproj = v.conj().dot(M.dot(E[:, :k])).conj()
- v = v - E[:, :k].dot(vproj)
- 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, A.shape[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])
+ self.PODEngine = PODEngine(self.energyNormMatrix)
+ R = self.PODEngine.QRHouseholder(A, only_R = True)
else:
if self.rho == np.inf:
R = self.RArnoldi[: self.E + 1, self.E - self.N : self.E + 1]
else:
R = self.RArnoldi[: self.E + 1,
self.M - self.N + 1 : self.E + 1]
if self.rho == np.inf:
_, s, V = np.linalg.svd(R, full_matrices = False)
else:
Rtower = np.zeros((R.shape[0] * (self.E - self.M), self.N + 1),
dtype = np.complex)
for k in range(self.E - self.M):
Rtower[k * R.shape[0] : (k + 1) * R.shape[0], :] = (
self.rho ** k
* R[:, self.M - self.N + 1 + k : self.M + 1 + k])
_, s, V = np.linalg.svd(Rtower, full_matrices = False)
return s[::-1], V.conj().T[:, ::-1]
def evalApprox(self, mu:complex) -> Np1D:
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Approximant as numpy complex vector.
"""
self.setupApprox()
powerlist = np.power(mu - self.mu0, 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.mu0
diff --git a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
index 064254e..8d727dd 100644
--- a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
+++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
@@ -1,307 +1,302 @@
# 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 .
#
import warnings
import numpy as np
from scipy.special import comb
import scipy.sparse.linalg as spla
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
+from rrompy.reduction_methods.base.pod_engine import PODEngine
from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng
from rrompy.utilities import purgeDict
__all__ = ['GenericApproximantTaylor']
class GenericApproximantTaylor(GenericApproximant):
"""
ROM single-point approximant computation for parametric problems
(ABSTRACT).
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 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.
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;
- '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;
- '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): What to plot of derivatives of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives 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.
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;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
RArnoldi: Right factor of derivatives QR factorization. If no Arnoldi,
set to None.
HArnoldi: Hessenberg factor of derivatives Arnoldi algorithm. If no
Arnoldi, set to None.
POD: Whether to compute QR factorization of derivatives.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
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.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
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.
"""
def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
approxParameters : DictAny = {}, plotDer : ListAny = [],
plotDerSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["E", "Emax", "sampleType"]
super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
approxParameters = approxParameters)
self.plotDer = plotDer
self.plotDerSpecs = plotDerSpecs
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
self.solDerivatives = None
if hasattr(self, "HArnoldi"):
del self.HArnoldi
if hasattr(self, "RArnoldi"):
del self.RArnoldi
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change E and Emax.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters")
approxParametersCopy = purgeDict(approxParameters,
["E", "Emax", "sampleType"],
True, True)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "E" in keyList:
self._E = approxParameters["E"]
self._approxParameters["E"] = self.E
if "Emax" in keyList:
self.Emax = approxParameters["Emax"]
else:
if not hasattr(self, "Emax"):
self.Emax = self.E
else:
self.Emax = self.Emax
else:
if "Emax" in keyList:
self._E = approxParameters["Emax"]
self._approxParameters["E"] = self.E
self.Emax = self.E
else:
if not (hasattr(self, "Emax") and hasattr(self, "E")):
raise Exception("At least one of E and Emax must be set.")
if "sampleType" in keyList:
self.sampleType = approxParameters["sampleType"]
elif hasattr(self, "sampleType"):
self.sampleType = self.sampleType
else:
self.sampleType = "KRYLOV"
@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
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warnings.warn("Prescribed E is too large. Updating Emax to E.",
stacklevel = 2)
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
self._Emax = Emax
if hasattr(self, "E") and self.Emax < self.E:
warnings.warn("Prescribed Emax is too small. Updating Emax to E.",
stacklevel = 2)
self.Emax = self.E
else:
self._approxParameters["Emax"] = self.Emax
if EmaxOld >= self.Emax and self.solDerivatives is not None:
self.solDerivatives = self.solDerivatives[:, : self.Emax + 1]
if hasattr(self, "HArnoldi"):
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 sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType
else: sampleTypeOld = -1
try:
sampleType = sampleType.upper().strip().replace(" ","")
if sampleType not in ["ARNOLDI", "KRYLOV"]: raise
if sampleType == "ARNOLDI" and self.HFEngine.Aterms > 2:
warnings.warn(("Arnoldi sampling not allowed for parameter "
"dependences with degree higher than 1. "
"Overriding to 'KRYLOV'."), stacklevel = 2)
sampleType = "KRYLOV"
self._sampleType = sampleType
except:
warnings.warn(("Prescribed sampleType not recognized. Overriding "
"to 'KRYLOV'."), stacklevel = 2)
self._sampleType = "KRYLOV"
self._approxParameters["sampleType"] = self.sampleType
if sampleTypeOld != self.sampleType:
self.resetSamples()
def computeDerivatives(self):
"""Compute derivatives of solution map starting from order 0."""
if self.solDerivatives is None:
lAs = len(self.As)
lbs = len(self.bs)
uOlds = None
for j in range(self.Emax + 1):
if j == 0:
self.HFEngine = self._HFEngine0
bs = self.bs
else:
bs = [np.zeros_like(uOlds[0])] * max(lAs - 1, lbs - j)
for ii in range(lbs - j):
bs[ii] = self.bs[j + ii]
for ii in range(lAs - 1):
for jj in range(1, min(j + 1, lAs - ii)):
bs[ii] =(bs[ii]
- comb(jj + ii, jj, exact = False)
* self.As[jj + ii].dot(uOlds[lAs - jj - 1]))
A, b = self.buildLS(self.mu0, bs = bs)
uj = spla.spsolve(A, b)
self.HSEngine.plot(uj, name = "u_{0}".format(j),
what = self.plotDer, **self.plotDerSpecs)
if j == 0:
uOlds = [np.zeros_like(uj)] * (lAs - 2)
else:
del uOlds[0]
uOlds = uOlds + [self.manageDerivatives(uj, j)]
def manageDerivatives(self, u:Np1D, pos:int) -> Np1D:
"""
Store derivatives of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
Returns:
Adjusted solution derivative as numpy complex vector.
"""
if pos == 0:
self.solDerivatives = np.empty((u.shape[0], self.Emax + 1),
dtype = np.complex)
self.solDerivatives[:, pos] = u
if self.sampleType == "ARNOLDI":
if pos == 0:
+ self.PODEngine = PODEngine(self.energyNormMatrix)
self.HArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
dtype = np.complex)
self.RArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
dtype = np.complex)
- self.HArnoldi[: pos, pos] = self.solDerivatives[:, pos].conj().dot(
- self.energyNormMatrix.dot(self.solDerivatives[:, : pos])
- ).conj()
- self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
- - self.solDerivatives[:, : pos].dot(
- self.HArnoldi[: pos, pos]))
- self.HArnoldi[pos, pos] = (self.solDerivatives[:, pos].conj().dot(
- self.energyNormMatrix.dot(self.solDerivatives[:, pos])))**.5
- self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
- / self.HArnoldi[pos, pos])
+ self.solDerivatives[:, pos], self.HArnoldi[: pos + 1, pos] = (
+ self.PODEngine.GS(self.solDerivatives[:, pos],
+ self.solDerivatives, pos))
if pos == 0:
self.RArnoldi[pos, pos] = self.HArnoldi[pos, pos]
else:
self.RArnoldi[:pos+1, pos] = self.HArnoldi[: pos+1, 1 : pos+1]\
.dot(self.RArnoldi[: pos, pos - 1])
return self.solDerivatives[:, pos]
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.solDerivatives is not None
and super().checkComputedApprox())