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())