diff --git a/.gitignore b/.gitignore index 3b3a0f1..fe8f754 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,8 @@ # Byte-compiled / optimized / DLL files __pycache__/ +/build/ +/dist/ +/RROMPy.egg-info/ *.py[cod] *.dat .ipynb_checkpoints/ diff --git a/main/__init__.py b/MANIFEST.in similarity index 100% rename from main/__init__.py rename to MANIFEST.in diff --git a/README.md b/README.md index f417887..2ae0ba7 100644 --- a/README.md +++ b/README.md @@ -1,27 +1,31 @@ -# HelmholtzMOR -Module for the solution and model order reduction of the Helmholtz problem. Coded in Python 3.6. +RROMPy +============ + +Module for the solution and rational model order reduction of parametric PDE-based problem. Coded in Python 3.6. + +## Prerequisites +**RBniCS** requires +* **numpy**. +* **scipy**. ## Installing -Cloning the repository is enough. If you want to use Fenics, you can install [Anaconda3/Miniconda3](http://anaconda.org/) and run the command (from the main directory ./) +Clone the repository ``` -conda env create --file conda-fenics.yml +git clone https://c4science.ch/source/RROMPy.git ``` - -This will create an environment where Fenics can be used. To activate the environment, it is sufficent to run (from any directory) +and install the package by typing ``` -source activate fenicsenv +python3 setup.py install ``` -To deactivate it, run +### Fenics +Most of the PDE engines already provided rely on [FEniCS](https://fenicsproject.org/). If you do not have FEniCS installed, you may want to create an [Anaconda3/Miniconda3](http://anaconda.org/) environment using the provided [file](./conda-fenics.yml) by running the command ``` -source deactivate +conda env create --file conda-fenics.yml ``` - -## Running the samples -Many examples can be found in the ./examples/ folder. -To run the Python examples you can use +This will create an environment where Fenics can be used. In order to use FEniCS, the environment must be activated through ``` -python examplename.py +source activate fenicsenv ``` ## License -This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](./LICENSE.md) file for details. +This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](./LICENSE) file for details. diff --git a/main/FenicsHelmholtzEngine.py b/main/FenicsHelmholtzEngine.py deleted file mode 100644 index e69c32b..0000000 --- a/main/FenicsHelmholtzEngine.py +++ /dev/null @@ -1,883 +0,0 @@ -#!/usr/bin/python - -from copy import copy -import warnings -import numpy as np -import sympy as sp -import sympy.printing as sprint -import scipy.sparse as scsp -import fenics as fen -from RROMPyTypes import DictAny, Np1D, Np2D, GenExpr, FenExpr, BfSExpr, Tuple, List - -fenZERO = fen.Constant(0) -fenZEROC = [fenZERO, fenZERO] - -class DomainError(Exception): - """ - Domain error exception. - - Args: - value: Human readable string describing the exception. - - Attributes: - value: Human readable string describing the exception. - """ - def __init__(self, value:str): - self.value = value - def __str__(self): - return self.value - -def CustomExpressionParser(value:GenExpr, degree:int)\ - -> Tuple[FenExpr, FenExpr]: - """ - Numerical and Fenics expressions parser. - - Args: - value: Expression to be parsed. Accepts 2-tuples composed of real and - imaginary parts. Available elementary formats are numbers, strings - and Fenics Expression's. - degree: Degree of Fenics FE interpolant of expression. - - Returns: - Fenics FE interpolant of expression. - """ - try: - if isinstance(value, (list, tuple)): - if len(value) == 1: - return CustomExpressionParser(value[0]) - elif len(value) != 2: - raise Exception("Parsing error") - if isinstance(value[0], (int, float, complex)): - if any([isinstance(y, complex) for y in value]): - raise Exception("Parsing error") - valueRe = fen.Constant(value[0]) - valueIm = fen.Constant(value[1]) - elif isinstance(value[0], str): - x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]", - real=True) - xyDict = {"x": x, "y": y, "z": z} - valueReParsed = value[0].replace("x[0]", "x")\ - .replace("x[1]", "y")\ - .replace("x[2]", "z") - valueImParsed = value[1].replace("x[0]", "x")\ - .replace("x[1]", "y")\ - .replace("x[2]", "z") - valueReSym = sp.sympify(valueReParsed, locals=xyDict)\ - .subs([(x, x0), (y, x1), (z, x2)]) - valueImSym = sp.sympify(valueImParsed, locals=xyDict)\ - .subs([(x, x0), (y, x1), (z, x2)]) - valueReStr = sprint.ccode(valueReSym).rpartition("\n")[-1] - valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1] - valueRe = fen.Expression(valueReStr, degree = degree) - valueIm = fen.Expression(valueImStr, degree = degree) - else: - valueRe = value[0] - valueIm = value[1] - else: - if isinstance(value, (int, float, complex)): - valueRe = fen.Constant(np.real(value)) - valueIm = fen.Constant(np.imag(value)) - elif isinstance(value, str): - x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]", - real=True) - xyDict = {"x": x, "y": y, "z": z} - valueParsed = value.replace("x[0]", "x").replace("x[1]", "y")\ - .replace("x[2]", "z") - valueSym = sp.sympify(valueParsed, locals=xyDict)\ - .subs([(x, x0), (y, x1), (z, x2)]) - valueReStr = sprint.ccode(sp.re(valueSym)).rpartition("\n")[-1] - valueImStr = sprint.ccode(sp.im(valueSym)).rpartition("\n")[-1] - valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1] - valueRe = fen.Expression(valueReStr, degree = degree) - valueIm = fen.Expression(valueImStr, degree = degree) - else: - valueRe = value - valueIm = fenZERO - except: - raise Exception("Parsing error") - return valueRe, valueIm - -class FenicsHelmholtzEngine: - """ - Fenics-based solver for Helmholtz problems. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu + h u = g2 on \Gamma_R - - Args: - mesh: Domain of Helmholtz problem. - FEDegree: FE degree. - wavenumber2: Value of k^2. - diffusivity(optional): Value of a. Defaults to identically 1. - refractionIndex(optional): Value of n. Defaults to identically 1. - forcingTerm(optional): Value of f. Defaults to identically 0. - DirichletBoundary(optional): Function flagging Dirichlet boundary as - True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are - accepted. Defaults to False everywhere. - NeumannBoundary(optional): Function flagging Neumann boundary as True, - in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - RobinBoundary(optional): Function flagging Robin boundary as True, in - Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - DirichletDatum(optional): Value of u0. Defaults to identically 0. - NeumannDatum(optional): Value of g1. Defaults to identically 0. - RobinDatum(optional): Value of h. Defaults to identically 0. - RobinDatum2(optional): Value of g2. Defaults to identically 0. - - Attributes: - FEDegree: FE degree. - wavenumber2: Copy of processed wavenumber2 parameter. - diffusivity: Copy of processed diffusivity parameter. - refractionIndex: Copy of processed refractionIndex parameter. - forcingTerm: Copy of processed forcingTerm parameter. - DirichletBoundary: Copy of processed DirichletBoundary parameter. - NeumannBoundary: Copy of processed NeumannBoundary parameter. - RobinBoundary: Copy of processed RobinBoundary parameter. - DirichletDatum: Copy of processed DirichletDatum parameter. - NeumannDatum: Copy of processed NeumannDatum parameter. - RobinDatum: Copy of processed RobinDatum parameter. - RobinDatum2: Copy of processed RobinDatum2 parameter. - LU: Whether to use LU solver for computation of derivatives. - V: Real FE space. - u: Helmholtz problem solution as numpy complex vector. - v1: Generic trial functions for variational form evaluation. - v2: Generic test functions for variational form evaluation. - ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). - nRe,nIm: Real and imaginary parts of n. - n2Re,n2Im: Real and imaginary parts of n^2. - fRe,fIm: Real and imaginary parts of f. - u0Re,u0Im: Real and imaginary parts of u0. - g1Re,g1Im: Real and imaginary parts of g1. - hRe,hIm: Real and imaginary parts of h. - g2Re,g2Im: Real and imaginary parts of g2. - DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary - parts of Dirichlet data. - A*: Scipy sparse array representation (in CSC format) of A*. - b*: Numpy array representation of b*. - """ - - def __init__(self, mesh:fen.mesh, - FEDegree:int, - wavenumber:complex, - diffusivity : GenExpr = 1, - refractionIndex : GenExpr = 1, - forcingTerm : GenExpr = fenZEROC, - DirichletBoundary : BfSExpr = None, - NeumannBoundary : BfSExpr = None, - RobinBoundary : BfSExpr = None, - DirichletDatum : GenExpr = fenZEROC, - NeumannDatum : GenExpr = fenZEROC, - RobinDatum : GenExpr = fenZEROC, - RobinDatum2 : GenExpr = fenZEROC): - self.FEDegree = FEDegree - - # process boundaries - boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary] - for i in range(3): - if isinstance(boundariesList[i], str): - boundariesList[i] = boundariesList[i].upper() - if boundariesList[i] == "NONE": boundariesList[i] = None - DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList - if boundariesList.count(None) == 3: - raise DomainError("At least one boundary must be prescribed.") - if boundariesList.count("ALL") + boundariesList.count("REST") >= 2: - raise DomainError("Only one keyword 'ALL'/'REST' can be used.") - def DirichletB(x, on_boundary): - if DirichletBoundary is None: - return False - elif DirichletBoundary == "ALL": - return on_boundary - elif DirichletBoundary == "REST": - return (on_boundary - and not self.NeumannBoundary(x, on_boundary) - and not self.RobinBoundary(x, on_boundary)) - else: - return DirichletBoundary(x, on_boundary) - def NeumannB(x, on_boundary): - if NeumannBoundary is None: - return False - elif NeumannBoundary == "ALL": - return on_boundary - elif NeumannBoundary == "REST": - return (on_boundary - and not self.DirichletBoundary(x, on_boundary) - and not self.RobinBoundary(x, on_boundary)) - else: - return NeumannBoundary(x, on_boundary) - def RobinB(x, on_boundary): - if RobinBoundary is None: - return False - elif RobinBoundary == "ALL": - return on_boundary - elif RobinBoundary == "REST": - return (on_boundary - and not self.DirichletBoundary(x, on_boundary) - and not self.NeumannBoundary(x, on_boundary)) - else: - return RobinBoundary(x, on_boundary) - self.DirichletBoundary = DirichletB - self.NeumannBoundary = NeumannB - self.RobinBoundary = RobinB - - # create boundary measure - class NBoundary(fen.SubDomain): - def inside(self, x, on_boundary): - return NeumannB(x, on_boundary) - class RBoundary(fen.SubDomain): - def inside(self, x, on_boundary): - return RobinB(x, on_boundary) - boundary_markers = fen.MeshFunction("size_t", mesh, - mesh.topology().dim() - 1) - NBoundary().mark(boundary_markers, 0) - RBoundary().mark(boundary_markers, 1) - self.ds = fen.Measure("ds", domain=mesh, - subdomain_data=boundary_markers) - - self.V = fen.FunctionSpace(mesh, "P", FEDegree) - self.v1 = fen.TrialFunction(self.V) - self.v2 = fen.TestFunction(self.V) - self.wavenumber2 = wavenumber ** 2 - self.diffusivity = diffusivity - self.refractionIndex = refractionIndex - self.forcingTerm = forcingTerm - self.DirichletDatum = DirichletDatum - self.NeumannDatum = NeumannDatum - self.RobinDatum = RobinDatum - self.RobinDatum2 = RobinDatum2 - - functional = lambda self, u: 0. - - def Vdim(self) -> int: - """Dimension of finite element space.""" - return self.V.dim() - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - normMatFen = fen.assemble(( - fen.dot(fen.grad(self.v1), fen.grad(self.v2)) - + w * fen.dot(self.v1, self.v2) - ) * fen.dx) - normMat = fen.as_backend_type(normMatFen).mat() - return scsp.csr_matrix(normMat.getValuesCSR()[::-1], - shape = normMat.size) - - def problemData(self) -> DictAny: - """List of HF problem data.""" - dataDict = {} - dataDict["forcingTerm"] = self.forcingTerm - dataDict["DirichletDatum"] = self.DirichletDatum - dataDict["NeumannDatum"] = self.NeumannDatum - dataDict["RobinDatum2"] = self.RobinDatum2 - return dataDict - - def setProblemData(self, dataDict:DictAny, k2:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k2: Parameter value. - """ - self.wavenumber2 = k2 - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "b0"): del self.b0 - - def liftDirichletData(self) -> Np1D: - """Lift Dirichlet datum.""" - solLRe = fen.interpolate(self.u0Re, self.V) - solLIm = fen.interpolate(self.u0Im, self.V) - return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector()) - - @property - def diffusivity(self): - """Value of a. Its assignment changes aRe and aIm.""" - return self._diffusivity - @diffusivity.setter - def diffusivity(self, diffusivity): - if hasattr(self, "A0"): del self.A0 - self.aRe, self.aIm = CustomExpressionParser(diffusivity, - degree = self.FEDegree) - self._diffusivity = [self.aRe, self.aIm] - - @property - def refractionIndex(self): - """Value of n. Its assignment changes nRe, nIm, n2Re and n2Im.""" - return self._refractionIndex - @refractionIndex.setter - def refractionIndex(self, refractionIndex): - if hasattr(self, "A1"): del self.A1 - self.nRe, self.nIm = CustomExpressionParser(refractionIndex, - degree = self.FEDegree) - self.n2Re = self.nRe*self.nRe - self.nIm*self.nIm - self.n2Im = 2 * self.nRe * self.nIm - self._refractionIndex = [self.nRe, self.nIm] - - @property - def forcingTerm(self): - """Value of f. Its assignment changes fRe and fIm.""" - return self._forcingTerm - @forcingTerm.setter - def forcingTerm(self, forcingTerm): - if hasattr(self, "b0"): del self.b0 - self.fRe, self.fIm = CustomExpressionParser(forcingTerm, - degree = self.FEDegree) - self._forcingTerm = [self.fRe, self.fIm] - - @property - def DirichletDatum(self): - """ - Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and - DirichletBCIm. - """ - return self._DirichletDatum - @DirichletDatum.setter - def DirichletDatum(self, DirichletDatum): - if hasattr(self, "b0"): del self.b0 - self.u0Re, self.u0Im = CustomExpressionParser(DirichletDatum, - degree = self.FEDegree) - self.DirichletBCRe = fen.DirichletBC(self.V, self.u0Re, - self.DirichletBoundary) - self.DirichletBCIm = fen.DirichletBC(self.V, self.u0Im, - self.DirichletBoundary) - self.DirichletBC0 = fen.DirichletBC(self.V, fenZERO, - self.DirichletBoundary) - self._DirichletDatum = [self.u0Re, self.u0Im] - - @property - def NeumannDatum(self): - """Value of g1. Its assignment changes g1Re and g1Im.""" - return self._NeumannDatum - @NeumannDatum.setter - def NeumannDatum(self, NeumannDatum): - if hasattr(self, "b0"): del self.b0 - self.g1Re, self.g1Im = CustomExpressionParser(NeumannDatum, - degree = self.FEDegree) - self._NeumannDatum = [self.g1Re, self.g1Im] - - @property - def RobinDatum(self): - """Value of h. Its assignment changes hRe and hIm.""" - return self._RobinDatum - @RobinDatum.setter - def RobinDatum(self, RobinDatum): - if hasattr(self, "A0"): del self.A0 - self.hRe, self.hIm = CustomExpressionParser(RobinDatum, - degree = self.FEDegree) - self._RobinDatum = [self.hRe, self.hIm] - - @property - def RobinDatum2(self): - """Value of g2. Its assignment changes g2Re and g2Im.""" - return self._RobinDatum2 - @RobinDatum2.setter - def RobinDatum2(self, RobinDatum2): - if hasattr(self, "b0"): del self.b0 - self.g2Re, self.g2Im = CustomExpressionParser(RobinDatum2, - degree = self.FEDegree) - self._RobinDatum2 = [self.g2Re, self.g2Im] - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - if not hasattr(self, "A0"): - a0Re = (self.aRe * fen.dot(fen.grad(self.v1), - fen.grad(self.v2)) * fen.dx - + self.hRe * fen.dot(self.v1, self.v2) * self.ds(1)) - a0Im = (self.aIm * fen.dot(fen.grad(self.v1), - fen.grad(self.v2)) * fen.dx - + self.hIm * fen.dot(self.v1, self.v2) * self.ds(1)) - A0Re = fen.assemble(a0Re) - A0Im = fen.assemble(a0Im) - self.DirichletBC0.apply(A0Re) - self.DirichletBC0.zero(A0Im) - A0ReMat = fen.as_backend_type(A0Re).mat() - A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR() - A0Imr, A0Imc, A0Imv = fen.as_backend_type( - A0Im).mat().getValuesCSR() - self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), - shape = A0ReMat.size) - + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), - shape = A0ReMat.size)) - if not hasattr(self, "A1"): - a1Re = - self.n2Re * fen.dot(self.v1, self.v2) * fen.dx - a1Im = - self.n2Im * fen.dot(self.v1, self.v2) * fen.dx - A1Re = fen.assemble(a1Re) - A1Im = fen.assemble(a1Im) - self.DirichletBC0.zero(A1Re) - self.DirichletBC0.zero(A1Im) - A1Rer, A1Rec, A1Rev = fen.as_backend_type(A1Re).mat()\ - .getValuesCSR() - A1Imr, A1Imc, A1Imv = fen.as_backend_type(A1Im).mat()\ - .getValuesCSR() - self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer), - shape = self.A0.shape) - + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr), - shape = self.A0.shape)) - - def assembleb(self): - """Assemble RHS blocks of linear system.""" - if not hasattr(self, "b0"): - L0Re = (fen.dot(self.fRe, self.v2) * fen.dx - + fen.dot(self.g1Re, self.v2) * self.ds(0) - + fen.dot(self.g2Re, self.v2) * self.ds(1)) - L0Im = (fen.dot(self.fIm, self.v2) * fen.dx - + fen.dot(self.g1Im, self.v2) * self.ds(0) - + fen.dot(self.g2Im, self.v2) * self.ds(1)) - b0Re = fen.assemble(L0Re) - b0Im = fen.assemble(L0Im) - self.DirichletBCRe.apply(b0Re) - self.DirichletBCIm.apply(b0Im) - self.b0 = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) - - parDegree = 1 - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1] - - def bs(self) -> List[Np1D]: - """RHS blocks of linear system.""" - self.assembleb() - return [self.b0] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber2 * self.A1 - - def b(self) -> Np1D: - """Assemble RHS of linear system.""" - self.assembleb() - return self.b0 - - def solve(self) -> Np1D: - """ - Find solution of linear system. - - Returns: - Solution vector. - """ - return scsp.linalg.spsolve(self.A(), self.b()) - - def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]: - """ - Get blocks of linear system. - - Returns: - Blocks of LS (\sum_{j=0}^J k^j A_j = \sum_{j=0}^K k^j b_j). - """ - return self.As(), self.bs() - - -class FenicsHelmholtzScatteringEngine(FenicsHelmholtzEngine): - """ - Fenics-based solver for Helmholtz scattering problems with Robin ABC. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu - i k u = g2 on \Gamma_R - - Args: - mesh: Domain of Helmholtz problem. - FEDegree: FE degree. - wavenumber: Value of k. - diffusivity(optional): Value of a. Defaults to identically 1. - refractionIndex(optional): Value of n. Defaults to identically 1. - forcingTerm(optional): Value of f. Defaults to identically 0. - DirichletBoundary(optional): Function flagging Dirichlet boundary as - True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are - accepted. Defaults to False everywhere. - NeumannBoundary(optional): Function flagging Neumann boundary as True, - in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - RobinBoundary(optional): Function flagging Robin boundary as True, in - Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - DirichletDatum(optional): Value of u0. Defaults to identically 0. - NeumannDatum(optional): Value of g1. Defaults to identically 0. - RobinDatum2(optional): Value of g2. Defaults to identically 0. - - Attributes: - FEDegree: FE degree. - wavenumber: Copy of processed wavenumber parameter. - diffusivity: Copy of processed diffusivity parameter. - refractionIndex: Copy of processed refractionIndex parameter. - forcingTerm: Copy of processed forcingTerm parameter. - DirichletBoundary: Copy of processed DirichletBoundary parameter. - NeumannBoundary: Copy of processed NeumannBoundary parameter. - RobinBoundary: Copy of processed RobinBoundary parameter. - DirichletDatum: Copy of processed DirichletDatum parameter. - NeumannDatum: Copy of processed NeumannDatum parameter. - RobinDatum: Value of h, i.e. (- i * k). - RobinDatum2: Copy of processed RobinDatum2 parameter. - V: Real FE space. - u: Helmholtz problem solution as numpy complex vector. - v1: Generic trial functions for variational form evaluation. - v2: Generic test functions for variational form evaluation. - ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). - nRe,nIm: Real and imaginary parts of n. - n2Re,n2Im: Real and imaginary parts of n^2. - fRe,fIm: Real and imaginary parts of f. - u0Re,u0Im: Real and imaginary parts of u0. - g1Re,g1Im: Real and imaginary parts of g1. - hRe,hIm: Real and imaginary parts of h. - h: Value of h. - g2Re,g2Im: Real and imaginary parts of g2. - DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary - parts of Dirichlet data. - A*: Scipy sparse array representation (in CSC format) of A*. - b*: Numpy array representation of b*. - """ - - def __init__(self, mesh:fen.mesh, - FEDegree:int, - wavenumber:complex, - diffusivity : GenExpr = 1, - refractionIndex : GenExpr = 1, - forcingTerm : GenExpr = fenZEROC, - DirichletBoundary : BfSExpr = None, - NeumannBoundary : BfSExpr = None, - RobinBoundary : BfSExpr = None, - DirichletDatum : GenExpr = fenZEROC, - NeumannDatum : GenExpr = fenZEROC, - RobinDatum2 : GenExpr = fenZEROC): - self.wavenumber = wavenumber - FenicsHelmholtzEngine.__init__(self, mesh = mesh, - FEDegree = FEDegree, - wavenumber = wavenumber, - diffusivity = diffusivity, - refractionIndex = refractionIndex, - forcingTerm = forcingTerm, - DirichletBoundary = DirichletBoundary, - NeumannBoundary = NeumannBoundary, - RobinBoundary = RobinBoundary, - DirichletDatum = DirichletDatum, - NeumannDatum = NeumannDatum, - RobinDatum = -1.j * wavenumber, - RobinDatum2 = RobinDatum2) - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - return FenicsHelmholtzEngine.energyNormMatrix(self, w ** 2.) - - def setProblemData(self, dataDict:DictAny, k:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k: Parameter value. - """ - self.wavenumber = k - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "A2"): del self.A2 - if hasattr(self, "b0"): del self.b0 - - @property - def wavenumber2(self): - """Value of k^2.""" - return self._wavenumber2 - @wavenumber2.setter - def wavenumber2(self, wavenumber2): - self._wavenumber2 = wavenumber2 - self._wavenumber = self.wavenumber2 ** .5 - if (not hasattr(self, "h") - or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)): - self.RobinDatum = -1.j * self.wavenumber - - @property - def wavenumber(self): - """Value of k.""" - return self._wavenumber - @wavenumber.setter - def wavenumber(self, wavenumber): - self.wavenumber2 = wavenumber ** 2. - - @property - def RobinDatum(self): - """ - Value of h. Its assignment changes hRe and hIm (and maybe kRe, kIm, k - and z). - """ - return super().RobinDatum - @RobinDatum.setter - def RobinDatum(self, RobinDatum): - self.h = RobinDatum - self.hRe, self.hIm = np.real(self.h), np.imag(self.h) - self._RobinDatum = [self.hRe, self.hIm] - if (not hasattr(self, "wavenumber") - or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)): - self.wavenumber = 1.j * self.h - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - if not hasattr(self, "A0"): - a0Re = self.aRe * fen.dot(fen.grad(self.v1), - fen.grad(self.v2)) * fen.dx - a0Im = self.aIm * fen.dot(fen.grad(self.v1), - fen.grad(self.v2)) * fen.dx - A0Re = fen.assemble(a0Re) - A0Im = fen.assemble(a0Im) - self.DirichletBC0.apply(A0Re) - self.DirichletBC0.zero(A0Im) - A0ReMat = fen.as_backend_type(A0Re).mat() - A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR() - A0Imr, A0Imc, A0Imv = fen.as_backend_type( - A0Im).mat().getValuesCSR() - self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), - shape = A0ReMat.size) - + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), - shape = A0ReMat.size)) - if hasattr(self, "A1"): - aux = copy(self.A1) - else: - aux = None - if not hasattr(self, "A2"): - if aux is not None: - del self.A1 - FenicsHelmholtzEngine.assembleA(self) - self.A2 = copy(self.A1) - if aux is None: - a1 = fen.dot(self.v1, self.v2) * self.ds(1) - A1 = fen.assemble(a1) - self.DirichletBC0.zero(A1) - A1r, A1c, A1v = fen.as_backend_type(A1).mat().getValuesCSR() - self.A1 = - 1.j * scsp.csr_matrix((A1v, A1c, A1r), - shape = self.A0.shape) - else: - self.A1 = aux - - parDegree = 2 - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1, self.A2] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber * self.A1 + self.wavenumber2 * self.A2 - - -class FenicsHelmholtzScatteringAugmentedEngine(FenicsHelmholtzScatteringEngine): - """ - Fenics-based solver for Helmholtz scattering problems with Robin ABC. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu - i k u = g2 on \Gamma_R - Linear dependence on k is achieved by introducing an additional variable. - - Args: - mesh: Domain of Helmholtz problem. - FEDegree: FE degree. - wavenumber: Value of k. - diffusivity(optional): Value of a. Defaults to identically 1. - refractionIndex(optional): Value of n. Defaults to identically 1. - forcingTerm(optional): Value of f. Defaults to identically 0. - DirichletBoundary(optional): Function flagging Dirichlet boundary as - True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are - accepted. Defaults to False everywhere. - NeumannBoundary(optional): Function flagging Neumann boundary as True, - in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - RobinBoundary(optional): Function flagging Robin boundary as True, in - Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted. - Defaults to False everywhere. - DirichletDatum(optional): Value of u0. Defaults to identically 0. - NeumannDatum(optional): Value of g1. Defaults to identically 0. - RobinDatum2(optional): Value of g2. Defaults to identically 0. - constraintType(optional): Type of augmentation. Keywords 'IDENTITY' and - 'MASS' are accepted. Defaults to 'IDENTITY'. - - Attributes: - FEDegree: FE degree. - wavenumber: Copy of processed wavenumber parameter. - diffusivity: Copy of processed diffusivity parameter. - refractionIndex: Copy of processed refractionIndex parameter. - forcingTerm: Copy of processed forcingTerm parameter. - DirichletBoundary: Copy of processed DirichletBoundary parameter. - NeumannBoundary: Copy of processed NeumannBoundary parameter. - RobinBoundary: Copy of processed RobinBoundary parameter. - DirichletDatum: Copy of processed DirichletDatum parameter. - NeumannDatum: Copy of processed NeumannDatum parameter. - RobinDatum: Value of h, i.e. (- i * k). - RobinDatum2: Copy of processed RobinDatum2 parameter. - constraintType: Type of augmentation. - V: Real FE space. - u: Helmholtz problem solution as numpy complex vector. - v1: Generic trial functions for variational form evaluation. - v2: Generic test functions for variational form evaluation. - ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). - nRe,nIm: Real and imaginary parts of n. - n2Re,n2Im: Real and imaginary parts of n^2. - fRe,fIm: Real and imaginary parts of f. - u0Re,u0Im: Real and imaginary parts of u0. - g1Re,g1Im: Real and imaginary parts of g1. - hRe,hIm: Real and imaginary parts of h. - g2Re,g2Im: Real and imaginary parts of g2. - DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary - parts of Dirichlet data. - A*: Scipy sparse array representation (in CSC format) of A*. - b*: Numpy array representation of b*. - """ - - def __init__(self, mesh:fen.mesh, - FEDegree:int, - wavenumber:complex, - diffusivity : GenExpr = 1, - refractionIndex : GenExpr = 1, - forcingTerm : GenExpr = fenZEROC, - DirichletBoundary : BfSExpr = None, - NeumannBoundary : BfSExpr = None, - RobinBoundary : BfSExpr = None, - DirichletDatum : GenExpr = fenZEROC, - NeumannDatum : GenExpr = fenZEROC, - RobinDatum2 : GenExpr = fenZEROC, - constraintType : str = "IDENTITY"): - self.constraintType = constraintType - FenicsHelmholtzScatteringEngine.__init__(self, mesh = mesh, - FEDegree = FEDegree, - wavenumber = wavenumber, - diffusivity = diffusivity, - refractionIndex = refractionIndex, - forcingTerm = forcingTerm, - DirichletBoundary = DirichletBoundary, - NeumannBoundary = NeumannBoundary, - RobinBoundary = RobinBoundary, - DirichletDatum = DirichletDatum, - NeumannDatum = NeumannDatum, - RobinDatum2 = RobinDatum2) - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - normMatFen = FenicsHelmholtzScatteringEngine.energyNormMatrix(self, w) - return scsp.block_diag((normMatFen, 1/w * normMatFen)) - - def liftDirichletData(self): - """Lift Dirichlet datum.""" - lift1 = FenicsHelmholtzScatteringEngine.liftDirichletData(self) - return np.pad(lift1, (0, len(lift1)), 'constant') - - def setProblemData(self, dataDict:dict, k:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k: Parameter value. - """ - self.wavenumber = k - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "b0"): del self.b0 - - @property - def constraintType(self): - """Value of constraintType.""" - return self._constraintType - @constraintType.setter - def constraintType(self, constraintType): - try: - constraintType = constraintType.upper().strip().replace(" ","") - if constraintType not in ["IDENTITY", "MASS"]: raise - self._constraintType = constraintType - except: - warnings.warn(("Prescribed constraintType not recognized. " - "Overriding to 'IDENTITY'."), stacklevel = 2) - self.constraintType = "IDENTITY" - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - builtA0 = hasattr(self, "A0") - builtA1 = hasattr(self, "A1") - if not (builtA0 and builtA1): - if builtA1 and self.constraintType == "IDENTITY": - self.A2 = None - FenicsHelmholtzScatteringEngine.assembleA(self) - if self.constraintType == "IDENTITY": - if not builtA0: - block = scsp.eye(self.A0.shape[0]) - else: - block = scsp.eye(self.A1.shape[0]) - else: #MASS - block = - self.A2 - I = list(self.DirichletBC0.get_boundary_values().keys()) - warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning) - block[I, I] = 1. - warnings.simplefilter("default", scsp.SparseEfficiencyWarning) - if not builtA0: - self.A0 = scsp.block_diag((self.A0, block)) - if not builtA1: - self.A1 = scsp.bmat([[self.A1, self.A2], [- block, None]]) - if hasattr(self, 'A2'): del self.A2 - - def assembleb(self): - """Assemble RHS of linear system.""" - if not hasattr(self, "b0"): - FenicsHelmholtzScatteringEngine.assembleb(self) - self.b0 = np.pad(self.b0, (0, self.b0.shape[0]), 'constant') - - parDegree = 1 - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber * self.A1 - - def solve(self) -> Tuple[Np1D, Np1D]: - """ - Find solution of linear system. - - Returns: - Solution vectors. - """ - u = FenicsHelmholtzScatteringEngine.solve(self) - lu2 = int(len(u) / 2) - return np.array([u[: lu2], u[lu2 :]]) diff --git a/main/FreeFemConversionTools.py b/main/FreeFemConversionTools.py deleted file mode 100644 index bcde571..0000000 --- a/main/FreeFemConversionTools.py +++ /dev/null @@ -1,198 +0,0 @@ -#!/usr/bin/python - -import numpy as np -import scipy.sparse as sp -import utilities -from RROMPyTypes import Np1D, Np2D, Tuple, List - -HSVjet = ("real[int] HSVjet=[0.6667,1.0000,0.5625,0.6667,1.0000,0.6250,0.6667," -"1.0000,0.6875,0.6667,1.0000,0.7500,0.6667,1.0000,0.8125,0.6667,1.0000,0.8750," -"0.6667,1.0000,0.9375,0.6667,1.0000,1.0000,0.6563,1.0000,1.0000,0.6458,1.0000," -"1.0000,0.6354,1.0000,1.0000,0.6250,1.0000,1.0000,0.6146,1.0000,1.0000,0.6042," -"1.0000,1.0000,0.5938,1.0000,1.0000,0.5833,1.0000,1.0000,0.5729,1.0000,1.0000," -"0.5625,1.0000,1.0000,0.5521,1.0000,1.0000,0.5417,1.0000,1.0000,0.5313,1.0000," -"1.0000,0.5208,1.0000,1.0000,0.5104,1.0000,1.0000,0.5000,1.0000,1.0000,0.4889," -"0.9375,1.0000,0.4762,0.8750,1.0000,0.4615,0.8125,1.0000,0.4444,0.7500,1.0000," -"0.4242,0.6875,1.0000,0.4000,0.6250,1.0000,0.3704,0.5625,1.0000,0.3333,0.5000," -"1.0000,0.2963,0.5625,1.0000,0.2667,0.6250,1.0000,0.2424,0.6875,1.0000,0.2222," -"0.7500,1.0000,0.2051,0.8125,1.0000,0.1905,0.8750,1.0000,0.1778,0.9375,1.0000," -"0.1667,1.0000,1.0000,0.1563,1.0000,1.0000,0.1458,1.0000,1.0000,0.1354,1.0000," -"1.0000,0.1250,1.0000,1.0000,0.1146,1.0000,1.0000,0.1042,1.0000,1.0000,0.0938," -"1.0000,1.0000,0.0833,1.0000,1.0000,0.0729,1.0000,1.0000,0.0625,1.0000,1.0000," -"0.0521,1.0000,1.0000,0.0417,1.0000,1.0000,0.0313,1.0000,1.0000,0.0208,1.0000," -"1.0000,0.0104,1.0000,1.0000,0,1.0000,1.0000,0,1.0000,0.9375,0,1.0000,0.8750," -"0,1.0000,0.8125,0,1.0000,0.7500,0,1.0000,0.6875,0,1.0000,0.6250,0,1.0000," -"0.5625,0,1.0000,0.5000];") - - -def prod2(v1 : str = "u", v2 : str = "v", compl : bool = True) -> str: - if compl: - return "({} * conj({}))".format(v1, v2) - else: - return "({} * {})".format(v1, v2) - -def diff(var : str = "u", dim : int = 2) -> str: - dstr = "(abs(dx({}))^2".format(var) - if dim > 1: - dstr = dstr + " + abs(dy({}))^2".format(var) - if dim > 2: - dstr = dstr + " + abs(dz({}))^2".format(var) - return dstr + ")" - -def diff2(v1 : str = "u", v2 : str = "v", dim : int = 2, - compl : bool = True) -> str: - dstr = prod2("(dx({})".format(v1), "dx({})".format(v2), compl) - if dim > 1: - dstr = dstr + " + " + prod2("dy({})".format(v1), - "dy({})".format(v2), compl) - if dim > 2: - dstr = dstr + " + " + prod2("dz({})".format(v1), - "dz({})".format(v2), compl) - return dstr + ")" - -def ffvec2npvec(filename:str) -> Np1D: - with open(filename) as f: - a = f.readline() - v = np.zeros(int(a), dtype = np.complex) - i = 0 - for line in f: - for x in line[: -1].split(): - if "(" in x: - y = x[1:-1].split(",") - v[i] = float(y[0]) + 1.j * float(y[1]) - else: - v[i] = complex(x) - i = i + 1 - return v - -def ffvec2npveccleanrows(filename:str, tgv : float = 1e30, - DBC : List[int] = []) -> Np1D: - jBC = 0 - with open(filename) as f: - a = f.readline() - v = np.zeros(int(a), dtype = np.complex) - i = 0 - for line in f: - for x in line[: -1].split(): - if "(" in x: - y = x[1:-1].split(",") - v[i] = float(y[0]) + 1.j * float(y[1]) - else: - v[i] = complex(x) - if i in DBC[jBC : jBC + 2]: - v[i] = v[i] / tgv - if jBC + 1 < len(DBC) and i == DBC[jBC + 1]: - jBC = jBC + 1 - i = i + 1 - return v - -def ffmat2spmat(filename:str) -> Np2D: - with open(filename) as f: - a = f.readline()[: -1] - while a[0] == "#": - a = f.readline()[: -1] - n, m, sym, nnz = [int(x) for x in a.split()] - row = np.zeros(nnz) - col = np.zeros(nnz) - data = np.zeros(nnz, dtype = np.complex) - for i, line in enumerate(f): - x = line[: -1].split() - row[i] = int(x[0]) - 1 - col[i] = int(x[1]) - 1 - if "(" in x[2]: - y = x[2][1:-1].split(",") - data[i] = float(y[0]) + 1.j * float(y[1]) - else: - data[i] = complex(x) - return sp.coo_matrix((data, (row, col)), shape = (n, m)) - -def ffmat2spmatnormalize(filename:str, tgv : float = 1e30)\ - -> Tuple[Np2D, List[int]]: - with open(filename) as f: - a = f.readline()[: -1] - while a[0] == "#": - a = f.readline()[: -1] - n, m, sym, nnz = [int(x) for x in a.split()] - DBC = [] - row = np.zeros(nnz) - col = np.zeros(nnz) - data = np.zeros(nnz, dtype = np.complex) - for i, line in enumerate(f): - x = line[: -1].split() - row[i] = int(x[0]) - 1 - col[i] = int(x[1]) - 1 - if "(" in x[2]: - y = x[2][1:-1].split(",") - data[i] = float(y[0]) + 1.j * float(y[1]) - else: - data[i] = complex(x) - if row[i] == col[i] and np.isclose(data[i], tgv): - DBC.append(row[i]) - jBC = 0 - idxkeep = [] - for i in range(nnz): - if row[i] in DBC[jBC : jBC + 2]: - if jBC + 1 < len(DBC) and row[i] == DBC[jBC + 1]: - jBC = jBC + 1 - if row[i] != col[i]: continue - data[i] = 1. - idxkeep.append(i) - row, col, data = row[idxkeep], col[idxkeep], data[idxkeep] - return sp.coo_matrix((data, (row, col)), shape = (n, m)), DBC - -def ffmat2spmatcleanrows(filename:str, DBC:List[int]) -> Np2D: - jBC = 0 - with open(filename) as f: - a = f.readline()[: -1] - while a[0] == "#": - a = f.readline()[: -1] - n, m, sym, nnz = [int(x) for x in a.split()] - idxkeep = [] - row = np.zeros(nnz) - col = np.zeros(nnz) - data = np.zeros(nnz, dtype = np.complex) - for i, line in enumerate(f): - x = line[: -1].split() - row[i] = int(x[0]) - 1 - if row[i] in DBC[jBC : jBC + 2]: - if jBC + 1 < len(DBC) and row[i] == DBC[jBC + 1]: - jBC = jBC + 1 - continue - col[i] = int(x[1]) - 1 - if "(" in x[2]: - y = x[2][1:-1].split(",") - data[i] = float(y[0]) + 1.j * float(y[1]) - else: - data[i] = complex(x) - idxkeep.append(i) - row, col, data = row[idxkeep], col[idxkeep], data[idxkeep] - return sp.coo_matrix((data, (row, col)), shape = (n, m)) - -def npvec2ffvec(v:Np1D) -> str: - filename = utilities.getNewFilename("v", "tmp") - with open(filename, 'w', buffering = 1) as f: - f.write(str(len(v)) + '\n') - for x in v: - if v.dtype == np.complex: - f.write("({0},{1})\n".format(np.real(x), np.imag(x))) - else: - f.write(str(x) + '\n') - return filename - -def spmat2ffmat(m:Np2D) -> str: - mEff = m.tocoo() - filename = utilities.getNewFilename("m", "tmp") - with open(filename, 'w', buffering = 1) as f: - f.write(("# Sparse Matrix (Morse)\n# first line: n m (is symmetic) " - "nbcoef\n# after for each nonzero coefficient: i j a_ij " - "where (i,j) \in {1,...,n}x{1,...,m}\n")) - f.write("{} {} 0 {}\n".format(mEff.shape[0], mEff.shape[1], mEff.nnz)) - for i in range(mEff.nnz): - r = mEff.row[i] - c = mEff.col[i] - x = mEff.data[i] - if mEff.dtype == np.complex: - f.write("{} {} ({},{})\n".format(r, c, np.real(x), np.imag(x))) - else: - f.write("{} {} {}\n".format(r, c, x)) - return filename diff --git a/main/FreeFemHSEngine.py b/main/FreeFemHSEngine.py deleted file mode 100644 index 15abeab..0000000 --- a/main/FreeFemHSEngine.py +++ /dev/null @@ -1,327 +0,0 @@ -#!/usr/bin/python - -import os -import subprocess -import numpy as np -import warnings -import utilities -import FreeFemConversionTools as FFCT -from RROMPyTypes import Np1D, strLst, N2FSExpr - -class FreeFemHSEngine: - """ - FreeFem++-based Hilbert space engine. - """ - - def __init__(self, V:str, compl : bool = True, dim : int = 2, - meshname : str = "Th", Vname : str = "V"): - if compl: - self.compl = "" - else: - self.compl = "" - self.V = V - self.dim = dim - self.meshname = meshname - self.Vname = Vname - - def name(self) -> str: - """Class label.""" - return self.__class__.__name__ - - def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> float: - """ - Compute general norm of complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - normType(optional): Target norm identifier. If matrix, target norm - is one induced by normType. If number, target norm is weighted - H^1 norm with given weight. If string, must be among 'H1', - 'L2', and 'H10'. Defaults to 'H10'. - - Returns: - Norm of the function (non-negative). - """ - if type(normType).__name__[-6:] == "matrix": - return np.abs(u.dot(normType.dot(u).conj())) ** .5 - filename = utilities.getNewFilename("temp", "edp") - ufile = FFCT.npvec2ffvec(u) - if normType == 'L2': - coeff1 = 0. - coeff2 = 1. - else: - coeff1 = 1. - if isinstance(normType, (int, float)): - coeff2 = normType - elif normType == 'H10': - coeff2 = 0. - elif normType == 'H1': - coeff2 = 1. - else: - raise Exception("Norm type not recognized.") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("ifstream fin(\"{}\");\n".format(ufile)) - f.write("{}{} u;\n".format(self.Vname, self.compl)) - f.write("fin >> u[];\n") - f.write("cout.precision(15); ") - f.write("cout.scientific << sqrt(") - if np.isclose(coeff1, 0.): - f.write("0") - else: - f.write("int{}d({})({})".format(self.dim, self.meshname, - FFCT.diff("u", self.dim))) - if np.isclose(coeff2, 0.): - f.write(");\n") - else: - f.write(" + {} * int{}d({})(abs(u)^2));\n".format( - coeff2, - self.dim, - self.meshname)) - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE) - output, error = process.communicate() - try: - nrm = np.float(output) - except: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee files {} and {}.")\ - .format(output.decode('utf-8'), filename, ufile)) - os.remove(filename) - os.remove(ufile) - return nrm - - def plot(self, u:Np1D, name : str = "u", save : bool = False, - what : strLst = "ALL", figspecs : str = "fill = 1, value = 1"): - """ - Do some nice plots of the complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - name(optional): Name to be shown as title of the plots. Defaults to - 'u'. - save(optional): Whether to save plot(s). Defaults to False. - figspecs(optional): Optional arguments for FreeFem++ figure - creation. Defaults to "fill = 1, value = 1". - """ - if isinstance(what, (str,)): - what = what.upper() - if what != 'ALL': - what = [what] - elif save: - warnings.warn(("Saving and 'ALL' are not compatible. " - "Overriding 'ALL' to ['ABS', 'PHASE', 'REAL', " - "'IMAG']."), stacklevel = 2) - what = ['ABS', 'PHASE', 'REAL', 'IMAG'] - if what != 'ALL': - what = utilities.purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], - listname = self.name() + ".what") - if len(what) == 0: return - - filename = utilities.getNewFilename("temp", "edp") - ufile = FFCT.npvec2ffvec(u) - - figspecs = figspecs.strip() - if len(figspecs) > 0 and figspecs[0] != ",": - figspecs = ", " + figspecs - - - with open(filename, 'w') as f: - f.write(FFCT.HSVjet + "\n") - f.write(self.V + "\n") - f.write("ifstream fin(\"{}\");\n".format(ufile)) - f.write("{}{} u;\n".format(self.Vname, self.compl)) - f.write("fin >> u[];\n") - if isinstance(what, (str,)): - f.write(("plot(u, wait = 1, cmm = \"{}\"{}, " - "hsv = HSVjet);\n").format(name, figspecs)) - else: - for i, w in enumerate(what): - if i == 0: - VnameEff = self.Vname - else: - VnameEff = "" - if w == 'ABS': - outspecs = "" - if save: - outspecs = (", ps = \"" - + utilities.getNewFilename("fig", "abs.eps") - + "\"") - f.write("{} v = abs(u);\n".format(VnameEff)) - f.write(("plot(v, wait = 1, cmm = \"|{}|\"{}{}, " - "hsv = HSVjet);\n").format(name, figspecs, - outspecs)) - if w == 'PHASE': - outspecs = "" - if save: - outspecs = (", ps = \"" - + utilities.getNewFilename("fig", "ph.eps") - + "\"") - f.write("{} v = arg(u);\n".format(VnameEff)) - f.write(("plot(v, wait = 1, cmm = \"angle({})\"{}{}, " - "hsv = HSVjet);\n").format(name, figspecs, - outspecs)) - if w == 'REAL': - outspecs = "" - if save: - outspecs = (", ps = \"" - + utilities.getNewFilename("fig", "re.eps") - + "\"") - f.write("{} v = real(u);\n".format(VnameEff)) - f.write(("plot(v, wait = 1, cmm = \"Re({})\"{}{}, " - "hsv = HSVjet);\n").format(name, figspecs, - outspecs)) - if w == 'IMAG': - outspecs = "" - if save: - outspecs = (", ps = \"" - + utilities.getNewFilename("fig", "im.eps") - + "\"") - f.write("{} v = imag(u);\n".format(VnameEff)) - f.write(("plot(v, wait = 1, cmm = \"Im({})\"{}{}, " - "hsv = HSVjet);\n").format(name, figspecs, - outspecs)) - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee files {} and {}.")\ - .format(output.decode('utf-8'), filename, ufile)) - os.remove(filename) - os.remove(ufile) - - def plotmesh(self, name : str = "Mesh", save : bool = False, - figspecs : str = ""): - """ - Do a nice plot of the mesh. - - Args: - name(optional): Name to be shown as title of the plots. Defaults to - 'Mesh'. - save(optional): Whether to save plot(s). Defaults to False. - figspecs(optional): Optional arguments for FreeFem++ figure - creation. - """ - filename = utilities.getNewFilename("temp", "edp") - - figspecs = figspecs.strip() - if len(figspecs) > 0 and figspecs[0] != ",": - figspecs = ", " + figspecs - outspecs = "" - if save: - outspecs = (", ps = \"" + utilities.getNewFilename("msh", "eps") - + "\"") - - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("plot({}, wait = 1, cmm=\"{}\"{}{});\n".format( - self.meshname, name, figspecs, outspecs)) - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - os.remove(filename) - - def convExpression(self, expr:str) -> Np1D: - """ - Evaluate dofs of general expression. - - Args: - expr: expression. - - Returns: - Numpy 1D array with dofs. - """ - filename = utilities.getNewFilename("temp", "edp") - efile = utilities.getNewFilename("e", "tmp") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} expr = {};\n".format(self.Vname, self.compl, - expr.replace("j", "i"))) - f.write("ofstream eout(\"{}\");\n".format(efile)) - f.write("eout << expr[];\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee files {} and {}.")\ - .format(output.decode('utf-8'), filename, efile)) - out = FFCT.ffvec2npvec(efile) - os.remove(filename) - os.remove(efile) - return out - - -class FreeFemHSAugmentedEngine(FreeFemHSEngine): - """ - FreeFem-based Hilbert space engine for augmented spaces. - """ - def __init__(self, V:str, d : int = 2, compl : bool = True, dim : int = 2, - meshname : str = "Th", Vname : str = "V"): - self.d = d - FreeFemHSEngine.__init__(self, V = V, compl = compl, dim = dim, - meshname = meshname, Vname = Vname) - - def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> Np1D: - """ - Compute general norm of complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - normType(optional): Target norm identifier. If matrix, target norm - is one induced by normType. If number, target norm is weighted - H^1 norm with given weight. If string, must be among 'H1', - 'L2', and 'H10'. Defaults to 'H10'. - - Returns: - Norms of the function (non-negative). - """ - if isinstance(normType, (int, float, str)): - u = utilities.listify(u, self.d) - norms = [None] * self.d - for j in range(self.d): - norms[j] = FreeFemHSEngine.norm(self, u[j], normType) - return norms - elif type(normType).__name__[-6:] == "matrix": - u = utilities.vectorify(u, self.d) - if normType.shape[0] % u[0] == 0: - N = int(normType.shape[0] / u[0]) - else: - raise Exception("Energy matrix dimension mismatch.") - return FreeFemHSEngine.norm(self, np.concatenate(u[: N]), normType) - raise Exception("Norm type not recognized.") - - def plot(self, u:Np1D, split : bool = True, name : str = "u", - save : bool = False, what : strLst = "ALL", - figspecs : str = "fill = 1, value = 1"): - """ - Do some nice plots of the complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - split: whether to split the given array. - name(optional): Name to be shown as title of the plots. Defaults to - 'u'. - save(optional): Whether to save plot(s). Defaults to False. - what(optional): Which plots to do. If list, can contain 'ABS', - 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. - Defaults to 'ALL'. - figspecs(optional): Optional arguments for FreeFem++ figure - creation. Defaults to "fill = 1, value = 1". - """ - if split: - u = utilities.listify(u, self.d) - - for j in range(self.d): - FreeFemHSEngine.plot(self, u[j], what = what, save = save, - name = "{0}, comp. {1}".format(name, j), - figspecs = figspecs) - else: - FreeFemHSEngine.plot(self, u, what = what, save = save, - name = name, figspecs = figspecs) - diff --git a/main/FreeFemHelmholtzEngine.py b/main/FreeFemHelmholtzEngine.py deleted file mode 100644 index 03b1223..0000000 --- a/main/FreeFemHelmholtzEngine.py +++ /dev/null @@ -1,739 +0,0 @@ -#!/usr/bin/python - -import warnings -import os -import subprocess -import numpy as np -import scipy.sparse as scsp -import utilities -import FreeFemConversionTools as FFCT -from RROMPyTypes import Np1D, Np2D, DictAny, Tuple, List - -class DomainError(Exception): - """ - Domain error exception. - - Args: - value: Human readable string describing the exception. - - Attributes: - value: Human readable string describing the exception. - """ - def __init__(self, value:str): - self.value = value - def __str__(self): - return self.value - -class FreeFemHelmholtzEngine: - """ - FreeFem++-based solver for Helmholtz problems. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu + h u = g2 on \Gamma_R - - Args: - - Attributes: - """ - - def __init__(self, V:str, - wavenumber:complex, - compl : bool = True, - dim : int = 2, - meshname : str = "Th", - Vname : str = "V", - diffusivity : str = "1", - refractionIndex : str = "1", - forcingTerm : str = "0", - DirichletBoundary : str = "", - NeumannBoundary : str = "", - RobinBoundary : str = "", - DirichletDatum : str = "0", - NeumannDatum : str = "0", - RobinDatum : str = "0", - RobinDatum2 : str = "0", - tgv : float = 1e30): - # process boundaries - boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary] - for i in range(3): - boundariesList[i] = boundariesList[i].strip() - DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList - if boundariesList.count("") == 3: - raise DomainError("At least one boundary must be prescribed.") - self.DirichletBoundary = DirichletBoundary - self.NeumannBoundary = NeumannBoundary - self.RobinBoundary = RobinBoundary - - if compl: - self.compl = "" - else: - self.compl = "" - self.V = V - self.dim = dim - self.meshname = meshname - self.Vname = Vname - self.wavenumber2 = wavenumber ** 2. - self.diffusivity = diffusivity - self.refractionIndex = refractionIndex - self.forcingTerm = forcingTerm - self.DirichletDatum = DirichletDatum - self.NeumannDatum = NeumannDatum - self.RobinDatum = RobinDatum - self.RobinDatum2 = RobinDatum2 - self.tgv = tgv - - functional = lambda self, u: 0. - - def Vdim(self) -> int: - """Dimension of finite element space.""" - if not hasattr(self, "liftDData"): self.liftDirichletData() - return len(self.liftDData) - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - filename = utilities.getNewFilename("temp", "edp") - Mfile = utilities.getNewFilename("m", "tmp") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} u, v;\n".format(self.Vname, self.compl)) - f.write("varf energy(u, v) = int{}d({})({})".format( - self.dim, self.meshname, - FFCT.diff2(dim = self.dim, - compl = not self.compl == ""))) - f.write(" + int{}d({})({} * {});\n".format( - self.dim, self.meshname, w, - FFCT.prod2(compl = not self.compl == ""))) - f.write("matrix{0} M = energy({1}, {1});\n".format(self.compl, - self.Vname)) - f.write("ofstream Mout(\"{}\");\n".format(Mfile)) - f.write("Mout << M;\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - M = FFCT.ffmat2spmat(Mfile) - os.remove(filename) - os.remove(Mfile) - return M.tocsr() - - def problemData(self) -> DictAny: - """List of HF problem data.""" - dataDict = {} - dataDict["forcingTerm"] = self.forcingTerm - dataDict["DirichletDatum"] = self.DirichletDatum - dataDict["NeumannDatum"] = self.NeumannDatum - dataDict["RobinDatum2"] = self.RobinDatum2 - return dataDict - - def setProblemData(self, dataDict:DictAny, k2:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k2: Parameter value. - """ - self.wavenumber2 = k2 - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "b0"): del self.b0 - - def liftDirichletData(self): - """Lift Dirichlet datum.""" - if not hasattr(self, "liftDData"): - filename = utilities.getNewFilename("temp", "edp") - u0file = utilities.getNewFilename("u0", "tmp") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} u0 = {};\n".format(self.Vname, self.compl, - self.DirichletDatum.replace('j', 'i'))) - f.write("ofstream u0out(\"{}\");\n".format(u0file)) - f.write("u0out << u0[];\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), - stdout=subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - self.liftDData = FFCT.ffvec2npvec(u0file) - os.remove(filename) - os.remove(u0file) - return self.liftDData - - @property - def diffusivity(self): - """Value of a.""" - return self._diffusivity - @diffusivity.setter - def diffusivity(self, diffusivity): - if hasattr(self, "A0"): del self.A0 - self._diffusivity = diffusivity - - @property - def refractionIndex(self): - """Value of n.""" - return self._refractionIndex - @refractionIndex.setter - def refractionIndex(self, refractionIndex): - if hasattr(self, "A1"): del self.A1 - self._refractionIndex = refractionIndex - - @property - def forcingTerm(self): - """Value of f.""" - return self._forcingTerm - @forcingTerm.setter - def forcingTerm(self, forcingTerm): - if hasattr(self, "b0"): del self.b0 - self._forcingTerm = forcingTerm - - @property - def DirichletDatum(self): - """Value of u0.""" - return self._DirichletDatum - @DirichletDatum.setter - def DirichletDatum(self, DirichletDatum): - if hasattr(self, "b0"): del self.b0 - if hasattr(self, "liftDData"): del self.liftDData - self._DirichletDatum = DirichletDatum - - @property - def NeumannDatum(self): - """Value of g1.""" - return self._NeumannDatum - @NeumannDatum.setter - def NeumannDatum(self, NeumannDatum): - if hasattr(self, "b0"): del self.b0 - self._NeumannDatum = NeumannDatum - - @property - def RobinDatum(self): - """Value of h.""" - return self._RobinDatum - @RobinDatum.setter - def RobinDatum(self, RobinDatum): - if hasattr(self, "A0"): del self.A0 - self._RobinDatum = RobinDatum - - @property - def RobinDatum2(self): - """Value of g2.""" - return self._RobinDatum2 - @RobinDatum2.setter - def RobinDatum2(self, RobinDatum2): - if hasattr(self, "b0"): del self.b0 - self._RobinDatum2 = RobinDatum2 - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - if not (hasattr(self, "A0") and hasattr(self, "A1")): - filename = utilities.getNewFilename("temp", "edp") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} u, v;\n".format(self.Vname, self.compl)) - if not hasattr(self, "A0"): - A0file = utilities.getNewFilename("A0", "tmp") - f.write("varf a0(u, v) = int{}d({})({} * {})".format( - self.dim, self.meshname, - self.diffusivity.replace('j', 'i'), - FFCT.diff2(dim = self.dim, - compl = not self.compl == ""))) - if len(self.RobinBoundary) > 0: - f.write(" + int{}d({}, {})({} * {})".format( - self.dim - 1, self.meshname, - self.RobinBoundary, - self.RobinDatum.replace('j', 'i'), - FFCT.prod2(compl = not self.compl == ""))) - if len(self.DirichletBoundary) > 0: - f.write(" + on({}, u = 0.);\n".format( - self.DirichletBoundary)) - else: - f.write(";\n") - f.write("matrix{0} A0 = a0({1}, {1}, tgv = {2});\n".format( - self.compl, self.Vname, self.tgv)) - f.write("ofstream A0out(\"{}\");\n".format(A0file)) - f.write("A0out << A0;\n") - if not hasattr(self, "A1"): - A1file = utilities.getNewFilename("A1", "tmp") - f.write("varf a1(u, v) = -int{}d({})(({})^2*{});\n".format( - self.dim, self.meshname, - self.refractionIndex.replace('j', 'i'), - FFCT.prod2(compl = not self.compl == ""))) - f.write("matrix{0} A1 = a1({1},{1});\n".format( - self.compl, self.Vname)) - f.write("ofstream A1out(\"{}\");\n".format(A1file)) - f.write("A1out << A1;\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), - stdout = subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - if not hasattr(self, "A0"): - self.A0, self.DBCRows = FFCT.ffmat2spmatnormalize(A0file, - self.tgv) - os.remove(A0file) - if not hasattr(self, "A1"): - self.A1 = FFCT.ffmat2spmatcleanrows(A1file, self.DBCRows) - os.remove(A1file) - os.remove(filename) - - def assembleb(self): - """Assemble RHS blocks of linear system.""" - if not hasattr(self, "b0"): - filename = utilities.getNewFilename("temp", "edp") - b0file = utilities.getNewFilename("b0", "tmp") - if self.compl == "": - complEff = "real" - else: - complEff = "complex" - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} u, v;\n".format(self.Vname, self.compl)) - f.write("varf b0(u, v) = int{}d({})({})".format( - self.dim, self.meshname, - FFCT.prod2(self.forcingTerm.replace('j', 'i'), - compl = not self.compl == ""))) - if len(self.NeumannBoundary) > 0: - f.write(" + int{}d({}, {})({})".format( - self.dim - 1, self.meshname, - self.NeumannBoundary, - FFCT.prod2(self.NeumannDatum.replace('j', 'i'), - compl = not self.compl == ""))) - if len(self.RobinBoundary) > 0: - f.write(" + int{}d({}, {})({})".format( - self.dim - 1, self.meshname, - self.RobinBoundary, - FFCT.prod2(self.RobinDatum2.replace('j', 'i'), - compl = not self.compl == ""))) - if len(self.DirichletBoundary) > 0: - f.write(" + on({}, u = {});\n".format( - self.DirichletBoundary, - self.DirichletDatum.replace('j', 'i'))) - else: - f.write(";\n") - f.write("{}[int] B0 = b0(0, {}, tgv = {});\n".format(complEff, - self.Vname, self.tgv)) - f.write("ofstream b0out(\"{}\");\n".format(b0file)) - f.write("b0out << B0;\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), - stdout = subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - if len(self.DirichletBoundary) > 0 and not hasattr(self,'DBCRows'): - warnings.warn(("Dirichlet boundary condition management for " - "RHS requires system matrix. Proceeding with " - "matrix assembly."), stacklevel = 2) - self.assembleA() - self.b0 = FFCT.ffvec2npveccleanrows(b0file, self.tgv, self.DBCRows) - os.remove(b0file) - os.remove(filename) - - parDegree = 1 - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1] - - def bs(self) -> List[Np1D]: - """RHS blocks of linear system.""" - self.assembleb() - return [self.b0] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber2 * self.A1 - - def b(self) -> Np1D: - """Assemble RHS of linear system.""" - self.assembleb() - return self.b0 - - def solve(self) -> Np1D: - """ - Find solution of linear system. - - Returns: - Solution vector. - """ - return scsp.linalg.spsolve(self.A(), self.b()) - - def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]: - """ - Get blocks of linear system. - - Returns: - Blocks of LS (\sum_{j=0}^J k^j A_j = \sum_{j=0}^K k^j b_j). - """ - return self.As(), self.bs() - - -class FreeFemHelmholtzScatteringEngine(FreeFemHelmholtzEngine): - """ - FreeFem++-based solver for Helmholtz scattering problems with Robin ABC. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu - i k u = g2 on \Gamma_R - - Args: - - Attributes: - """ - - def __init__(self, V:str, - wavenumber:complex, - compl : bool = True, - dim : int = 2, - meshname : str = "Th", - Vname : str = "V", - diffusivity : str = "1", - refractionIndex : str = "1", - forcingTerm : str = "0", - DirichletBoundary : str = "", - NeumannBoundary : str = "", - RobinBoundary : str = "", - DirichletDatum : str = "0", - NeumannDatum : str = "0", - RobinDatum2 : str = "0", - tgv : float = 1e30): - self.wavenumber = wavenumber - FreeFemHelmholtzEngine.__init__(self, V = V, - wavenumber = wavenumber, - compl = compl, - dim = dim, - meshname = meshname, - Vname = Vname, - diffusivity = diffusivity, - refractionIndex = refractionIndex, - forcingTerm = forcingTerm, - DirichletBoundary = DirichletBoundary, - NeumannBoundary = NeumannBoundary, - RobinBoundary = RobinBoundary, - DirichletDatum = DirichletDatum, - NeumannDatum = NeumannDatum, - RobinDatum = -1.j * wavenumber, - RobinDatum2 = RobinDatum2, - tgv = tgv) - if len(self.RobinBoundary) == 0: - raise Exception(("Empty Robin boundary not allowed in scattering " - "engine. Use standard engine instead.")) - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - return FreeFemHelmholtzEngine.energyNormMatrix(self, w ** 2.) - - def setProblemData(self, dataDict:DictAny, k:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k: Parameter value. - """ - self.wavenumber = k - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "A2"): del self.A2 - if hasattr(self, "b0"): del self.b0 - - @property - def wavenumber2(self): - """Value of k^2.""" - return self._wavenumber2 - @wavenumber2.setter - def wavenumber2(self, wavenumber2): - self._wavenumber2 = wavenumber2 - self._wavenumber = self.wavenumber2 ** .5 - self.RobinDatum = "{}".format(-1.j*self.wavenumber) - - @property - def wavenumber(self): - """Value of k.""" - return self._wavenumber - @wavenumber.setter - def wavenumber(self, wavenumber): - self.wavenumber2 = wavenumber ** 2. - - @property - def RobinDatum(self): - """Value of h.""" - return super().RobinDatum - @RobinDatum.setter - def RobinDatum(self, RobinDatum): - try: - self._RobinDatum = "{}".format(complex(RobinDatum)) - except: - raise Exception(("Value of RobinDatum not recognized. Must be " - "complex number of str convertible to complex.")) - if not (hasattr(self, 'wavenumber') - and np.isclose(self.wavenumber, 1.j * complex(self.RobinDatum))): - self.wavenumber = 1.j * complex(self.RobinDatum) - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - if not (hasattr(self, "A0") and hasattr(self, "A1") - and hasattr(self, "A2")): - filename = utilities.getNewFilename("temp", "edp") - with open(filename, 'w') as f: - f.write(self.V + "\n") - f.write("{}{} u, v;\n".format(self.Vname, self.compl)) - if not hasattr(self, "A0"): - A0file = utilities.getNewFilename("A0", "tmp") - f.write("varf a0(u, v) = int{}d({})({} * {})".format( - self.dim, self.meshname, - self.diffusivity.replace('j', 'i'), - FFCT.diff2(dim = self.dim, - compl = not self.compl == ""))) - if len(self.DirichletBoundary) > 0: - f.write(" + on({}, u = 0.);".format( - self.DirichletBoundary)) - f.write(";\n") - f.write("matrix{0} A0 = a0({1}, {1}, tgv = {2});\n".format( - self.compl, self.Vname, self.tgv)) - f.write("ofstream A0out(\"{}\");\n".format(A0file)) - f.write("A0out << A0;\n") - if not hasattr(self, "A1"): - A1file = utilities.getNewFilename("A1", "tmp") - f.write("varf a1(u, v) = -int{}d({},{})(1.i*{});\n".format( - self.dim - 1, self.meshname, - self.RobinBoundary, - FFCT.prod2(compl = not self.compl == ""))) - f.write("matrix{0} A1 = a1({1},{1});\n".format(self.compl, - self.Vname)) - f.write("ofstream A1out(\"{}\");\n".format(A1file)) - f.write("A1out << A1;\n") - if not hasattr(self, "A2"): - A2file = utilities.getNewFilename("A2", "tmp") - f.write("varf a2(u, v) = - int{}d({})({} * {});\n".format( - self.dim, self.meshname, - self.refractionIndex.replace('j', 'i'), - FFCT.prod2(compl = not self.compl == ""))) - f.write("matrix{0} A2 = a2({1},{1});\n".format(self.compl, - self.Vname)) - f.write("ofstream A2out(\"{}\");\n".format(A2file)) - f.write("A2out << A2;\n") - bashCommand = "FreeFem++ -v 0 {}".format(filename) - process = subprocess.Popen(bashCommand.split(), - stdout = subprocess.PIPE) - output, error = process.communicate() - if len(output) > 0: - raise Exception(("Error was encountered in the execution of " - "FreeFem++ code:\n{}\nSee file {}.")\ - .format(output.decode('utf-8'), filename)) - if not hasattr(self, "A0"): - self.A0, self.DBCRows = FFCT.ffmat2spmatnormalize(A0file, - self.tgv) - os.remove(A0file) - if not hasattr(self, "A1"): - self.A1 = FFCT.ffmat2spmatcleanrows(A1file, self.DBCRows) - os.remove(A1file) - if not hasattr(self, "A2"): - self.A2 = FFCT.ffmat2spmatcleanrows(A2file, self.DBCRows) - os.remove(A2file) - os.remove(filename) - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1, self.A2] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber * self.A1 + self.wavenumber2 * self.A2 - - -class FreeFemHelmholtzScatteringAugmentedEngine( - FreeFemHelmholtzScatteringEngine): - """ - FreeFem++-based solver for Helmholtz scattering problems with Robin ABC. - - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega - u = u0 on \Gamma_D - \partial_nu = g1 on \Gamma_N - \partial_nu - i k u = g2 on \Gamma_R - Linear dependence on k is achieved by introducing an additional variable. - - Args: - - Attributes: - """ - - def __init__(self, V:str, - wavenumber:complex, - compl : bool = True, - dim : int = 2, - meshname : str = "Th", - Vname : str = "V", - diffusivity : str = "1", - refractionIndex : str = "1", - forcingTerm : str = "0", - DirichletBoundary : str = "", - NeumannBoundary : str = "", - RobinBoundary : str = "", - DirichletDatum : str = "0", - NeumannDatum : str = "0", - RobinDatum2 : str = "0", - tgv : float = 1e30, - constraintType : str = "IDENTITY"): - self.constraintType = constraintType - FreeFemHelmholtzScatteringEngine.__init__(self, V = V, - wavenumber = wavenumber, - compl = compl, - dim = dim, - meshname = meshname, - Vname = Vname, - diffusivity = diffusivity, - refractionIndex = refractionIndex, - forcingTerm = forcingTerm, - DirichletBoundary = DirichletBoundary, - NeumannBoundary = NeumannBoundary, - RobinBoundary = RobinBoundary, - DirichletDatum = DirichletDatum, - NeumannDatum = NeumannDatum, - RobinDatum2 = RobinDatum2, - tgv = tgv) - - def energyNormMatrix(self, w:float) -> Np2D: - """ - Sparse matrix (in CSR format) assoociated to w-weighted H10 norm. - - Args: - w: Weight. - - Returns: - Sparse matrix in CSR format. - """ - normMatFen = FreeFemHelmholtzScatteringEngine.energyNormMatrix(self, w) - return scsp.block_diag((normMatFen, 1/w * normMatFen)) - - def liftDirichletData(self): - """Lift Dirichlet datum.""" - lift1 = FreeFemHelmholtzScatteringEngine.liftDirichletData(self) - return np.pad(lift1, (0, len(lift1)), 'constant') - - def setProblemData(self, dataDict:DictAny, k:complex): - """ - Set problem data. - - Args: - dataDict: Dictionary of problem data. - k: Parameter value. - """ - self.wavenumber = k - self.forcingTerm = dataDict["forcingTerm"] - self.DirichletDatum = dataDict["DirichletDatum"] - self.NeumannDatum = dataDict["NeumannDatum"] - self.RobinDatum2 = dataDict["RobinDatum2"] - if hasattr(self, "A0"): del self.A0 - if hasattr(self, "A1"): del self.A1 - if hasattr(self, "b0"): del self.b0 - - @property - def constraintType(self): - """Value of constraintType.""" - return self._constraintType - @constraintType.setter - def constraintType(self, constraintType): - try: - constraintType = constraintType.upper().strip().replace(" ","") - if constraintType not in ["IDENTITY", "MASS"]: raise - self._constraintType = constraintType - except: - warnings.warn(("Prescribed constraintType not recognized. " - "Overriding to 'IDENTITY'."), stacklevel = 2) - self.constraintType = "IDENTITY" - - def assembleA(self): - """Assemble matrix blocks of linear system.""" - builtA0 = hasattr(self, "A0") - builtA1 = hasattr(self, "A1") - if not (builtA0 and builtA1): - if builtA1 and self.constraintType == "IDENTITY": - self.A2 = None - FreeFemHelmholtzScatteringEngine.assembleA(self) - if self.constraintType == "IDENTITY": - if not builtA0: - block = scsp.eye(self.A0.shape[0]) - else: - block = scsp.eye(self.A1.shape[0]) - else: #MASS - block = - self.A2 - I = list(self.DirichletBC0.get_boundary_values().keys()) - warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning) - block[I, I] = self.A0[I[0], I[0]] - warnings.simplefilter("default", scsp.SparseEfficiencyWarning) - if not builtA0: - self.A0 = scsp.block_diag((self.A0, block)) - if not builtA1: - self.A1 = scsp.bmat([[self.A1, self.A2], [- block, None]]) - if hasattr(self, 'A2'): del self.A2 - - def assembleb(self): - """Assemble RHS of linear system.""" - if not hasattr(self, "b0"): - FreeFemHelmholtzScatteringEngine.assembleb(self) - self.b0 = np.pad(self.b0, (0, self.b0.shape[0]), 'constant') - - parDegree = 1 - - def As(self) -> List[Np2D]: - """Matrix blocks of linear system.""" - self.assembleA() - return [self.A0, self.A1] - - def A(self) -> Np2D: - """Assemble matrix of linear system.""" - self.assembleA() - return self.A0 + self.wavenumber * self.A1 - - def solve(self) -> Tuple[Np1D, Np1D]: - """ - Find solution of linear system. - - Returns: - Solution vector. - """ - u = FreeFemHelmholtzScatteringEngine.solve(self) - lu2 = int(len(u) / 2) - return np.array([u[: lu2], u[lu2 :]]) - - diff --git a/main/ROMApproximantLagrangeVF.py b/main/ROMApproximantLagrangeVF.py deleted file mode 100644 index ea9d627..0000000 --- a/main/ROMApproximantLagrangeVF.py +++ /dev/null @@ -1,312 +0,0 @@ -#!/usr/bin/python - -from copy import copy -import warnings -import numpy as np -import utilities -from ROMApproximantLagrange import ROMApproximantLagrange -from RROMPyTypes import Np1D, ListAny - -class ROMApproximantLagrangeVF(ROMApproximantLagrange): - """ - ROM Lagrange VF interpolant computation for parametric problems. - - Args: - HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - - energyNormMatrix: assemble sparse matrix (in CSC format) - associated to weighted H10 norm; - - problemData: list of HF problem data (excluding mu); - - setProblemData: set HF problem data and mu to prescribed values; - - getLSBlocks: get algebraic system blocks. - HSEngine: Hilbert space general purpose engine. Should have members: - - norm: compute Hilbert norm of expression represented as complex - numpy vector; - - plot: plot Hilbert expression represented as complex numpy - vector. - mus(optional): Boundaries of snapshot positions. Defaults to - np.array([0, 1]). - ws(optional): Array of snapshot weigths. Defaults to - np.ones_like(self.mus). - w(optional): Weight for norm computation. If None, set to - Re(np.mean(self.mus)). Defaults to None. - approxParameters(optional): Dictionary containing values for main - parameters of approximant. Recognized keys are: - - 'POD': whether to compute POD of snapshots; defaults to True; - - 'S': total number of samples current approximant relies upon; - defaults to 2; - - 'M': degree of VF interpolant numerator; defaults to 0; - - 'N': degree of VF interpolant denominator; defaults to 0. - Defaults to empty dict. - plotSnap(optional): What to plot of snapshots of the Helmholtz - solution map. See plot method in HSEngine. Defaults to []. - plotSnapSpecs(optional): How to plot snapshots of the Helmholtz - solution map. See plot method in HSEngine. Defaults to {}. - - Attributes: - HFEngine: HF problem solver. - HSEngine: Hilbert space general purpose engine. - mu0: Default parameter. - mus: Boundaries of snapshot positions. - w: Weight for norm computation. - approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: - - 'POD': whether to compute POD of snapshots; - - 'S': total number of samples current approximant relies upon; - - 'M': degree of VF interpolant numerator; - - 'N': degree of VF interpolant denominator; - extraApproxParameters: List of approxParameters keys in addition to - mother class's. - autoNodeTypes: List of possible values for autoNode. - M: Numerator degree of approximant. - N: Denominator degree of approximant. - S: Number of solution snapshots over which current approximant is - based upon. - plotSnap: What to plot of snapshots of the Helmholtz solution map. - plotSnapSpecs: How to plot snapshots of the Helmholtz solution map. - solSnapshots: Array whose columns are FE dofs of snapshots of solution - map. - RPOD: Right factor of snapshots POD. If no POD, set to None. - POD: Whether to compute POD of snapshots. - autoNode: Type of nodes, if automatically generated. Otherwise None. - energyNormMatrix: Sparse matrix (in CSC format) assoociated to - w-weighted H10 norm. - uHF: High fidelity solution with wavenumber lastSolvedHF as numpy - complex vector. - lastSolvedHF: Wavenumber corresponding to last computed high fidelity - solution. - Q: Numpy 1D vector containing complex coefficients of approximant - denominator. - P: Numpy 2D vector whose columns are FE dofs of coefficients of - approximant numerator. - uApp: Last evaluated approximant as numpy complex vector. - lastApproxParameters: List of parameters corresponding to last - computed approximant. - As: List of sparse matrices (in CSC format) representing coefficients - of linear system matrix wrt mu. - bs: List of numpy vectors representing coefficients of linear system - RHS wrt mu. - """ - - extraApproxParameters = ["M", "N"] - - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximantLagrange.parameterList(self) - + ROMApproximantLagrangeVF.extraApproxParameters) - - @property - def approxParameters(self): - """ - Value of approximant parameters. Its assignment may change M, N and S. - """ - return self._approxParameters - @approxParameters.setter - def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantLagrangeVF.parameterList(self), - dictname = self.name() + ".approxParameters") - keyList = list(approxParameters.keys()) - if "M" in keyList: - self.M = 0 #to avoid warnings if M and S are changed simultaneously - if "N" in keyList: - self.N = 0 #to avoid warnings if N and S are changed simultaneously - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantLagrangeVF.extraApproxParameters, - True, True) - ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy) - if "M" in keyList: - self.M = approxParameters["M"] - elif hasattr(self, "M"): - self.M = self.M - else: - self.M = 0 - if "N" in keyList: - self.N = approxParameters["N"] - elif hasattr(self, "N"): - self.N = self.N - else: - self.N = 0 - - @property - def M(self): - """Value of M. Its assignment may change S.""" - return self._M - @M.setter - def M(self, M): - if M < 0: raise ArithmeticError("M must be non-negative.") - self._M = M - self._approxParameters["M"] = self.M - if hasattr(self, "S") and self.S < self.M + 1: - warnings.warn("Prescribed S is too small. Updating S to M + 1.", - stacklevel = 2) - self.S = self.M + 1 - - @property - def N(self): - """Value of N. Its assignment may change S.""" - return self._N - @N.setter - def N(self, N): - if N < 0: raise ArithmeticError("N must be non-negative.") - self._N = N - self._approxParameters["N"] = self.N - if hasattr(self, "S") and self.S < self.N + 1: - warnings.warn("Prescribed S is too small. Updating S to N + 1.", - stacklevel = 2) - self.S = self.N + 1 - - @property - def S(self): - """Value of S.""" - return self._S - @S.setter - def S(self, S): - if S <= 0: raise ArithmeticError("S must be positive.") - if hasattr(self, "S"): Sold = self.S - else: Sold = -1 - vals, label = [0] * 2, {0:"M", 1:"N"} - if hasattr(self, "M"): vals[0] = self.M - if hasattr(self, "N"): vals[1] = self.N - idxmax = np.argmax(vals) - if vals[idxmax] + 1 > S: - warnings.warn("Prescribed S is too small. Updating S to {} + 1."\ - .format(label[idxmax]), stacklevel = 2) - self.Emax = vals[idxmax] + 1 - else: - self._S = S - self._approxParameters["S"] = self.S - if Sold != self.S: - self.resetSamples() - - def approxRadius(self) -> float: - """Sample radius.""" - return np.max(self.mu0 - self.mus) - - def setupApprox(self): - """Compute VF interpolant. SVD-based robust eigenvalue management.""" - if not self.checkComputedApprox(): - S0 = self.S - M1 = self.M + 1 - N1 = self.N + 1 - if self.autoNode not in self.autoNodeTypes[1 :]: - raise Exception("Manual sample selection not yet implemented.") - self.computeSnapshots() - Phi = np.zeros((S0, max(M1, N1)), dtype = np.complex) - Phi[:, 0] = np.ones((S0,)) / 2**.5 - for j in range(1, max(M1, N1)): - c = np.zeros((j + 1,)) - c[-1] = 1. - if self.polyBasis == "CHEBYSHEV": - Phi[:, j] = np.polynomial.chebyshev.chebval( - self.radiusVF(self.mus), c) - elif self.polyBasis == "GAUSSLEGENDRE": - Phi[:, j] = (j + .5) ** .5 * np.polynomial.legendre.legval( - self.radiusVF(self.mus), c) - - if self.POD: - II = np.array(np.arange(0, S0**3, S0) - + np.kron(np.arange(S0), np.ones(S0)), - dtype = np.int) - Rtall = np.zeros(S0**3, dtype = np.complex) - Rtall[II] = np.reshape(self.RPOD.T, (S0**2,)) - Rtall = np.reshape(Rtall, (S0, S0**2)).T - - B = Rtall.dot(Phi[:, : N1]) - Z = copy(B) - Z = np.reshape(Z.T, (S0 * (N1), S0)).T - - for j in range(2): #twice is enough! - Z = Z - Phi[:, : M1].dot(Phi[:, : M1].conj().T.dot( - np.multiply(self.ws, Z))) - Z = np.reshape(Z.T, (N1, S0**2)).T - else: - ker = self.solSnapshots.conj().T.dot(self.energyNormMat.dot( - self.solSnapshots)) - WPhi = np.reshape(np.multiply(self.ws, Phi[:, : M1]).T, - (S0 * M1)).conj()[:, None] - Y = np.multiply(WPhi, np.kron(np.ones((M1, 1)), Phi[:, : N1])) - Ylarge = np.reshape(Y.T, (M1 * N1, S0)).T - - B = np.reshape(ker.dot(Ylarge).T, (N1, M1 * S0)).T - D = np.multiply(self.ws, np.diag(ker)[:, None]) - D = Phi[:, : N1].conj().T.dot(np.multiply(D, Phi[:, : N1])) - - Z = B.conj().T.dot(Y) - D - _, ev, eV = np.linalg.svd(Z, full_matrices = False) - phi = eV[-1, :].conj() - - if self.POD: - c = Phi[:, : M1].conj().T.dot(np.multiply(self.ws, - np.reshape(B.dot(phi).T, (S0, S0)).T)) - else: - c = np.reshape(Y.dot(phi).T, (M1, S0)) - - polybase = np.zeros((max(M1, N1), max(M1, N1))) - polybase[0, 0] = 1 - polybase[1, 1] = 1 - for j in range(2, max(M1, N1)): - if self.polyBasis == "CHEBYSHEV": - polybase[1 : j + 1, j] = 2 * polybase[: j, j - 1] - polybase[: j + 1, j] = (polybase[: j + 1, j] - - polybase[: j + 1, j - 2]) - elif self.polyBasis == "GAUSSLEGENDRE": - polybase[1 : j + 1, j] = ((2 * j - 1) / j - * polybase[: j, j - 1]) - polybase[: j + 1, j] = (polybase[: j + 1, j] - - ((j - 1) / j * polybase[: j + 1, j - 2])) - if self.polyBasis == "CHEBYSHEV": - polybase[0, 0] = .5 ** .5 - elif self.polyBasis == "GAUSSLEGENDRE": - polybase = np.multiply(np.power( - np.arange(.5, max(M1, N1)), .5), polybase) - self.P = self.solSnapshots.dot(polybase[: M1, : M1].dot(c).T) - self.Q = polybase[: N1, : N1].dot(phi) - - self.approxParameters = {"N" : self.N, "M" : self.M, "S" : S0} - self.lastApproxParameters = copy(self.approxParameters) - - def radiusVF(self, mu:complex) -> float: - """ - Compute translated radius to be plugged into VF approximant. - - Args: - mu: Target parameter. - - Returns: - Translated radius to be plugged into VF approximant. - """ - return (mu - self.mu0) / self.approxRadius() - - def evalApprox(self, mu:complex) -> Np1D: - """ - Evaluate VF approximant at arbitrary wavenumber. - - Args: - mu: Target parameter. - - Returns: - Real and imaginary parts of approximant. - """ - self.setupApprox() - powerlist = np.power(self.radiusVF(mu), range(max(self.M, self.N) + 1)) - self.uApp = (self.P.dot(powerlist[:self.M + 1]) - / self.Q.dot(powerlist[:self.N + 1])) - return self.uApp - - def getPoles(self) -> Np1D: - """ - Obtain approximant poles. - - Returns: - Numpy complex vector of poles. - """ - self.setupApprox() - return np.roots(self.Q[::-1]) * self.approxRadius() + self.mu0 diff --git a/main/RROMPyTypes.py b/main/RROMPyTypes.py deleted file mode 100644 index 3e0c93b..0000000 --- a/main/RROMPyTypes.py +++ /dev/null @@ -1,21 +0,0 @@ -#!/usr/bin/python - -from typing import TypeVar, List, Tuple, Dict, Any - -TupleAny = Tuple[Any] -ListAny = List[Any] -DictAny = Dict[Any, Any] -strLst = TypeVar("str or list of str") -Np1D = TypeVar("NumPy 1D array") -Np2D = TypeVar("NumPy 2D array") -Np1DLst = TypeVar("NumPy 1D array or list of NumPy 1D array") -GenExpr = TypeVar("Generic expression") -FenSpace = TypeVar("FEniCS FESpace") -FenExpr = TypeVar("FEniCS expression") -FenFunc = TypeVar("FEniCS function") -BfSExpr = TypeVar("Boolean function or string") -N2FSExpr = TypeVar("NumPy 2D array, float or str") - -HFEng = TypeVar("High fidelity engine") -HSEng = TypeVar("Hilbert space engine") -ROMEng = TypeVar("ROM engine") diff --git a/main/utilities.py b/main/utilities.py deleted file mode 100644 index f4bdbe0..0000000 --- a/main/utilities.py +++ /dev/null @@ -1,162 +0,0 @@ -#!/usr/bin/python - -import os -import warnings -import numpy as np -from RROMPyTypes import Any, ListAny, DictAny, Np1DLst, List, Np1D - -def findDictStrKey(key:Any, keyList:ListAny): - for akey in keyList: - if isinstance(key, str) and key.lower() == akey.lower(): - return akey - return None - -def purgeList(lst:ListAny, allowedEntries : ListAny = [], - silent : bool = False, complement : bool = False, - listname : str = ""): - if listname != "": - listname = " in " + listname - lstcp = [] - for x in lst: - ax = findDictStrKey(x, allowedEntries) - if (ax is None) != complement: - if not silent: - warnings.warn("Ignoring entry {0}{1}.".format(x, listname), - stacklevel = 2) - else: - lstcp = lstcp + [ax] - return lstcp - -def purgeDict(dct:DictAny, allowedKeys : ListAny = [], silent : bool = False, - complement : bool = False, dictname : str = ""): - if dictname != "": - dictname = " in " + dictname - dctcp = {} - for key in dct.keys(): - akey = findDictStrKey(key, allowedKeys) - if (akey is None) != complement: - if not silent: - warnings.warn("Ignoring key {0}{2} with value {1}."\ - .format(key, dct[key], dictname), stacklevel = 2) - else: - if akey is None: - akey = key - dctcp[akey] = dct[key] - return dctcp - -def getNewFilename(prefix : str = "", extension : str = "dat") -> str: - n = len(extension) + 1 - filename = "{}{}.{}".format(prefix, np.random.randint(0, 10), extension) - while os.path.exists(filename): - filename = filename[: - n] + "{}.{}".format(np.random.randint(10), - extension) - return filename - -prime_v = [] #memoization vector - -def squareResonances(a:int, b:int, zero_terms : bool = True): - spectrum = [] - a = max(int(np.floor(a)), 0) - b = max(int(np.ceil(b)), 0) - global prime_v - if len(prime_v) == 0: - prime_v = [2, 3] - if a > prime_v[-1]: - for i in range(prime_v[-1], a, 2): - get_next_prime_factor(i) - for i in range(a, b + 1): - spectrum = spectrum + [i] * count_square_sums(i, zero_terms) - return np.array(spectrum) - -def get_next_prime_factor(n:int): - global prime_v - for x in prime_v: - if n % x == 0: - return x - if prime_v[-1] < n: - prime_v = prime_v + [n] - return n - -def prime_factorize(n:int): - factors = [] - number = n - while number > 1: - factor = get_next_prime_factor(number) - factors.append(factor) - number = int(number / factor) - if n < -1: - factors[0] = -factors[0] - return list(factors) - -def count_square_sums(n:int, zero_terms : bool = True): - if n < 2: return (n + 1) * zero_terms - factors = prime_factorize(n) - funique, fcounts = np.unique(factors, return_counts = True) - count = 1 - for fac, con in zip(funique, fcounts): #using number theory magic - if fac % 4 == 1: - count = count * (con + 1) - elif fac % 4 == 3 and con % 2 == 1: - return 0 - return count + (2 * zero_terms - 1) * (int(n ** .5) ** 2 == n) - - -def listify(u:Np1DLst, d:int) -> List[Np1D]: - """ - Convert numpy array into list-like format. - - Args: - u: numpy complex array with function dofs or list of such. - d: secondary dimension. - - Returns: - List-like numpy array of numpy arrays. - """ - if isinstance(u, (list,)) or len(u) == d: - if len(u) != d: - raise Exception(("Input list lenght must be equal to " - "augmentation dimension d.")) - N = len(u[0]) - for x in u[1 :]: - if len(x) != N: - raise Exception(("Elements of input list must all have " - "the same lenght.")) - return np.array(u) - else: - N = int(len(u) / d) - if not np.isclose(len(u), d * N): - raise Exception(("Input vector lenght must be multiple of " - "augmentation dimension d.")) - v = [None] * d - for i in range(d): - v[i] = u[i * N : (i + 1) * N] - return np.array(v) - -def vectorify(u:Np1DLst, d:int) -> Np1D: - """ - Convert list-like object into numpy array. - - Args: - u: numpy complex array with function dofs or list of such. - d: secondary dimension. - - Returns: - Numpy arrays. - """ - if isinstance(u, (list,)) or len(u) == d: - if len(u) != d: - raise Exception(("Input list lenght must be equal to " - "augmentation dimension d.")) - N = len(u[0]) - for x in u[1 :]: - if len(x) != N: - raise Exception(("Elements of input list must all have " - "the same lenght.")) - return np.concatenate(u) - else: - N = int(len(u) / d) - if not np.isclose(len(u), d * N): - raise Exception(("Input vector lenght must be multiple of" - "augmentation dimension d.")) - return u - diff --git a/rrompy/__init__.py b/rrompy/__init__.py new file mode 100644 index 0000000..892c0e6 --- /dev/null +++ b/rrompy/__init__.py @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 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 . +# + diff --git a/rrompy/hfengines/__init__.py b/rrompy/hfengines/__init__.py new file mode 100644 index 0000000..892c0e6 --- /dev/null +++ b/rrompy/hfengines/__init__.py @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 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 . +# + diff --git a/rrompy/hfengines/fenics/__init__.py b/rrompy/hfengines/fenics/__init__.py new file mode 100644 index 0000000..52463b6 --- /dev/null +++ b/rrompy/hfengines/fenics/__init__.py @@ -0,0 +1,39 @@ +#!/usr/bin/python +# 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.hfengines.fenics.generic_problem_engine import GenericProblemEngine +from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine +from rrompy.hfengines.fenics.helmholtz_box_scattering_problem_engine import HelmholtzBoxScatteringProblemEngine +from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine +from rrompy.hfengines.fenics.helmholtz_square_bubble_problem_engine import HelmholtzSquareBubbleProblemEngine +from rrompy.hfengines.fenics.helmholtz_square_transmission_problem_engine import HelmholtzSquareTransmissionProblemEngine +from rrompy.hfengines.fenics.scattering_problem_engine import ScatteringProblemEngine + +__all__ = [ + 'GenericProblemEngine', + 'HelmholtzBaseProblemEngine', + 'HelmholtzBoxScatteringProblemEngine', + 'HelmholtzProblemEngine', + 'HelmholtzSquareBubbleProblemEngine', + 'HelmholtzSquareTransmissionProblemEngine', + 'ScatteringProblemEngine' + ] + + + diff --git a/rrompy/hfengines/fenics/generic_problem_engine.py b/rrompy/hfengines/fenics/generic_problem_engine.py new file mode 100644 index 0000000..154dec1 --- /dev/null +++ b/rrompy/hfengines/fenics/generic_problem_engine.py @@ -0,0 +1,81 @@ +#!/usr/bin/python +# 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 abc import abstractmethod +import numpy as np +import scipy.sparse as scsp +from rrompy.utilities.types import Np1D, Np2D, Tuple, List + +__all__ = ['GenericProblemEngine'] + +class GenericProblemEngine: + """ + Generic solver for parametric problems. + + Attributes: + As: Scipy sparse array representation (in CSC format) of As. + bs: Numpy array representation of bs. + """ + Aterms, bterms = 1, 1 + + functional = lambda self, u: 0. + + @abstractmethod + def energyNormMatrix(self) -> Np2D: + """Sparse matrix (in CSR format) representative of scalar product.""" + pass + + @abstractmethod + def assembleA(self): + """Assemble matrix blocks of linear system.""" + self.As = [None] * 1 + pass + + @abstractmethod + def assembleb(self): + """Assemble matrix blocks of linear system.""" + self.bs = [None] * 1 + pass + + def A(self, mu:complex) -> Np2D: + """Assemble matrix of linear system.""" + self.assembleA() + A = self.As[0] + for j in range(1, len(self.As)): + A = A + np.power(mu, j) * self.As[j] + return A + + def b(self, mu:complex) -> Np1D: + """Assemble RHS of linear system.""" + self.assembleb() + b = self.bs[0] + for j in range(1, len(self.bs)): + b = b + np.power(mu, j) * self.bs[j] + return b + + def solve(self, mu:complex) -> Np1D: + """Find solution of linear system.""" + return scsp.linalg.spsolve(self.A(mu), self.b(mu)) + + def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]: + """Get blocks of linear system.""" + self.assembleA() + self.assembleb() + return self.As, self.bs + diff --git a/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py new file mode 100644 index 0000000..5e9f67f --- /dev/null +++ b/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py @@ -0,0 +1,308 @@ +#!/usr/bin/python +# 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 +import scipy.sparse as scsp +import fenics as fen +from rrompy.hfengines.fenics.generic_problem_engine import GenericProblemEngine +from rrompy.utilities.types import Np1D, Np2D +from rrompy.utilities.fenics import fenZERO, fenONE, bdrTrue, bdrFalse + +__all__ = ['HelmholtzBaseProblemEngine'] + +class HelmholtzBaseProblemEngine(GenericProblemEngine): + """ + ABSTRACT + Solver for generic Helmholtz problems with parametric wavenumber. + - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega + u = u0 on \Gamma_D + \partial_nu = g1 on \Gamma_N + \partial_nu + h u = g2 on \Gamma_R + + Attributes: + omega: Value of omega. + diffusivity: Value of a. + refractionIndex: Value of n. + forcingTerm: Value of f. + DirichletDatum: Value of u0. + NeumannDatum: Value of g1. + RobinDatumG: Value of g2. + DirichletBoundary: Function handle to \Gamma_D. + NeumannBoundary: Function handle to \Gamma_N. + RobinBoundary: Function handle to \Gamma_R. + V: Real FE space. + u: Generic trial functions for variational form evaluation. + v: Generic test functions for variational form evaluation. + ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). + A0: Scipy sparse array representation (in CSC format) of A0. + A1: Scipy sparse array representation (in CSC format) of A1. + b0: Numpy array representation of b0. + As: Scipy sparse array representation (in CSC format) of As. + bs: Numpy array representation of bs. + """ + Aterms = 2 + omega = 1. + _diffusivity = [fenONE, fenZERO] + _refractionIndex = [fenONE, fenZERO] + _forcingTerm = [fenZERO, fenZERO] + _DirichletDatum = [fenZERO, fenZERO] + _NeumannDatum = [fenZERO, fenZERO] + _RobinDatumG = [fenZERO, fenZERO] + + def __init__(self): + self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1) + self.DirichletBoundary = "ALL" + + @property + def V(self): + """Value of V.""" + return self._V + @V.setter + def V(self, V): + if hasattr(self, "A0"): del self.A0 + if hasattr(self, "A1"): del self.A1 + if hasattr(self, "b0"): del self.b0 + if not type(V).__name__ == 'FunctionSpace': + raise Exception("V type not recognized.") + self._V = V + self.u = fen.TrialFunction(V) + self.v = fen.TestFunction(V) + self.dsToBeSet = True + + @property + def diffusivity(self): + """Value of a.""" + return self._diffusivity + @diffusivity.setter + def diffusivity(self, diffusivity): + if hasattr(self, "A0"): del self.A0 + if not isinstance(diffusivity, (list, tuple,)): + diffusivity = [diffusivity, fenZERO] + self._diffusivity = diffusivity + + @property + def refractionIndex(self): + """Value of n.""" + return self._refractionIndex + @refractionIndex.setter + def refractionIndex(self, refractionIndex): + if hasattr(self, "A1"): del self.A1 + if not isinstance(refractionIndex, (list, tuple,)): + refractionIndex = [refractionIndex, fenZERO] + self._refractionIndex = refractionIndex + + @property + def forcingTerm(self): + """Value of f.""" + return self._forcingTerm + @forcingTerm.setter + def forcingTerm(self, forcingTerm): + if hasattr(self, "b0"): del self.b0 + if not isinstance(forcingTerm, (list, tuple,)): + forcingTerm = [forcingTerm, fenZERO] + self._forcingTerm = forcingTerm + + @property + def DirichletDatum(self): + """ + Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and + DirichletBCIm. + """ + return self._DirichletDatum + @DirichletDatum.setter + def DirichletDatum(self, DirichletDatum): + if hasattr(self, "b0"): del self.b0 + if not isinstance(DirichletDatum, (list, tuple,)): + DirichletDatum = [DirichletDatum, fenZERO] + self._DirichletDatum = DirichletDatum + + @property + def NeumannDatum(self): + """Value of g1.""" + return self._NeumannDatum + @NeumannDatum.setter + def NeumannDatum(self, NeumannDatum): + if hasattr(self, "b0"): del self.b0 + if not isinstance(NeumannDatum, (list, tuple,)): + NeumannDatum = [NeumannDatum, fenZERO] + self._NeumannDatum = NeumannDatum + + @property + def RobinDatumG(self): + """Value of g2.""" + return self._RobinDatumG + @RobinDatumG.setter + def RobinDatumG(self, RobinDatumG): + if hasattr(self, "b0"): del self.b0 + if not isinstance(RobinDatumG, (list, tuple,)): + RobinDatumG = [RobinDatumG, fenZERO] + self._RobinDatumG = RobinDatumG + + @property + def DirichletBoundary(self): + """Function handle to DirichletBoundary.""" + return self._DirichletBoundary + @DirichletBoundary.setter + def DirichletBoundary(self, DirichletBoundary): + if hasattr(self, "A0"): del self.A0 + if hasattr(self, "A1"): del self.A1 + if isinstance(DirichletBoundary, (str,)): + if DirichletBoundary.upper() == "NONE": + self._DirichletBoundary = bdrFalse + self.DREST = 0 + elif DirichletBoundary.upper() == "ALL": + self._DirichletBoundary = bdrTrue + self.NeumannBoundary = "NONE" + self.RobinBoundary = "NONE" + self.DREST = 0 + elif DirichletBoundary.upper() == "REST": + if self.NREST + self.RREST > 0: + raise Exception("Only 1 'REST' wildcard can be specified.") + self._DirichletBoundary = lambda x, on_boundary : (on_boundary + and not self.NeumannBoundary(x, on_boundary) + and not self.RobinBoundary(x, on_boundary)) + self.DREST = 1 + else: + raise Exception("DirichletBoundary label not recognized.") + elif callable(DirichletBoundary): + self._DirichletBoundary = DirichletBoundary + self.DREST = 0 + else: + raise Exception("DirichletBoundary type not recognized.") + + @property + def NeumannBoundary(self): + """Function handle to NeumannBoundary.""" + return self._NeumannBoundary + @NeumannBoundary.setter + def NeumannBoundary(self, NeumannBoundary): + if hasattr(self, "A1"): del self.A1 + self.dsToBeSet = True + if isinstance(NeumannBoundary, (str,)): + if NeumannBoundary.upper() == "NONE": + self._NeumannBoundary = bdrFalse + self.NREST = 0 + elif NeumannBoundary.upper() == "ALL": + self._NeumannBoundary = bdrTrue + self.DirichletBoundary = "NONE" + self.RobinBoundary = "NONE" + self.NREST = 0 + elif NeumannBoundary.upper() == "REST": + if self.DREST + self.RREST > 0: + raise Exception("Only 1 'REST' wildcard can be specified.") + self._NeumannBoundary = lambda x, on_boundary : (on_boundary + and not self.DirichletBoundary(x, on_boundary) + and not self.RobinBoundary(x, on_boundary)) + self.NREST = 1 + else: + raise Exception("NeumannBoundary label not recognized.") + elif callable(NeumannBoundary): + self._NeumannBoundary = NeumannBoundary + self.NREST = 0 + else: + raise Exception("DirichletBoundary type not recognized.") + + @property + def RobinBoundary(self): + """Function handle to RobinBoundary.""" + return self._RobinBoundary + @RobinBoundary.setter + def RobinBoundary(self, RobinBoundary): + if hasattr(self, "A0"): del self.A0 + if hasattr(self, "A1"): del self.A1 + self.dsToBeSet = True + if isinstance(RobinBoundary, (str,)): + if RobinBoundary.upper() == "NONE": + self._RobinBoundary = bdrFalse + self.RREST = 0 + elif RobinBoundary.upper() == "ALL": + self._RobinBoundary = bdrTrue + self.DirichletBoundary = "NONE" + self.NeumannBoundary = "NONE" + self.RREST = 0 + elif RobinBoundary.upper() == "REST": + if self.DREST + self.NREST > 0: + raise Exception("Only 1 'REST' wildcard can be specified.") + self._RobinBoundary = lambda x, on_boundary : (on_boundary + and not self.DirichletBoundary(x, on_boundary) + and not self.NeumannBoundary(x, on_boundary)) + self.RREST = 1 + else: + raise Exception("RobinBoundary label not recognized.") + return + elif callable(RobinBoundary): + self._RobinBoundary = RobinBoundary + self.RREST = 0 + else: + raise Exception("RobinBoundary type not recognized.") + + def autoSetDS(self): + """Set FEniCS boundary measure based on boundary function handles.""" + if self.dsToBeSet: + NB = self.NeumannBoundary + RB = self.RobinBoundary + class NBoundary(fen.SubDomain): + def inside(self, x, on_boundary): + return NB(x, on_boundary) + class RBoundary(fen.SubDomain): + def inside(self, x, on_boundary): + return RB(x, on_boundary) + boundary_markers = fen.MeshFunction("size_t", self.V.mesh(), + self.V.mesh().topology().dim() - 1) + NBoundary().mark(boundary_markers, 0) + RBoundary().mark(boundary_markers, 1) + self.ds = fen.Measure("ds", domain = self.V.mesh(), + subdomain_data = boundary_markers) + self.dsToBeSet = False + + def energyNormMatrix(self) -> Np2D: + """Sparse matrix (in CSR format) representative of scalar product.""" + normMatFen = fen.assemble((fen.dot(fen.grad(self.u), fen.grad(self.v)) + + np.abs(self.omega)**2 * fen.dot(self.u, self.v)) * fen.dx) + normMat = fen.as_backend_type(normMatFen).mat() + return scsp.csr_matrix(normMat.getValuesCSR()[::-1], + shape = normMat.size) + + def liftDirichletData(self) -> Np1D: + """Lift Dirichlet datum.""" + solLRe = fen.interpolate(self.DirichletDatum[0], self.V) + solLIm = fen.interpolate(self.DirichletDatum[1], self.V) + return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector()) + + def assembleb(self): + """Assemble matrix blocks of linear system.""" + if not hasattr(self, "b0"): + fRe, fIm = self.forcingTerm + g1Re, g1Im = self.NeumannDatum + g2Re, g2Im = self.RobinDatumG + L0Re = (fen.dot(fRe, self.v) * fen.dx + + fen.dot(g1Re, self.v) * self.ds(0) + + fen.dot(g2Re, self.v) * self.ds(1)) + L0Im = (fen.dot(fIm, self.v) * fen.dx + + fen.dot(g1Im, self.v) * self.ds(0) + + fen.dot(g2Im, self.v) * self.ds(1)) + b0Re = fen.assemble(L0Re) + b0Im = fen.assemble(L0Im) + fen.DirichletBC(self.V, self.DirichletDatum[0], + self.DirichletBoundary).apply(b0Re) + fen.DirichletBC(self.V, self.DirichletDatum[1], + self.DirichletBoundary).apply(b0Im) + self.b0 = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) + self.bs = [self.b0] + diff --git a/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py new file mode 100644 index 0000000..a5dbde0 --- /dev/null +++ b/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py @@ -0,0 +1,58 @@ +#!/usr/bin/python +# 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 +import fenics as fen +from rrompy.hfengines.fenics.scattering_problem_engine import ScatteringProblemEngine + +__all__ = ['HelmholtzBoxScatteringProblemEngine'] + +class HelmholtzBoxScatteringProblemEngine(ScatteringProblemEngine): + """ + Solver for scattering problem outside a box with parametric wavenumber. + - \Delta u - omega^2 * n^2 * u = 0 in \Omega + u = 0 on \Gamma_D + \partial_nu - i k u = 0 on \Gamma_R + with exact solution a transmitted plane wave. + """ + def __init__(self, R:float, kappa:float, theta:float, n:int): + super().__init__(self) + + import mshr + scatterer = mshr.Polygon([fen.Point(-1, -.5), fen.Point(1, -.5), + fen.Point(1, .5), fen.Point(.8, .5), + fen.Point(.8, -.3), fen.Point(-.8, -.3), + fen.Point(-.8, .5), fen.Point(-1, .5),]) + mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R)-scatterer, n) + self.V = fen.FunctionSpace(mesh, "P", 3) + + self.DirichletBoundary = (lambda x, on_boundary: + on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R) + self.RobinBoundary = (lambda x, on_boundary: + on_boundary and (x[0]**2+x[1]**2)**.5 > .95 * R) + + import sympy as sp + x, y = sp.symbols('x[0] x[1]', real=True) + phiex = kappa * (x * np.cos(theta) + y * np.sin(theta)) + u0ex = - sp.exp(1.j * phiex) + u0 = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), + sp.printing.ccode(sp.simplify(sp.im(u0ex)))] + self.DirichletDatum = [fen.Expression(x, degree = 3) for x in u0] + + diff --git a/rrompy/hfengines/fenics/helmholtz_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_problem_engine.py new file mode 100644 index 0000000..27e44a4 --- /dev/null +++ b/rrompy/hfengines/fenics/helmholtz_problem_engine.py @@ -0,0 +1,119 @@ +#!/usr/bin/python +# 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 scipy.sparse as scsp +import fenics as fen +from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine, fenZERO + +__all__ = ['HelmholtzProblemEngine'] + +class HelmholtzProblemEngine(HelmholtzBaseProblemEngine): + """ + Solver for Helmholtz problems with parametric wavenumber. + - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega + u = u0 on \Gamma_D + \partial_nu = g1 on \Gamma_N + \partial_nu + h u = g2 on \Gamma_R + + Attributes: + omega: Value of omega. + diffusivity: Value of a. + refractionIndex: Value of n. + forcingTerm: Value of f. + DirichletDatum: Value of u0. + NeumannDatum: Value of g1. + RobinDatumH: Value of h. + RobinDatumG: Value of g2. + DirichletBoundary: Function handle to \Gamma_D. + NeumannBoundary: Function handle to \Gamma_N. + RobinBoundary: Function handle to \Gamma_R. + V: Real FE space. + u: Generic trial functions for variational form evaluation. + v: Generic test functions for variational form evaluation. + ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). + A0: Scipy sparse array representation (in CSC format) of A0. + A1: Scipy sparse array representation (in CSC format) of A1. + b0: Numpy array representation of b0. + As: Scipy sparse array representation (in CSC format) of As. + bs: Numpy array representation of bs. + """ + _RobinDatumH = [fenZERO, fenZERO] + + def __init__(self): + super().__init__(self) + self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1) + self.DirichletBoundary = "ALL" + + @property + def RobinDatumH(self): + """Value of h.""" + return self._RobinDatumH + @RobinDatumH.setter + def RobinDatumH(self, RobinDatumH): + if hasattr(self, "A0"): del self.A0 + if not isinstance(RobinDatumH, (list, tuple,)): + RobinDatumH = [RobinDatumH, fenZERO] + self._RobinDatumH = RobinDatumH + + def assembleA(self): + """Assemble matrix blocks of linear system.""" + self.autoSetDS() + if not hasattr(self, "A0"): + DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + aRe, aIm = self.diffusivity + hRe, hIm = self.RobinDatumH + a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx + + hRe * fen.dot(self.u, self.v) * self.ds(1)) + a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx + + hIm * fen.dot(self.u, self.v) * self.ds(1)) + A0Re = fen.assemble(a0Re) + A0Im = fen.assemble(a0Im) + DirichletBC0.apply(A0Re) + DirichletBC0.zero(A0Im) + A0ReMat = fen.as_backend_type(A0Re).mat() + A0ImMat = fen.as_backend_type(A0Im).mat() + A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR() + A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR() + self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), + shape = A0ReMat.size) + + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), + shape = A0ImMat.size)) + if not hasattr(self, "A1"): + DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + nRe, nIm = self.refractionIndex + n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm + a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx + a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx + A1Re = fen.assemble(a1Re) + A1Im = fen.assemble(a1Im) + DirichletBC0.zero(A1Re) + DirichletBC0.zero(A1Im) + A1ReMat = fen.as_backend_type(A1Re).mat() + A1ImMat = fen.as_backend_type(A1Im).mat() + A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR() + A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR() + self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer), + shape = A1ReMat.size) + + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr), + shape = A1ImMat.size)) + self.As = [self.A0, self.A1] + + diff --git a/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py new file mode 100644 index 0000000..f0565ab --- /dev/null +++ b/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py @@ -0,0 +1,48 @@ +#!/usr/bin/python +# 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 +import fenics as fen +from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine + +__all__ = ['HelmholtzSquareBubbleProblemEngine'] + +class HelmholtzSquareBubbleProblemEngine(HelmholtzProblemEngine): + """ + Solver for square bubble Helmholtz problems with parametric wavenumber. + - \Delta u - omega^2 * u = f in \Omega + u = 0 on \Gamma_D + with exact solution square bubble times plane wave. + """ + def __init__(self, kappa:float, theta:float, n:int): + super().__init__(self) + + mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n) + self.V = fen.FunctionSpace(mesh, "P", 3) + + import sympy as sp + x, y = sp.symbols('x[0] x[1]', real=True) + wex = 16/np.pi**4 * x * y * (x - np.pi) * (y - np.pi) + phiex = kappa * (x * np.cos(theta) + y * np.sin(theta)) + uex = wex * sp.exp(-1.j * phiex) + fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex + forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), + sp.printing.ccode(sp.simplify(sp.im(fex)))] + self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm] + diff --git a/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py new file mode 100644 index 0000000..687a092 --- /dev/null +++ b/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py @@ -0,0 +1,61 @@ +#!/usr/bin/python +# 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 +import fenics as fen +from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine + +__all__ = ['HelmholtzSquareTransmissionProblemEngine'] + +class HelmholtzSquareTransmissionProblemEngine(HelmholtzProblemEngine): + """ + Solver for square transmission Helmholtz problems with parametric + wavenumber. + - \Delta u - omega^2 * n^2 * u = 0 in \Omega + u = 0 on \Gamma_D + with exact solution a transmitted plane wave. + """ + def __init__(self, nT:float, nB:float, kappa:float, theta:float, n:int): + super().__init__(self) + + mesh = fen.RectangleMesh(fen.Point(-np.pi/2, -np.pi/2), + fen.Point(np.pi/2, np.pi/2), n, n) + self.V = fen.FunctionSpace(mesh, "P", 3) + + import sympy as sp + dx, dy = np.cos(theta), np.sin(theta) + Kx = kappa * nB * dx + Ky = kappa * (nT**2. - (nB * dx)**2. + 0.j)**.5 + T = 2 * kappa * nB * dy / (Ky + kappa * nB * dy) + x, y = sp.symbols('x[0] x[1]', real=True) + uT = T * sp.exp(1.j * (Kx*x + Ky*y)) + uB = ( sp.exp(1.j * kappa * nB * (dx*x + dy*y)) + + (T - 1) * sp.exp(1.j * kappa * nB * (dx*x - dy*y))) + uRe = fen.Expression('x[1]>=0 ? {} : {}'.format( + sp.printing.ccode(sp.re(uT)), + sp.printing.ccode(sp.re(uB))), + degree = 3) + uIm = fen.Expression('x[1]>=0 ? {} : {}'.format( + sp.printing.ccode(sp.im(uT)), + sp.printing.ccode(sp.im(uB))), + degree = 3) + + self.refractionIndex = fen.Expression('x[1] >= 0 ? nT : nB', + nT = nT, nB = nB, degree = 0) + self.DirichletDatum = [uRe, uIm] diff --git a/rrompy/hfengines/fenics/scattering_problem_engine.py b/rrompy/hfengines/fenics/scattering_problem_engine.py new file mode 100644 index 0000000..c53fc90 --- /dev/null +++ b/rrompy/hfengines/fenics/scattering_problem_engine.py @@ -0,0 +1,114 @@ +#!/usr/bin/python +# 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 scipy.sparse as scsp +import fenics as fen +from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine, fenZERO + +__all__ = ['ScatteringProblemEngine'] + +class ScatteringProblemEngine(HelmholtzBaseProblemEngine): + """ + Solver for scattering problems with parametric wavenumber. + - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega + u = u0 on \Gamma_D + \partial_nu = g1 on \Gamma_N + \partial_nu +- i k u = g2 on \Gamma_R + + Attributes: + signR: Sign in ABC. + omega: Value of omega. + diffusivity: Value of a. + refractionIndex: Value of n. + forcingTerm: Value of f. + DirichletDatum: Value of u0. + NeumannDatum: Value of g1. + RobinDatumG: Value of g2. + DirichletBoundary: Function handle to \Gamma_D. + NeumannBoundary: Function handle to \Gamma_N. + RobinBoundary: Function handle to \Gamma_R. + V: Real FE space. + u: Generic trial functions for variational form evaluation. + v: Generic test functions for variational form evaluation. + ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). + A0: Scipy sparse array representation (in CSC format) of A0. + A1: Scipy sparse array representation (in CSC format) of A1. + A2: Scipy sparse array representation (in CSC format) of A1. + b0: Numpy array representation of b0. + As: Scipy sparse array representation (in CSC format) of As. + bs: Numpy array representation of bs. + """ + Aterms = 3 + signR = - 1. + + def __init__(self): + self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1) + self.DirichletBoundary = "ALL" + + def assembleA(self): + """Assemble matrix blocks of linear system.""" + self.autoSetDS() + if not hasattr(self, "A0"): + DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + aRe, aIm = self.diffusivity + a0Re = aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx + a0Im = aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx + A0Re = fen.assemble(a0Re) + A0Im = fen.assemble(a0Im) + DirichletBC0.apply(A0Re) + DirichletBC0.zero(A0Im) + A0ReMat = fen.as_backend_type(A0Re).mat() + A0ImMat = fen.as_backend_type(A0Im).mat() + A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR() + A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR() + self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), + shape = A0ReMat.size) + + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), + shape = A0ImMat.size)) + if not hasattr(self, "A1"): + DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + a1 = fen.dot(self.u, self.v) * self.ds(1) + A1 = fen.assemble(a1) + DirichletBC0.zero(A1) + A1Mat = fen.as_backend_type(A1).mat() + A1r, A1c, A1v = A1Mat.getValuesCSR() + self.A1 = self.signR * 1.j * scsp.csr_matrix((A1v, A1c, A1r), + shape = A1Mat.size) + if not hasattr(self, "A2"): + DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + nRe, nIm = self.refractionIndex + n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm + a2Re = - n2Re * fen.dot(self.u, self.v) * fen.dx + a2Im = - n2Im * fen.dot(self.u, self.v) * fen.dx + A2Re = fen.assemble(a2Re) + A2Im = fen.assemble(a2Im) + DirichletBC0.zero(A2Re) + DirichletBC0.zero(A2Im) + A2ReMat = fen.as_backend_type(A2Re).mat() + A2ImMat = fen.as_backend_type(A2Im).mat() + A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR() + A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR() + self.A2 = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer), + shape = A2ReMat.size) + + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr), + shape = A2ImMat.size)) + self.As = [self.A0, self.A1, self.A2] diff --git a/rrompy/hsengines/__init__.py b/rrompy/hsengines/__init__.py new file mode 100644 index 0000000..892c0e6 --- /dev/null +++ b/rrompy/hsengines/__init__.py @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 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 . +# + diff --git a/rrompy/hsengines/fenics/__init__.py b/rrompy/hsengines/fenics/__init__.py new file mode 100644 index 0000000..45f8c5f --- /dev/null +++ b/rrompy/hsengines/fenics/__init__.py @@ -0,0 +1,26 @@ +#!/usr/bin/python +# 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.hsengines.fenics.hsengine import HSEngine + +__all__ = [ + 'HSEngine' + ] + + diff --git a/main/FenicsHSEngine.py b/rrompy/hsengines/fenics/hsengine.py similarity index 57% rename from main/FenicsHSEngine.py rename to rrompy/hsengines/fenics/hsengine.py index 386678b..a1b0751 100644 --- a/main/FenicsHSEngine.py +++ b/rrompy/hsengines/fenics/hsengine.py @@ -1,198 +1,152 @@ #!/usr/bin/python +# 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 -import utilities from matplotlib import pyplot as plt import fenics as fen -from RROMPyTypes import Np1D, strLst, FenSpace, N2FSExpr, List +from rrompy.utilities.types import Np1D, strLst, FenSpace, N2FSExpr +from rrompy.utilities import purgeList, getNewFilename -class FenicsHSEngine: +__all__ = ['HSEngine'] + +class HSEngine: """ Fenics-based Hilbert space engine. """ def __init__(self, V:FenSpace): self.V = V def name(self) -> str: """Class label.""" return self.__class__.__name__ def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> float: """ Compute general norm of complex-valued function with given dofs. Args: u: numpy complex array with function dofs. normType(optional): Target norm identifier. If matrix, target norm is one induced by normType. If number, target norm is weighted H^1 norm with given weight. If string, must be recognizable by Fenics norm command. Defaults to 'H10'. Returns: Norm of the function (non-negative). """ if type(normType).__name__[-6:] == "matrix": return np.abs(u.dot(normType.dot(u).conj())) ** .5 if isinstance(normType, (int, float)): - return (FenicsHSEngine.norm(self, u, "H10")**2 - + (normType * FenicsHSEngine.norm(self, u, "L2"))**2)**.5 + return ( HSEngine.norm(self, u, "H10")**2 + + (normType * HSEngine.norm(self, u, "L2"))**2)**.5 uAbs = fen.Function(self.V) uAbs.vector()[:] = np.array(np.abs(u), dtype = float) return fen.norm(uAbs, normType) def plot(self, u:Np1D, name : str = "u", save : bool = False, what : strLst = 'all', **figspecs): """ Do some nice plots of the complex-valued function with given dofs. Args: u: numpy complex array with function dofs. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Whether to save plot(s). Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ if isinstance(what, (str,)): if what.upper() == 'ALL': what = ['ABS', 'PHASE', 'REAL', 'IMAG'] else: what = [what] - what = utilities.purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], - listname = self.name() + ".what") + what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], + listname = self.name() + ".what") if len(what) == 0: return if 'figsize' not in figspecs.keys(): figspecs['figsize'] = (13. * len(what) / 4, 3) subplotcode = 100 + len(what) * 10 if 'REAL' in what and 'IMAG' in what: combined_data_ReIm = np.concatenate([np.real(u), np.imag(u)]) _min = np.amin(combined_data_ReIm) _max = np.amax(combined_data_ReIm) pars = {'vmin' : _min, 'vmax' : _max} else: pars = {} plt.figure(**figspecs) plt.jet() if 'ABS' in what: uAb = fen.Function(self.V) uAb.vector()[:] = np.array(np.abs(u), dtype = float) subplotcode = subplotcode + 1 plt.subplot(subplotcode) p = fen.plot(uAb, title = "|{0}|".format(name)) plt.colorbar(p) if 'PHASE' in what: uPh = fen.Function(self.V) uPh.vector()[:] = np.array(np.angle(u), dtype = float) subplotcode = subplotcode + 1 plt.subplot(subplotcode) p = fen.plot(uPh, title = "phase({0})".format(name)) plt.colorbar(p) if 'REAL' in what: uRe = fen.Function(self.V) uRe.vector()[:] = np.array(np.real(u), dtype = float) subplotcode = subplotcode + 1 plt.subplot(subplotcode) p = fen.plot(uRe, title = "Re({0})".format(name), **pars) plt.colorbar(p) if 'IMAG' in what: uIm = fen.Function(self.V) uIm.vector()[:] = np.array(np.imag(u), dtype = float) subplotcode = subplotcode + 1 plt.subplot(subplotcode) p = fen.plot(uIm, title = "Im({0})".format(name), **pars) plt.colorbar(p) if save: - plt.savefig(utilities.getNewFilename("fig", "eps"), - format='eps', dpi=1000) + plt.savefig(getNewFilename("fig", "eps"), format='eps', dpi=1000) plt.show() plt.close() def plotmesh(self, name : str = "Mesh", save : bool = False, **figspecs): """ Do a nice plot of the mesh. Args: u: numpy complex array with function dofs. name(optional): Name to be shown as title of the plots. Defaults to 'u'. save(optional): Whether to save plot(s). Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ plt.figure(**figspecs) fen.plot(self.V.mesh()) if save: - plt.savefig(utilities.getNewFilename("msh", "eps"), - format='eps', dpi=1000) + plt.savefig(getNewFilename("msh", "eps"), format='eps', dpi=1000) plt.show() plt.close() -class FenicsHSAugmentedEngine(FenicsHSEngine): - """ - Fenics-based Hilbert space engine for augmented spaces. - """ - def __init__(self, V:FenSpace, d : int = 2): - self.d = d - FenicsHSEngine.__init__(self, V = V) - - def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> List[float]: - """ - Compute general norm of complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - normType(optional): Target norm identifier. If matrix, target norm - is one induced by normType. If number, target norm is weighted - H^1 norm with given weight. If string, must be recognizable by - Fenics norm command. Defaults to 'H10'. - - Returns: - Norms of the function (non-negative). - """ - if isinstance(normType, (int, float, str)): - u = utilities.listify(u, self.d) - norms = [None] * self.d - for j in range(self.d): - norms[j] = FenicsHSEngine.norm(self, u[j], normType) - return norms - elif type(normType).__name__[-6:] == "matrix": - u = utilities.vectorify(u, self.d) - if normType.shape[0] % u[0] == 0: - N = int(normType.shape[0] / u[0]) - else: - raise Exception("Energy matrix dimension mismatch.") - return FenicsHSEngine.norm(self, np.concatenate(u[: N]), normType) - raise Exception("Norm type not recognized.") - - def plot(self, u:Np1D, split : bool = True, name : str = "u", - save : bool = False, what : strLst = 'all', **figspecs): - """ - Do some nice plots of the complex-valued function with given dofs. - - Args: - u: numpy complex array with function dofs. - split: whether to split the given array. - name(optional): Name to be shown as title of the plots. Defaults to - 'u'. - save(optional): Whether to save plot(s). Defaults to False. - what(optional): Which plots to do. If list, can contain 'ABS', - 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. - Defaults to 'ALL'. - figspecs(optional key args): Optional arguments for matplotlib - figure creation. - """ - if split: - u = utilities.listify(u, self.d) - for j in range(self.d): - FenicsHSEngine.plot(self, u[j], what = what, - name = "{0}, comp. {1}".format(name, j), - **figspecs) - else: - FenicsHSEngine.plot(self, u, what = what, name = name, **figspecs) - diff --git a/rrompy/reduction_methods/__init__.py b/rrompy/reduction_methods/__init__.py new file mode 100644 index 0000000..892c0e6 --- /dev/null +++ b/rrompy/reduction_methods/__init__.py @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 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 . +# + diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py new file mode 100644 index 0000000..694f87c --- /dev/null +++ b/rrompy/reduction_methods/base/__init__.py @@ -0,0 +1,28 @@ +#!/usr/bin/python +# 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 + +__all__ = [ + 'GenericApproximant', + 'ParameterSweeper' + ] + + diff --git a/main/ROMApproximant.py b/rrompy/reduction_methods/base/generic_approximant.py similarity index 86% rename from main/ROMApproximant.py rename to rrompy/reduction_methods/base/generic_approximant.py index 4d1b9e4..9ece2f8 100644 --- a/main/ROMApproximant.py +++ b/rrompy/reduction_methods/base/generic_approximant.py @@ -1,328 +1,336 @@ #!/usr/bin/python +# 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 abc import abstractmethod import numpy as np -import utilities from copy import copy import scipy.sparse.linalg as spla -from RROMPyTypes import Np1D, Np2D, ListAny, DictAny, N2FSExpr, HFEng, HSEng,\ - Tuple, List +from rrompy.utilities.types import Np1D, Np2D, DictAny, N2FSExpr +from rrompy.utilities.types import HFEng, HSEng, Tuple, List +from rrompy.utilities import purgeDict -class ROMApproximant: +__all__ = ['GenericApproximant'] + +class GenericApproximant: """ + ABSTRACT ROM approximant computation for parametric problems. Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding k); - 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. - w(optional): Weight for norm computation. If None, set to Re(k0). - Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True. Defaults to empty dict. Attributes: HFEngine: HF problem solver. HSEngine: Hilbert space general purpose engine. mu0: Default parameter. w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + approximant. Recognized keys are in parameterList. + parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots. - initialHFData: HF problem initial data. 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, - w : float = None, approxParameters : DictAny = {}): + approxParameters : DictAny = {}): self.HFEngine = HFEngine self.HSEngine = HSEngine - self.initialHFData = self.HFEngine.problemData() + self._HFEngine0 = copy(HFEngine) + self.parameterList = ["POD"] + self.solLifting = self.HFEngine.liftDirichletData() self.mu0 = mu0 - if w is None: - w = np.real(self.mu0) - self.w = w self.approxParameters = approxParameters - self.energyNormMatrix = self.HFEngine.energyNormMatrix(self.w) + self.energyNormMatrix = self.HFEngine.energyNormMatrix() self.As, self.bs = self.HFEngine.getLSBlocks() def name(self) -> str: """Approximant label.""" return self.__class__.__name__ - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return ["POD"] - @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): if not hasattr(self, "approxParameters"): self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximant.parameterList(self), + approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters") keyList = list(approxParameters.keys()) if "POD" in keyList: self.POD = approxParameters["POD"] elif hasattr(self, "POD"): self.POD = self.POD else: self.POD = True + return approxParameters @property def POD(self): """Value of POD.""" return self._POD @POD.setter def POD(self, POD): if hasattr(self, "POD"): PODold = self.POD else: PODold = -1 self._POD = POD self._approxParameters["POD"] = self.POD if PODold != self.POD: self.resetSamples() def solveHF(self, mu : complex = None): """ Find high fidelity solution with original parameters and arbitrary parameter. Args: mu: Target parameter. """ if mu is None: mu = self.mu0 if (not hasattr(self, "lastSolvedHF") or not np.isclose(self.lastSolvedHF, mu)): A, b = self.buildLS(mu) self.uHF = spla.spsolve(A, b) self.lastSolvedHF = mu return self.uHF def buildLS(self, mu:complex, As : List[Np2D] = None, bs : List[Np1D] = None) -> Tuple[Np2D, Np1D]: """ Build linear system from polynomial blocks. Args: mu: Parameter value. As: Matrix blocks. If None, set to self.As. Defaults to None. bs: RHS blocks. If None, set to self.bs. Defaults to None. Returns: Matrix and RHS of system. """ if As is None: As = self.As if bs is None: bs = self.bs A = copy(As[0]) b = copy(bs[0]) for j in range(1, len(As)): A = A + np.power(mu, j) * As[j] for j in range(1, len(bs)): b = b + np.power(mu, j) * bs[j] return A, b @abstractmethod def resetSamples(self): """Reset samples. (ABSTRACT)""" pass @abstractmethod def setupApprox(self): """ Setup approximant. (ABSTRACT) Any specialization should include something like self.computeDerivatives() if not self.checkComputedApprox(): ... self.lastApproxParameters = copy(self.approxParameters) """ pass def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (hasattr(self, "lastApproxParameters") and self.approxParameters == self.lastApproxParameters) @abstractmethod def evalApprox(self, mu:complex) -> Np1D: """ Evaluate approximant at arbitrary parameter. (ABSTRACT) Any specialization should include something like self.setupApprox() self.uApp = ... Args: mu: Target parameter. Returns: Approximant as numpy complex vector. """ pass @abstractmethod def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ pass def HFNorm(self, mu:complex, normType : N2FSExpr = "H10") -> float: """ Compute norm of HF solution at arbitrary wavenumber. Args: mu: Target parameter. normType(optional): Target norm identifier. Must be recognizable by HSEngine norm command. If None, set to w. Defaults to None. Returns: Target norm of HFsolution. """ self.solveHF(mu) return self.HSEngine.norm(self.uHF, normType) def approxNorm(self, mu:complex, normType : N2FSExpr = "H10") -> float: """ Compute norm of (translated) approximant at arbitrary wavenumber. Args: mu: Target parameter. normType(optional): Target norm identifier. Must be recognizable by HSEngine norm command. If None, set to w. Defaults to None. Returns: Target norm of approximant. """ self.evalApprox(mu) return self.HSEngine.norm(self.uApp, normType) def approxError(self, mu:complex, normType : N2FSExpr = "H10") -> float: """ Compute norm of approximant at arbitrary wavenumber. Args: mu: Target parameter. normType(optional): Target norm identifier. Must be recognizable by HSEngine norm command. If None, set to w. Defaults to None. Returns: Target norm of (approximant - HFsolution). """ self.evalApprox(mu) self.solveHF(mu) self.evalApprox(mu) return self.HSEngine.norm(self.uApp - self.uHF, normType) def getHF(self, mu:complex) -> Np1D: """ Get HF solution at arbitrary wavenumber. Args: mu: Target parameter. Returns: HFsolution as numpy complex vector. """ self.solveHF(mu) return self.uHF def getApp(self, mu:complex) -> Np1D: """ Get approximant at arbitrary wavenumber. Args: mu: Target parameter. Returns: Approximant as numpy complex vector. """ self.evalApprox(mu) return self.uApp def plotHF(self, mu:complex, name : str = "u", **figspecs): """ Do some nice plots of the HF solution at arbitrary wavenumber. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ self.solveHF(mu) self.HSEngine.plot(self.uHF, name = name, **figspecs) def plotApp(self, mu:complex, name : str = "u", **figspecs): """ Do some nice plots of approximant at arbitrary wavenumber. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ self.evalApprox(mu) self.HSEngine.plot(self.uApp, name = name, **figspecs) def plotErr(self, mu:complex, name : str = "u", **figspecs): """ Do some nice plots of approximation error at arbitrary wavenumber. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ self.evalApprox(mu) self.solveHF(mu) self.HSEngine.plot(self.uApp - self.uHF, name = name, **figspecs) diff --git a/main/ROMApproximantSweeper.py b/rrompy/reduction_methods/base/parameter_sweeper.py similarity index 91% rename from main/ROMApproximantSweeper.py rename to rrompy/reduction_methods/base/parameter_sweeper.py index ac9b068..47d9b23 100644 --- a/main/ROMApproximantSweeper.py +++ b/rrompy/reduction_methods/base/parameter_sweeper.py @@ -1,275 +1,294 @@ #!/usr/bin/python +# 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 os import csv import warnings import numpy as np -from context import utilities -from RROMPyTypes import Np1D, DictAny, List, ROMEng, Tuple +from rrompy.utilities.types import Np1D, DictAny, List, ROMEng, Tuple +from rrompy.utilities import purgeList -class ROMApproximantSweeper: +__all__ = ['ParameterSweeper'] + +class ParameterSweeper: """ - ROM approximant sweeper. + ROM approximant parameter sweeper. Args: ROMEngine(optional): ROMApproximant class. Defaults to None. mutars(optional): Array of parameter values to sweep. Defaults to empty array. params(optional): List of parameter settings (each as a dict) to explore. Defaults to single empty set. mostExpensive(optional): String containing label of most expensive step, to be executed fewer times. Allowed options are 'HF' and 'Approx'. Defaults to 'HF'. Attributes: ROMEngine: ROMApproximant class. mutars: Array of parameter values to sweep. params: List of parameter settings (each as a dict) to explore. mostExpensive: String containing label of most expensive step, to be executed fewer times. """ allowedOutputs = ["HFNorm", "HFFunc", "AppNorm", "AppFunc", "ErrNorm", "ErrFunc"] def __init__(self, ROMEngine : ROMEng = None, mutars : Np1D = np.array([]), params : List[DictAny] = [{}], mostExpensive : str = "HF"): self.ROMEngine = ROMEngine self.mutars = mutars self.params = params self.mostExpensive = mostExpensive def name(self) -> str: """Approximant label.""" return self.__class__.__name__ @property def mostExpensive(self): """Value of mostExpensive.""" return self._mostExpensive @mostExpensive.setter def mostExpensive(self, mostExpensive:str): mostExpensive = mostExpensive.upper() if mostExpensive not in ["HF", "APPROX"]: warnings.warn(("Value of mostExpensive not recognized. Overriding " "to 'HF'."), stacklevel = 2) mostExpensive = "HF" self._mostExpensive = mostExpensive def checkValues(self) -> bool: """Check if sweep can be performed.""" if self.ROMEngine is None: warnings.warn("ROMEngine is missing. Aborting.", stacklevel = 2) return False if len(self.mutars) == 0: warnings.warn("Empty target parameter vector. Aborting.", stacklevel = 2) return False if len(self.params) == 0: warnings.warn("Empty method parameters vector. Aborting.", stacklevel = 2) return False return True def sweep(self, filename : str = "./out", outputs : List[str] = [], verbose : int = 1): if not self.checkValues(): return try: if outputs.upper() == "ALL": outputs = self.allowedOutputs + ["poles"] except: if len(outputs) == 0: outputs = self.allowedOutputs - outputs = utilities.purgeList(outputs, self.allowedOutputs + ["poles"], - listname = self.name() + ".outputs") + outputs = purgeList(outputs, self.allowedOutputs + ["poles"], + listname = self.name() + ".outputs") poles = ("poles" in outputs) if len(outputs) == 0: warnings.warn("Empty outputs. Aborting.", stacklevel = 2) return - outParList = self.ROMEngine.parameterList() + outParList = self.ROMEngine.parameterList Nparams = len(self.params) - allowedParams = self.ROMEngine.parameterList() + allowedParams = self.ROMEngine.parameterList while os.path.exists(filename): filename = filename + "{}".format(np.random.randint(10)) append_write = "w" initial_row = (outParList + ["muRe", "muIm"] + [x for x in self.allowedOutputs if x in outputs] + ["type"] + ["poles"] * poles) with open(filename, append_write, buffering = 1) as fout: writer = csv.writer(fout, delimiter=",") writer.writerow(initial_row) if self.mostExpensive == "HF": outerSet = self.mutars innerSet = self.params elif self.mostExpensive == "APPROX": outerSet = self.params innerSet = self.mutars for outerIdx, outerPar in enumerate(outerSet): if self.mostExpensive == "HF": i, mutar = outerIdx, outerPar elif self.mostExpensive == "APPROX": j, par = outerIdx, outerPar self.ROMEngine.approxParameters = {k: par[k] for k in\ par.keys() & allowedParams} self.ROMEngine.setupApprox() for innerIdx, innerPar in enumerate(innerSet): if self.mostExpensive == "APPROX": i, mutar = innerIdx, innerPar elif self.mostExpensive == "HF": j, par = innerIdx, innerPar self.ROMEngine.approxParameters = {k: par[k] for k in\ par.keys() & allowedParams} self.ROMEngine.setupApprox() if verbose >= 1: print("Set {}/{}\tmu_{} = {:.10f}{}"\ .format(j+1, Nparams, i, mutar, " " * 15), end="\r") outData = [] if "HFNorm" in outputs: val = self.ROMEngine.HFNorm(mutar) if isinstance(val, (list,)): val = val[0] outData = outData + [val] if "HFFunc" in outputs: outData = outData +[self.ROMEngine.HFEngine.functional( self.ROMEngine.solveHF(mutar))] if "AppNorm" in outputs: val = self.ROMEngine.approxNorm(mutar) if isinstance(val, (list,)): val = val[0] outData = outData + [val] if "AppFunc" in outputs: outData = outData +[self.ROMEngine.HFEngine.functional( self.ROMEngine.evalApprox(mutar))] if "ErrNorm" in outputs: val = self.ROMEngine.approxError(mutar) if isinstance(val, (list,)): val = val[0] outData = outData + [val] if "ErrFunc" in outputs: outData = outData +[self.ROMEngine.HFEngine.functional( self.ROMEngine.evalApprox(mutar) - self.ROMEngine.solveHF(mutar))] writeData = [] for parn in outParList: writeData = (writeData + [self.ROMEngine.approxParameters[parn]]) writeData = (writeData + [mutar.real, mutar.imag] + outData + [self.ROMEngine.name()]) if poles: writeData = writeData + list(self.ROMEngine.getPoles()) writer.writerow(str(x) for x in writeData) if verbose >= 1: if self.mostExpensive == "APPROX": print("Set {}/{}\tdone{}".format(j+1, Nparams, " "*25)) elif self.mostExpensive == "HF": print("Point mu_{} = {:.10f}\tdone{}".format(i, mutar, " " * 25)) self.filename = filename return self.filename def read(self, filename:str, restrictions : DictAny = {}, outputs : List[str] = []) -> DictAny: """ Execute a query on a custom format CSV. Args: filename: CSV filename. restrictions(optional): Parameter configurations to output. Defaults to empty dictionary, i.e. output all. outputs(optional): Values to output. Defaults to empty list, i.e. no output. Returns: Dictionary of desired results, with a key for each entry of outputs, and a numpy 1D array as corresponding value. """ def pairMuFromCSV(filename:str, muTarget:complex) -> Tuple[str, str]: """ Find complex point in CSV closer to a prescribed value. Args: filename: CSV filename. muTarget: Target complex value. Returns: Strings containing real and imaginary part of desired value, in the same format as in the CSV file. """ mutarsF = np.array([], dtype = complex) muRetarsF = np.array([], dtype = complex) muImtarsF = np.array([], dtype = complex) with open(filename, 'r') as f: reader = csv.reader(f, delimiter=',') header = next(reader) muReindex = header.index('muRe') muImindex = header.index('muIm') for row in reader: try: if row[muReindex] not in [" ", ""]: muRetarsF = np.append(muRetarsF, row[muReindex]) muImtarsF = np.append(muImtarsF, row[muImindex]) mutarsF = np.append(mutarsF, float(row[muReindex]) + 1.j * float(row[muImindex])) except: pass optimalIndex = np.argmin(np.abs(mutarsF - muTarget)) return [muRetarsF[optimalIndex], muImtarsF[optimalIndex]] with open(filename, 'r') as f: reader = csv.reader(f, delimiter=',') header = next(reader) restrIndices, outputIndices, outputData = {}, {}, {} for key in restrictions.keys(): try: restrIndices[key] = header.index(key) if not isinstance(restrictions[key], list): restrictions[key] = [restrictions[key]] restrictions[key] = [str(x) for x in restrictions[key]] except: warnings.warn("Ignoring key {} from restrictions"\ .format(key), stacklevel = 2) if 'muRe' in restrIndices.keys() or 'muIm' in restrIndices.keys(): if 'muRe' not in restrIndices.keys(): restrIndices['muRe'] = header.index('muRe') restrictions['muRe'] = [0.] * len(restrictions['muIm']) elif 'muIm' not in restrIndices.keys(): restrIndices['muIm'] = header.index('muIm') restrictions['muIm'] = [0.] * len(restrictions['muRe']) elif len(restrictions['muRe']) != len(restrictions['muIm']): raise Exception(("The lists of values for muRe and muIm " "must have the same length.")) for i in range(len(restrictions['muRe'])): mu = (1.0 * float(restrictions['muRe'][i]) + 1.j * float(restrictions['muIm'][i])) restrictions['muRe'][i], restrictions['muIm'][i] =\ pairMuFromCSV(filename, mu) for key in outputs: try: outputIndices[key] = header.index(key) outputData[key] = np.array([]) except: warnings.warn("Ignoring key {} from outputs".format(key), stacklevel = 2) for row in reader: if all([row[restrIndices[key]] in restrictions[key]\ for key in restrictions.keys()]): for key in outputIndices.keys(): try: val = row[outputIndices[key]] val = int(val) except: try: val = float(val) except: val = np.nan finally: outputData[key] = np.append(outputData[key], val) return outputData diff --git a/rrompy/reduction_methods/hermite/__init__.py b/rrompy/reduction_methods/hermite/__init__.py new file mode 100644 index 0000000..892c0e6 --- /dev/null +++ b/rrompy/reduction_methods/hermite/__init__.py @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 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 . +# + diff --git a/rrompy/reduction_methods/lagrange/__init__.py b/rrompy/reduction_methods/lagrange/__init__.py new file mode 100644 index 0000000..cb10b9b --- /dev/null +++ b/rrompy/reduction_methods/lagrange/__init__.py @@ -0,0 +1,30 @@ +#!/usr/bin/python +# 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.lagrange.generic_approximant_lagrange import GenericApproximantLagrange +from rrompy.reduction_methods.lagrange.approximant_lagrange_pade import ApproximantLagrangePade +from rrompy.reduction_methods.lagrange.approximant_lagrange_rb import ApproximantLagrangeRB + +__all__ = [ + 'GenericApproximantLagrange', + 'ApproximantLagrangePade', + 'ApproximantLagrangeRB' + ] + + diff --git a/main/ROMApproximantLagrangePade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py similarity index 84% rename from main/ROMApproximantLagrangePade.py rename to rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py index 55d543a..b3720f4 100644 --- a/main/ROMApproximantLagrangePade.py +++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py @@ -1,272 +1,273 @@ #!/usr/bin/python +# 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 utilities -from ROMApproximantLagrange import ROMApproximantLagrange -from RROMPyTypes import Np1D, ListAny +from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import GenericApproximantLagrange +from rrompy.utilities.types import Np1D, ListAny, DictAny, HSEng, HFEng -class ROMApproximantLagrangePade(ROMApproximantLagrange): +__all__ = ['ApproximantLagrangePade'] + +class ApproximantLagrangePade(GenericApproximantLagrange): """ ROM Lagrange Pade' interpolant computation for parametric problems. Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and mu to prescribed values; - getLSBlocks: get algebraic system blocks. HSEngine: Hilbert space general purpose engine. Should have members: - norm: compute Hilbert norm of expression represented as complex numpy vector; - plot: plot Hilbert expression represented as complex numpy vector. mus(optional): Boundaries of snapshot positions. Defaults to np.array([0, 1]). ws(optional): Array of snapshot weigths. Defaults to np.ones_like(self.mus). - w(optional): Weight for norm computation. If None, set to - Re(np.mean(self.mus)). Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; defaults to 2; - 'M': degree of Pade' interpolant numerator; defaults to 0; - 'N': degree of Pade' interpolant denominator; defaults to 0. Defaults to empty dict. plotSnap(optional): What to plot of snapshots of the Helmholtz solution map. See plot method in HSEngine. Defaults to []. plotSnapSpecs(optional): How to plot snapshots of the Helmholtz solution map. See plot method in HSEngine. Defaults to {}. Attributes: HFEngine: HF problem solver. HSEngine: Hilbert space general purpose engine. mu0: Default parameter. mus: Boundaries of snapshot positions. ws: Array of snapshot weigths. - w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + 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; - 'M': degree of Pade' interpolant numerator; - 'N': degree of Pade' interpolant denominator; extraApproxParameters: List of approxParameters keys in addition to mother class's. autoNodeTypes: List of possible values for autoNode. M: Numerator degree of approximant. N: Denominator degree of approximant. S: Number of solution snapshots over which current approximant is based upon. plotSnap: What to plot of snapshots of the Helmholtz solution map. plotSnapSpecs: How to plot snapshots of the Helmholtz solution map. solSnapshots: Array whose columns are FE dofs of snapshots of solution map. RPOD: Right factor of snapshots POD. If no POD, set to None. POD: Whether to compute POD of snapshots. autoNode: Type of nodes, if automatically generated. Otherwise None. Unused. 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. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. uApp: Last evaluated approximant as numpy complex vector. lastApproxParameters: List of parameters corresponding to last computed approximant. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix wrt mu. bs: List of numpy vectors representing coefficients of linear system RHS wrt mu. """ - extraApproxParameters = ["M", "N"] - - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximantLagrange.parameterList(self) - + ROMApproximantLagrangePade.extraApproxParameters) + def __init__(self, HFEngine:HFEng, HSEngine:HSEng, + mus : Np1D = np.array([0]), approxParameters : DictAny = {}, + plotSnap : ListAny = [], plotSnapSpecs : DictAny = {}): + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mus = mus, approxParameters = approxParameters, + plotSnap = plotSnap, plotSnapSpecs = plotSnapSpecs) + self.parameterList += ["M", "N"] @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change M, N and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantLagrangePade.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) keyList = list(approxParameters.keys()) - if "M" in keyList: - self.M = 0 #to avoid warnings if M and S are changed simultaneously - if "N" in keyList: - self.N = 0 #to avoid warnings if N and S are changed simultaneously - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantLagrangePade.extraApproxParameters, - True, True) - ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy) if "M" in keyList: self.M = approxParameters["M"] elif hasattr(self, "M"): self.M = self.M else: self.M = 0 if "N" in keyList: self.N = approxParameters["N"] elif hasattr(self, "N"): self.N = self.N else: self.N = 0 + return approxParameters @property def M(self): """Value of M. Its assignment may change S.""" return self._M @M.setter def M(self, M): if M < 0: raise ArithmeticError("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if hasattr(self, "S") and self.S < self.M + 1: warnings.warn("Prescribed S is too small. Updating S to M + 1.", stacklevel = 2) self.S = self.M + 1 @property def N(self): """Value of N. Its assignment may change S.""" return self._N @N.setter def N(self, N): if N < 0: raise ArithmeticError("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N if hasattr(self, "S") and self.S < self.N + 1: warnings.warn("Prescribed S is too small. Updating S to N + 1.", stacklevel = 2) self.S = self.N + 1 @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): if S <= 0: raise ArithmeticError("S must be positive.") if hasattr(self, "S"): Sold = self.S else: Sold = -1 vals, label = [0] * 2, {0:"M", 1:"N"} if hasattr(self, "M"): vals[0] = self.M if hasattr(self, "N"): vals[1] = self.N idxmax = np.argmax(vals) if vals[idxmax] + 1 > S: warnings.warn("Prescribed S is too small. Updating S to {} + 1."\ .format(label[idxmax]), stacklevel = 2) self.S = vals[idxmax] + 1 else: self._S = S self._approxParameters["S"] = self.S if Sold != self.S: self.resetSamples() def approxRadius(self) -> float: """Sample radius.""" return np.max(self.mu0 - self.mus) def setupApprox(self): """ Compute Pade' interpolant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): self.computeSnapshots() TN = np.vander(self.radiusPade(self.mus), N = self.N + 1, increasing = True) if self.POD: data = self.RPOD.T else: data = self.solSnapshots.T RHSFull = np.empty((TN.shape[0], data.shape[1] * TN.shape[1]), dtype = np.complex) for j in range(TN.shape[0]): RHSFull[j, :] = np.kron(TN[j, :], data[j, :]) G = np.polyfit(self.radiusPade(self.mus), RHSFull, deg = self.N, w = self.ws)[0, :].reshape((TN.shape[1], data.shape[1])).T if self.POD: _, _, V = np.linalg.svd(G, full_matrices = False) self.Q = V[-1, :].conj() else: G2 = G.conj().T.dot(self.energyNormMatrix.dot(G)) _, eV = np.linalg.eigh(G2) self.Q = eV[:, 0] TNQ = TN.dot(self.Q).reshape((self.S, 1)) if self.POD: data = self.solSnapshots.dot(self.RPOD) else: data = self.solSnapshots RHS = np.multiply(data.T, TNQ) self.P = np.polyfit(self.radiusPade(self.mus), RHS, deg = self.M, w = self.ws)[::-1, :].T self.approxParameters = {"N" : self.N, "M" : self.M, "S" : self.S} self.lastApproxParameters = copy(self.approxParameters) def radiusPade(self, mu:Np1D, mu0 : float = None) -> float: """ Compute translated radius to be plugged into Pade' approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Translated radius to be plugged into Pade' approximant. """ if mu0 is None: mu0 = self.mu0 return (mu - mu0) / self.approxRadius() def evalApprox(self, mu:complex) -> Np1D: """ Evaluate Pade' approximant at arbitrary wavenumber. Args: mu: Target parameter. Returns: Real and imaginary parts of approximant. """ self.setupApprox() powerlist = np.power(self.radiusPade(mu), range(max(self.M, self.N) + 1)) self.uApp = (self.P.dot(powerlist[:self.M + 1]) / self.Q.dot(powerlist[:self.N + 1])) return self.uApp def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() return np.roots(self.Q[::-1]) * self.approxRadius() + self.mu0 diff --git a/main/ROMApproximantLagrangeRB.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py similarity index 76% rename from main/ROMApproximantLagrangeRB.py rename to rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py index 576335c..2a44d59 100644 --- a/main/ROMApproximantLagrangeRB.py +++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py @@ -1,272 +1,265 @@ #!/usr/bin/python +# 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 -import utilities -from ROMApproximantLagrange import ROMApproximantLagrange -from RROMPyTypes import Np1D, ListAny, DictAny, HFEng, HSEng +from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import GenericApproximantLagrange +from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng -class ROMApproximantLagrangeRB(ROMApproximantLagrange): +__all__ = ['ApproximantLagrangeRB'] + +class ApproximantLagrangeRB(GenericApproximantLagrange): """ ROM RB approximant computation for parametric problems. Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and mu to prescribed values; - getLSBlocks: get algebraic system blocks; - 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. mus(optional): Array of snapshot parameters. Defaults to np.array([0]). ws(optional): Array of snapshot weigths (unused). Defaults to np.ones_like(self.mus). - w(optional): Weight for norm computation. If None, set to - Re(np.mean(self.mus)). Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; defaults to 2; - '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. - w: Weight for norm computation. 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: + 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; - 'R': rank for Galerkin projection. extraApproxParameters: List of approxParameters keys in addition to mother class's. autoNodeTypes: List of possible values for autoNode. S: Number of solution snapshots over which current approximant is based upon. 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. autoNode: Type of nodes, if automatically generated. Otherwise None. Unused. 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. """ - extraApproxParameters = ["R"] - - def __init__(self, HFEngine:HFEng, HSEngine:HSEng, - mus : Np1D = np.array([0]), w : float = None, - approxParameters : DictAny = {}, plotSnap : ListAny = [], - plotSnapSpecs : DictAny = {}): - ROMApproximantLagrange.__init__(self, HFEngine = HFEngine, - HSEngine = HSEngine, mus = mus, w = w, - approxParameters = approxParameters, - plotSnap = plotSnap, - plotSnapSpecs = plotSnapSpecs) - self.solLifting = self.HFEngine.liftDirichletData() + def __init__(self, HFEngine:HFEng, HSEngine:HSEng, + mus : Np1D = np.array([0]), approxParameters : DictAny = {}, + plotSnap : ListAny = [], plotSnapSpecs : DictAny = {}): + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mus = mus, approxParameters = approxParameters, + plotSnap = plotSnap, plotSnapSpecs = plotSnapSpecs) + self.parameterList += ["R"] def resetSamples(self): """Reset samples.""" - ROMApproximantLagrange.resetSamples(self) + super().resetSamples(self) self.projMat = None - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximantLagrange.parameterList(self) - + ROMApproximantLagrangeRB.extraApproxParameters) - @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change M, N and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantLagrangeRB.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) keyList = list(approxParameters.keys()) - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantLagrangeRB.extraApproxParameters, - True, True) - ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy) if "R" in keyList: self.R = approxParameters["R"] elif hasattr(self, "R"): self.R = self.R else: self.R = self.S + return approxParameters @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 manageSnapshots(self, u:ListAny, pos:int) -> Np1D: - """ - Post-process snapshots of solution map. - - Args: - u: Numpy 1D array of FE dofs of snapshot. - pos: Derivative index. - - Returns: - Numpy 1D array with adjusted snapshot dofs. - """ - return ROMApproximantLagrange.manageSnapshots(self, - u - self.solLifting, pos) +# def manageSnapshots(self, u:ListAny, pos:int) -> Np1D: +# """ +# Post-process snapshots of solution map. +# +# Args: +# u: Numpy 1D array of FE dofs of snapshot. +# pos: Derivative index. +# +# Returns: +# Numpy 1D array with adjusted snapshot dofs. +# """ +# return super().manageSnapshots(self, u - self.solLifting, 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.projMat is not None - and ROMApproximantLagrange.checkComputedApprox(self)) + and super().checkComputedApprox(self)) 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] * max(len(self.As), len(self.bs)) for j in range(len(self.As)): self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat)) if j < len(self.bs): self.bRBs[j] = projMatH.dot(self.bs[j] - self.As[j].dot(self.solLifting)) else: self.bRBs[j] = - projMatH.dot(self.As[j].dot(self.solLifting)) for j in range(len(self.As), 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.solLifting + self.solveReducedSystem(mu) +# self.uApp = self.solLifting + self.solveReducedSystem(mu) + 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/main/ROMApproximantLagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py similarity index 85% rename from main/ROMApproximantLagrange.py rename to rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py index 0e3d6a8..455756b 100644 --- a/main/ROMApproximantLagrange.py +++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py @@ -1,284 +1,282 @@ #!/usr/bin/python +# 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 -import utilities import warnings -from ROMApproximant import ROMApproximant -from RROMPyTypes import Np1D, ListAny, DictAny, HFEng, HSEng +from rrompy.reduction_methods.base.generic_approximant import GenericApproximant +from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng -class ROMApproximantLagrange(ROMApproximant): +__all__ = ['GenericApproximantLagrange'] + +class GenericApproximantLagrange(GenericApproximant): """ ROM Lagrange interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and 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. mus(optional): Array of snapshot parameters. Defaults to np.array([0]). ws(optional): Array of snapshot weigths (possibly unused). Defaults to np.ones_like(self.mus). - w(optional): Weight for norm computation. If None, set to - Re(np.mean(self.mus)). Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon. Defaults to 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 (possibly unused). - w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + 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. extraApproxParameters: List of approxParameters keys in addition to mother class's. autoNodeTypes: List of possible values for autoNode. S: Number of solution snapshots over which current approximant is based upon. plotSnap: What to plot of snapshots of the Helmholtz solution map. plotSnapSpecs: How to plot snapshots of the Helmholtz solution map. solSnapshots: Array whose columns are FE dofs of snapshots of solution map. RPOD: Right factor of snapshots POD. If no POD, set to None. POD: Whether to compute POD of snapshots. autoNode: Type of nodes, if automatically generated. Otherwise None. Possibly unused. 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. """ - extraApproxParameters = ["S"] autoNodeTypes = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"] def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mus : Np1D = np.array([0]), ws : Np1D = None, - w : float = None, approxParameters : DictAny = {}, - plotSnap : ListAny = [], plotSnapSpecs : DictAny = {}): + approxParameters : DictAny = {}, plotSnap : ListAny = [], + plotSnapSpecs : DictAny = {}): self.mus = mus if ws is None: ws = np.ones_like(self.mus) self.ws = ws - ROMApproximant.__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, - mu0 = np.mean(mus), w = w, - approxParameters = approxParameters) + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mu0 = self.mu0, approxParameters = approxParameters) + self.parameterList += ["S"] self.plotSnap = plotSnap self.plotSnapSpecs = plotSnapSpecs self.resetSamples() - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximant.parameterList(self) - + ROMApproximantLagrange.extraApproxParameters) - @property def mu0(self): """Dummy center of approximant (i.e. mu0).""" self.mu0 = np.mean(self.mus) return self._mu0 @mu0.setter def mu0(self, mu0:bool): self._mu0 = mu0 @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 ws(self): """Value of ws.""" return self._mus @ws.setter def ws(self, ws): if hasattr(self, 'ws'): wsOld = self.ws else: wsOld = None ws = np.array(ws) self._ws = np.abs(ws) if not np.array_equal(self._ws, ws): warnings.warn(("Negative weights not allowed. Overriding to " "absolute values."), stacklevel = 2) if (wsOld is None or len(self.ws) != len(wsOld) or not np.allclose(self.ws, wsOld, 1e-14)): self.autoNode = None def checkNodesWeightsConsistency(self) -> bool: """Check if mus and ws have consistent lengths.""" return self.mus.shape == self.ws.shape @property def approxParameters(self): """Value of approximant parameters. Its assignment may change S.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantLagrange.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) keyList = list(approxParameters.keys()) - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantLagrange.extraApproxParameters, - True, True) - ROMApproximant.approxParameters.fset(self, approxParametersCopy) if "S" in keyList: self.S = approxParameters["S"] elif hasattr(self, "S"): self.S = self.S else: self.S = 2 + return approxParameters @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() def resetSamples(self): """Reset samples.""" self.solSnapshots = None self.RPOD = None def setNodesWeights(self, xs : Np1D = np.array([0, 1]), nodeTypes : str = "UNIFORM"): """ Assign sampling nodes and weights based on a quadrature formula. Args: xs: extremes of the sampling interval; nodeTypes: type of nodes/weights; allowed values are in self.autoNodeTypes. """ if len(xs) != 2: raise Exception(("Value of xs not recognized. Should be numpy 1D " "vector or list with length 2.")) try: nodeTypes = nodeTypes.upper().strip().replace(" ","") if nodeTypes not in self.autoNodeTypes: raise except: warnings.warn(("Prescribed nodeTypes not recognized. Overriding " "to 'UNIFORM'."), stacklevel = 2) nodeTypes = "UNIFORM" if nodeTypes == "UNIFORM": nodes = np.linspace(1., -1., self.S) weights = np.ones((self.S, 1)) elif nodeTypes == "CHEBYSHEV": nodes, weights = np.polynomial.chebyshev.chebgauss(self.S) weights = weights[:, None] * 2. / np.pi else: #if nodeTypes == "GAUSSLEGENDRE": nodes, weights = np.polynomial.legendre.leggauss(self.S) weights = weights[:, None] mu0 = np.mean(xs) r = (xs[0] - xs[1]) / 2. self.mus = r * nodes + mu0 self.ws = np.abs(r) * weights self.autoNode = nodeTypes def computeSnapshots(self): """ Compute snapshots of solution map. """ if self.solSnapshots is None: if len(self.mus) != self.S: warnings.warn(("Number of prescribed nodes different from " "S. Overriding S to len(mus)."), stacklevel = 2) self.S = len(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.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] 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 ROMApproximant.checkComputedApprox(self)) + and super().checkComputedApprox(self)) diff --git a/rrompy/reduction_methods/taylor/__init__.py b/rrompy/reduction_methods/taylor/__init__.py new file mode 100644 index 0000000..5e416ae --- /dev/null +++ b/rrompy/reduction_methods/taylor/__init__.py @@ -0,0 +1,30 @@ +#!/usr/bin/python +# 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.taylor.generic_approximant_taylor import GenericApproximantTaylor +from rrompy.reduction_methods.taylor.approximant_taylor_pade import ApproximantTaylorPade +from rrompy.reduction_methods.taylor.approximant_taylor_rb import ApproximantTaylorRB + +__all__ = [ + 'GenericApproximantTaylor', + 'ApproximantTaylorPade', + 'ApproximantTaylorRB' + ] + + diff --git a/main/ROMApproximantTaylorPade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py similarity index 90% rename from main/ROMApproximantTaylorPade.py rename to rrompy/reduction_methods/taylor/approximant_taylor_pade.py index caf18cf..1b0a3dc 100644 --- a/main/ROMApproximantTaylorPade.py +++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py @@ -1,450 +1,450 @@ #!/usr/bin/python +# 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 utilities -from ROMApproximantTaylor import ROMApproximantTaylor -from RROMPyTypes import Np1D, Np2D, ListAny, Tuple +from rrompy.reduction_methods.taylor.generic_approximant_taylor import GenericApproximantTaylor +from rrompy.utilities.types import Np1D, Np2D, ListAny, Tuple, DictAny +from rrompy.utilities.types import HFEng, HSEng -class ROMApproximantTaylorPade(ROMApproximantTaylor): +__all__ = ['ApproximantTaylorPade'] + +class ApproximantTaylorPade(GenericApproximantTaylor): """ ROM single-point fast Pade' approximant computation for parametric problems. Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and mu to prescribed values; - 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. - w(optional): Weight for norm computation. If None, set to Re(self.mu0). - Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'rho': weight for computation of original Pade' approximant; defaults to 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. - w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + 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. - extraApproxParameters: List of approxParameters keys in addition to - mother class's. 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. """ - extraApproxParameters = ["M", "N", "robustTol"]#, "rho"] - - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. + def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0, + approxParameters : DictAny = {}, plotDer : ListAny = [], + plotDerSpecs : DictAny = {}): + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mu0 = mu0, approxParameters = approxParameters, + plotDer = plotDer, plotDerSpecs = plotDerSpecs) + self.parameterList += ["M", "N", "robustTol", "rho"] - Returns: - List of strings. - """ - return (ROMApproximantTaylor.parameterList(self) - + ROMApproximantTaylorPade.extraApproxParameters) - @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantTaylorPade.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) 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 - if "M" in keyList: - self.M = 0 #to avoid warnings if M and E are changed simultaneously - if "N" in keyList: - self.N = 0 #to avoid warnings if N and E are changed simultaneously - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantTaylorPade.extraApproxParameters, - True, True) - ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy) if "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 + return approxParameters @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]) 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/main/ROMApproximantTaylorRB.py b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py similarity index 82% rename from main/ROMApproximantTaylorRB.py rename to rrompy/reduction_methods/taylor/approximant_taylor_rb.py index 0708116..c2d2ed2 100644 --- a/main/ROMApproximantTaylorRB.py +++ b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py @@ -1,331 +1,313 @@ #!/usr/bin/python +# 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 -import utilities -from ROMApproximantTaylor import ROMApproximantTaylor -from RROMPyTypes import Np1D, ListAny +from rrompy.reduction_methods.taylor.generic_approximant_taylor import GenericApproximantTaylor +from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng -class ROMApproximantTaylorRB(ROMApproximantTaylor): +__all__ = ['ApproximantTaylorRB'] + +class ApproximantTaylorRB(GenericApproximantTaylor): """ ROM single-point fast RB approximant computation for parametric problems with polynomial dependence up to degree 2. Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: sparse matrix (in CSC format) associated to w-weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and mu to prescribed values; - getLSBlocks: get blocks of HF linear system as sparse matrices in CSC format; - liftDirichletData: perform lifting of Dirichlet BC as numpy complex vector; - 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. - w(optional): Weight for norm computation. If None, set to Re(self.mu0). - Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'R': rank for Galerkin projection; defaults to E + 1; - '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. - w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + approximant. Recognized keys are in parameterList. + parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'R': rank for Galerkin projection; - '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. - extraApproxParameters: List of approxParameters keys in addition to - mother class's. 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. R: Rank for Galerkin projection. 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, i.e. 'KRYLOV'. plotDer: What to plot of derivatives of the Helmholtz solution map. plotDerSpecs: How to plot derivatives of the Helmholtz solution map. energyNormMatrix: Sparse matrix (in CSC format) associated 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. uApp: Last evaluated approximant as numpy complex vector. lastApproxParameters: List of parameters corresponding to last computed approximant. solLifting: Numpy complex vector with lifting of real part of Dirichlet boundary datum. projMat: Numpy matrix representing projection onto RB space. projMat: Numpy matrix representing projection onto RB space. 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 RB coefficients of linear system matrix wrt mu. bRBs: List of numpy vectors representing RB coefficients of linear system RHS wrt mu. """ - extraApproxParameters = ["R"] + def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0, + approxParameters : DictAny = {}, plotDer : ListAny = [], + plotDerSpecs : DictAny = {}): + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mu0 = mu0, approxParameters = approxParameters, + plotDer = plotDer, plotDerSpecs = plotDerSpecs) + self.parameterList += ["R"] def resetSamples(self): """Reset samples.""" - ROMApproximantTaylor.resetSamples(self) + super().resetSamples(self) self.projMat = None - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximantTaylor.parameterList(self) - + ROMApproximantTaylorRB.extraApproxParameters) - @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change M, N and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantTaylorRB.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) keyList = list(approxParameters.keys()) - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantTaylorRB.extraApproxParameters, - True, True) - ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy) if "R" in keyList: self.R = approxParameters["R"] elif hasattr(self, "R"): self.R = self.R else: self.R = self.E + 1 + return approxParameters @property def POD(self): """Value of POD.""" return self._POD @POD.setter def POD(self, POD): - ROMApproximantTaylor.POD.fset(self, POD) + super().POD.fset(self, POD) if (hasattr(self, "sampleType") and self.sampleType == "ARNOLDI" and not self.POD): warnings.warn(("Arnoldi sampling implicitly forces POD-type " "derivative management."), stacklevel = 2) @property def sampleType(self): """Value of sampleType.""" return self._sampleType @sampleType.setter def sampleType(self, sampleType): - ROMApproximantTaylor.sampleType.fset(self, sampleType) + super().sampleType.fset(self, sampleType) if (hasattr(self, "POD") and not self.POD and self.sampleType == "ARNOLDI"): warnings.warn(("Arnoldi sampling implicitly forces POD-type " "derivative management."), stacklevel = 2) @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, "E") and self.E + 1 < self.R: warnings.warn("Prescribed E is too small. Updating E to R - 1.", stacklevel = 2) self.E = self.R - 1 - 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: - u = u - self.solLifting - return ROMApproximantTaylor.manageDerivatives(self, u, pos) - def setupApprox(self): """ Setup RB system. For usage of correlation matrix in POD see Section 6.3.1 in 'Reduced Basis Methods for PDEs. An introduction', A. Quarteroni, A. Manzoni, F. Negri, Springer, 2016. """ if not self.checkComputedApprox(): self.computeDerivatives() if self.POD and not self.sampleType == "ARNOLDI": A = copy(self.solDerivatives) M = self.energyNormMatrix V = np.zeros_like(A, dtype = np.complex) Q = 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]): Q[:, k] = (np.random.randn(Q.shape[0]) + 1.j * np.random.randn(Q.shape[0])) if k > 0: for times in range(2): #twice is enough! Q[:, k] = Q[:, k] - Q[:, :k].dot( (Q[:, k].conj().dot(M.dot(Q[:, :k]))).conj()) Q[:, k] = Q[:, k]/(Q[:, k].conj().dot(M.dot(Q[:, k])))**.5 R[k, k] = np.abs(A[:, k].conj().dot(M.dot(A[:, k]))) ** .5 alpha = Q[:, k].conj().dot(M.dot(A[:, k])) if np.isclose(np.abs(alpha), 0.): s = 1. else: s = - alpha / np.abs(alpha) Q[:, k] = s * Q[:, k] v = R[k, k] * Q[:, k] - A[:, k] for times in range(2): #twice is enough! v = v - Q[:, :k].dot((v.conj().dot(M.dot(Q[:, :k])) ).conj()) sigma = np.abs(v.conj().dot(M.dot(v))) ** .5 if np.isclose(sigma, 0.): v = Q[:, k] else: v = v / sigma V[:, k] = v 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] = Q[:, k].conj().dot(M.dot(A[:, J])) A[:, J] = A[:, J] - np.outer(Q[:, k], R[k, J]) for k in range(A.shape[1] - 1, -1, -1): v = V[:, k] J = np.arange(k, A.shape[1]) vtQ = v.conj().dot(M.dot(Q[:, J])) Q[:, J] = Q[:, J] - 2 * np.outer(v, vtQ) self.projMatQ = Q self.projMatR = R if self.POD and not self.sampleType == "ARNOLDI": U, _, _ = np.linalg.svd(self.projMatR[: self.R, : self.R]) self.projMat = self.projMatQ[:, : self.R].dot(U) else: self.projMat = self.solDerivatives[:, : self.R] self.assembleReducedSystem() self.lastApproxParameters = copy(self.approxParameters) def assembleReducedSystem(self): """Build affine blocks of RB linear system through projections.""" projMatH = self.projMat.T.conjugate() self.ARBs = [None] * len(self.As) - self.bRBs = [None] * max(len(self.As), len(self.bs)) + self.bRBs = [None] * len(self.bs) for j in range(len(self.As)): self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat)) - if j < len(self.bs): - self.bRBs[j] = projMatH.dot(self.bs[j] - - self.As[j].dot(self.solLifting)) - else: - self.bRBs[j] = - projMatH.dot(self.As[j].dot(self.solLifting)) - for j in range(len(self.As), len(self.bs)): + 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: Approximant as numpy 1D array. """ self.setupApprox() - self.uApp = self.solLifting + self.solveReducedSystem(mu) + 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/main/ROMApproximantTaylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py similarity index 85% rename from main/ROMApproximantTaylor.py rename to rrompy/reduction_methods/taylor/generic_approximant_taylor.py index d839996..d6b9b61 100644 --- a/main/ROMApproximantTaylor.py +++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py @@ -1,304 +1,301 @@ #!/usr/bin/python +# 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 utilities import scipy.sparse.linalg as spla -from ROMApproximant import ROMApproximant -from RROMPyTypes import Np1D, ListAny, DictAny, HFEng, HSEng +from rrompy.reduction_methods.base.generic_approximant import GenericApproximant +from rrompy.utilities.types import Np1D, ListAny, DictAny, HFEng, HSEng -class ROMApproximantTaylor(ROMApproximant): +__all__ = ['GenericApproximantTaylor'] + +class GenericApproximantTaylor(GenericApproximant): """ ROM single-point approximant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. Should have members: - - Vdim: domain dimension; - energyNormMatrix: assemble sparse matrix (in CSC format) associated to weighted H10 norm; - problemData: list of HF problem data (excluding mu); - setProblemData: set HF problem data and mu to prescribed values; - getLSBlocks: get algebraic system blocks. HSEngine: Hilbert space general purpose engine. Should have members: - norm: compute Hilbert norm of expression represented as complex numpy vector; - plot: plot Hilbert expression represented as complex numpy vector. mu0(optional): Default parameter. Defaults to 0. - w(optional): Weight for norm computation. If None, set to Re(self.mu0). - Defaults to None. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - '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. - w: Weight for norm computation. approxParameters: Dictionary containing values for main parameters of - approximant. Recognized keys are: + 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. - extraApproxParameters: List of approxParameters keys in addition to - mother class's. 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. """ - extraApproxParameters = ["E", "Emax", "sampleType"] - def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0, - w : float = None, approxParameters : DictAny = {}, - plotDer : ListAny = [], plotDerSpecs : DictAny = {}): - ROMApproximant.__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, - mu0 = mu0, w = w, - approxParameters = approxParameters) - self.solLifting = self.HFEngine.liftDirichletData() + approxParameters : DictAny = {}, plotDer : ListAny = [], + plotDerSpecs : DictAny = {}): + super().__init__(self, HFEngine = HFEngine, HSEngine = HSEngine, + mu0 = mu0, approxParameters = approxParameters) + self.parameterList += ["E", "Emax", "sampleType"] 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 - def parameterList(self) -> ListAny: - """ - List containing keys which are allowed in approxParameters. - - Returns: - List of strings. - """ - return (ROMApproximant.parameterList(self) - + ROMApproximantTaylor.extraApproxParameters) - @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change E and Emax. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): - if not hasattr(self, "approxParameters"): - self._approxParameters = {} - approxParameters = utilities.purgeDict(approxParams, - ROMApproximantTaylor.parameterList(self), - dictname = self.name() + ".approxParameters") + approxParameters = super().approxParameters.fset(self, approxParams) keyList = list(approxParameters.keys()) - approxParametersCopy = utilities.purgeDict(approxParameters, - ROMApproximantTaylor.extraApproxParameters, - True, True) - ROMApproximant.approxParameters.fset(self, approxParametersCopy) 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" + return approxParameters @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.parDegree > 1: + 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 = [np.zeros(self.HFEngine.Vdim())] * (lAs - 1) + uOlds = None for j in range(self.Emax + 1): if j == 0: - self.HFEngine.setProblemData(self.initialHFData, self.mu0) + 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) - del uOlds[0] + 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.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]) 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 ROMApproximant.checkComputedApprox(self)) diff --git a/rrompy/sampling/__init__.py b/rrompy/sampling/__init__.py new file mode 100644 index 0000000..fab6df7 --- /dev/null +++ b/rrompy/sampling/__init__.py @@ -0,0 +1,30 @@ +#!/usr/bin/python +# 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.sampling.generic_sampler import GenericSampler +from rrompy.sampling.quadrature_sampler import QuadratureSampler +from rrompy.sampling.random_sampler import RandomSampler + +__all__ = [ + 'GenericSampler', + 'QuadratureSampler', + 'RandomSampler' + ] + + diff --git a/rrompy/sampling/generic_sampler.py b/rrompy/sampling/generic_sampler.py new file mode 100644 index 0000000..4858b03 --- /dev/null +++ b/rrompy/sampling/generic_sampler.py @@ -0,0 +1,52 @@ +#!/usr/bin/python +# 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 abc import abstractmethod +import numpy as np +from rrompy.utilities.types import Np1D, Tuple, GenExpr + +__all__ = ['GenericSampler'] + +class GenericSampler: + """ABSTRACT. Generic generator of sample points.""" + + def __init__(self, lims:Tuple[Np1D, Np1D]): + self.lims = lims + + @property + def lims(self): + """Value of lims.""" + return self._lims + @lims.setter + def lims(self, lims): + if len(lims) != 2: + raise Exception("2 limits must be specified.") + for j in range(2): + try: + len(lims[j]) + except: + lims[j] = np.array([lims[j]]) + if len(lims[0]) != len(lims[1]): + raise Exception("The limits must have the same length.") + self._lims = lims + + @abstractmethod + def generatePoints(self, n:GenExpr) -> Tuple[Np1D, Np1D]: + """Array of points and array of weights.""" + pass diff --git a/rrompy/sampling/quadrature_sampler.py b/rrompy/sampling/quadrature_sampler.py new file mode 100644 index 0000000..b4761ce --- /dev/null +++ b/rrompy/sampling/quadrature_sampler.py @@ -0,0 +1,80 @@ +#!/usr/bin/python +# 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.sampling.generic_sampler import GenericSampler +from rrompy.utilities.types import Np1D, Tuple, List + +__all__ = ['QuadratureSampler'] + +class QuadratureSampler(GenericSampler): + """Generator of quadrature sample points.""" + + allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"] + + def __init__(self, lims:Tuple[Np1D, Np1D], kind:str): + super().__init__(self, lims = lims) + self.kind = kind + + @property + def kind(self): + """Value of kind.""" + return self._kind + @kind.setter + def kind(self, kind): + if kind.upper() not in self.allowedKinds: + raise Exception("Generator kind not recognized.") + self._kind = kind.upper() + + def generatePoints(self, n:Np1D) -> List[Np1D]: + """Array of quadrature points and array of weights.""" + d = len(self.lims[0]) + try: + len(n) + except: + n = np.array([n]) + if len(n) != d: + raise Exception(("Numbers of points must have same dimension as" + "limits.")) + for j in range(d): + a, b = self.lims[0][j], self.lims[1][j] + if self.kind == "UNIFORM": + xj = np.linspace(a, b, n[j])[:, None] + wj = np.abs(a - b) / n[j] * np.ones(n[j]) + elif self.kind == "CHEBYSHEV": + nodes, weights = np.polynomial.chebyshev.chebgauss(n[j]) + xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None] + wj = np.abs(a - b) / np.pi * weights + elif self.kind == "GAUSSLEGENDRE": + nodes, weights = np.polynomial.legendre.leggauss(n[j]) + xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None] + wj = np.abs(a - b) * weights + if j == 0: + x = xj + w = wj + xsize = n[0] + else: + x = np.concatenate((np.kron(np.ones(n[j])[:, None], x), + np.kron(xj, np.ones(xsize)[:, None])), + axis = 1) + w = np.multiply(np.kron(np.ones(n[j]), w), + np.kron(wj, np.ones(xsize))) + xsize = xsize * n[j] + return [y.flatten() for y in np.split(x, xsize)], w + diff --git a/rrompy/sampling/random_sampler.py b/rrompy/sampling/random_sampler.py new file mode 100644 index 0000000..50151dc --- /dev/null +++ b/rrompy/sampling/random_sampler.py @@ -0,0 +1,58 @@ +#!/usr/bin/python +# 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.utilities.sobol import sobolGenerate +from rrompy.sampling.generic_sampler import GenericSampler +from rrompy.utilities.types import Np1D, Tuple, List + +__all__ = ['RandomSampler'] + +class RandomSampler(GenericSampler): + """Generator of quadrature sample points.""" + + allowedKinds = ["UNIFORM", "SOBOL"] + + def __init__(self, lims:Tuple[Np1D, Np1D], kind:str): + super().__init__(self, lims = lims) + self.kind = kind + + @property + def kind(self): + """Value of kind.""" + return self._kind + @kind.setter + def kind(self, kind): + if kind.upper() not in self.allowedKinds: + raise Exception("Generator kind not recognized.") + self._kind = kind.upper() + + def generatePoints(self, n:int, seed : int = 0) -> List[Np1D]: + """Array of quadrature points and array of weights.""" + d = len(self.lims[0]) + if self.kind == "UNIFORM": + np.random.seed(seed) + x = np.random.uniform(size = (n, d)) + else: + x = sobolGenerate(d, n, seed) + for j in range(d): + a, b = self.lims[0][j], self.lims[1][j] + x[:, j] = a + (b - a) * x[:, j] + return [y.flatten() for y in np.split(x, n)], np.ones(n) + diff --git a/rrompy/utilities/__init__.py b/rrompy/utilities/__init__.py new file mode 100644 index 0000000..3c0e23a --- /dev/null +++ b/rrompy/utilities/__init__.py @@ -0,0 +1,42 @@ +#!/usr/bin/python +# 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.utilities.find_dict_str_key import findDictStrKey +from rrompy.utilities.get_new_filename import getNewFilename +from rrompy.utilities.purge_dict import purgeDict +from rrompy.utilities.purge_list import purgeList +from rrompy.utilities.number_theory import squareResonances, primeFactorize, getLowestPrimeFactor +from rrompy.utilities.sobol import sobolGenerate +from rrompy.utilities import types as Types +from rrompy.utilities import fenics as Fenics + +__all__ = [ + 'findDictStrKey', + 'getNewFilename', + 'purgeDict', + 'purgeList', + 'squareResonances', + 'primeFactorize', + 'getLowestPrimeFactor', + 'sobolGenerate', + 'Types', + 'Fenics' + ] + + diff --git a/rrompy/utilities/fenics.py b/rrompy/utilities/fenics.py new file mode 100644 index 0000000..9eab226 --- /dev/null +++ b/rrompy/utilities/fenics.py @@ -0,0 +1,31 @@ +#!/usr/bin/python +# 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 fenics as fen + +__all__ = ['fenZERO','fenONE','bdrTrue','bdrFalse'] + +# CONSTANTS +fenZERO = fen.Constant(0.) +fenONE = fen.Constant(1.) + +# BOUNDARY BOOLS +bdrTrue = lambda x, on_boundary : on_boundary +bdrFalse = lambda x, on_boundary : False + diff --git a/rrompy/utilities/find_dict_str_key.py b/rrompy/utilities/find_dict_str_key.py new file mode 100644 index 0000000..5a5666d --- /dev/null +++ b/rrompy/utilities/find_dict_str_key.py @@ -0,0 +1,29 @@ +#!/usr/bin/python +# 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.utilities.types import Any, ListAny + +__all__ = ['findDictStrKey'] + +def findDictStrKey(key:Any, keyList:ListAny): + for akey in keyList: + if isinstance(key, str) and key.lower() == akey.lower(): + return akey + return None + diff --git a/rrompy/utilities/get_new_filename.py b/rrompy/utilities/get_new_filename.py new file mode 100644 index 0000000..9fc6aa5 --- /dev/null +++ b/rrompy/utilities/get_new_filename.py @@ -0,0 +1,31 @@ +#!/usr/bin/python +# 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 os +import numpy as np + +__all__ = ['getNewFilename'] + +def getNewFilename(prefix : str = "", extension : str = "dat") -> str: + n = len(extension) + 1 + filename = "{}{}.{}".format(prefix, np.random.randint(0, 10), extension) + while os.path.exists(filename): + filename = filename[: - n] + "{}.{}".format(np.random.randint(10), + extension) + return filename diff --git a/rrompy/utilities/number_theory.py b/rrompy/utilities/number_theory.py new file mode 100644 index 0000000..b54ef97 --- /dev/null +++ b/rrompy/utilities/number_theory.py @@ -0,0 +1,71 @@ +#!/usr/bin/python +# 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 + +__all__ = ['squareResonances','getLowestPrimeFactor','primeFactorize'] + +prime_v = [] + +def squareResonances(a:int, b:int, zero_terms : bool = True): + spectrum = [] + a = max(int(np.floor(a)), 0) + b = max(int(np.ceil(b)), 0) + global prime_v + if len(prime_v) == 0: + prime_v = [2, 3] + if a > prime_v[-1]: + for i in range(prime_v[-1], a, 2): + getLowestPrimeFactor(i) + for i in range(a, b + 1): + spectrum = spectrum + [i] * countSquareSums(i, zero_terms) + return np.array(spectrum) + +def getLowestPrimeFactor(n:int): + global prime_v + for x in prime_v: + if n % x == 0: + return x + if prime_v[-1] < n: + prime_v = prime_v + [n] + return n + +def primeFactorize(n:int): + factors = [] + number = n + while number > 1: + factor = getLowestPrimeFactor(number) + factors.append(factor) + number = int(number / factor) + if n < -1: + factors[0] = -factors[0] + return list(factors) + +def countSquareSums(n:int, zero_terms : bool = True): + if n < 2: return (n + 1) * zero_terms + factors = primeFactorize(n) + funique, fcounts = np.unique(factors, return_counts = True) + count = 1 + for fac, con in zip(funique, fcounts): #using number theory magic + if fac % 4 == 1: + count = count * (con + 1) + elif fac % 4 == 3 and con % 2 == 1: + return 0 + return count + (2 * zero_terms - 1) * (int(n ** .5) ** 2 == n) + diff --git a/rrompy/utilities/purge_dict.py b/rrompy/utilities/purge_dict.py new file mode 100644 index 0000000..e480405 --- /dev/null +++ b/rrompy/utilities/purge_dict.py @@ -0,0 +1,41 @@ +#!/usr/bin/python +# 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 +from rrompy.utilities.find_dict_str_key import findDictStrKey +from rrompy.utilities.types import ListAny, DictAny + +__all__ = ['purgeDict'] + +def purgeDict(dct:DictAny, allowedKeys : ListAny = [], silent : bool = False, + complement : bool = False, dictname : str = ""): + if dictname != "": + dictname = " in " + dictname + dctcp = {} + for key in dct.keys(): + akey = findDictStrKey(key, allowedKeys) + if (akey is None) != complement: + if not silent: + warnings.warn("Ignoring key {0}{2} with value {1}."\ + .format(key, dct[key], dictname), stacklevel = 2) + else: + if akey is None: + akey = key + dctcp[akey] = dct[key] + return dctcp diff --git a/rrompy/utilities/purge_list.py b/rrompy/utilities/purge_list.py new file mode 100644 index 0000000..d92613a --- /dev/null +++ b/rrompy/utilities/purge_list.py @@ -0,0 +1,41 @@ +#!/usr/bin/python +# 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 +from rrompy.utilities.find_dict_str_key import findDictStrKey +from rrompy.utilities.types import ListAny + +__all__ = ['purgeList'] + +def purgeList(lst:ListAny, allowedEntries : ListAny = [], + silent : bool = False, complement : bool = False, + listname : str = ""): + if listname != "": + listname = " in " + listname + lstcp = [] + for x in lst: + ax = findDictStrKey(x, allowedEntries) + if (ax is None) != complement: + if not silent: + warnings.warn("Ignoring entry {0}{1}.".format(x, listname), + stacklevel = 2) + else: + lstcp = lstcp + [ax] + return lstcp + diff --git a/rrompy/utilities/sobol.py b/rrompy/utilities/sobol.py new file mode 100644 index 0000000..01c8c3a --- /dev/null +++ b/rrompy/utilities/sobol.py @@ -0,0 +1,80 @@ +#!/usr/bin/python +# 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 + +__all__ = ['sobolGenerate'] + +def low2bit(n): + if n < 0: + raise Exception("Only positive integers allowed.") + j = 2 + while (n % j) // (j / 2) == 1: j = j * 2 + return int(np.log2(j)) - 1 + +def sobolGenerate(dim_num, n, seed = 0): + SEEDBITS = 30 + v = np.zeros((dim_num, SEEDBITS), dtype = int) + v[:, 0] = 1 + v[2:, 1] = np.array([1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3,3,1,3,1,3, + 1,3,1,1,3,1,3,1,3,1,3])[:max(0,dim_num-2)] + v[3:, 2] = np.array([7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3,1,7,5,1,3,3,7, + 5,1,1,5,7,7,5,1,3,3])[:max(0,dim_num-3)] + v[5:, 3] = np.array([1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3,15,7,9,13,9,1, + 11,7,5,15,1,15,11,5,3,1,7,9])[:max(0,dim_num-5)] + v[7:, 4] = np.array([9,3,27,15,29,21,23,19,11,25,7,13,17,1,25,29,3,31,11,5, + 23,27,19,21,5,1,17,13,7,15,9,31,9])[:max(0,dim_num-7)] + v[13:, 5] = np.array([37,33,7,5,11,39,63,27,17,15,23,29,3,21,13,31,25,9,49, + 33,19,29,11,19,27,15,25])[:max(0,dim_num-13)] + v[19:, 6] = np.array([13,33,115,41,79,17,29,119,75,73,105,7,59,65,21,3,113, + 61,89,45,107])[:max(0,dim_num-19)] + v[37:, 7] = np.array([7,23,39])[:max(0,dim_num-37)] + poly = [1,3,7,11,13,19,25,37,59,47,61,55,41,67,97,91,109,103,115,131,193, + 137,145,143,241,157,185,167,229,171,213,191,253,203,211,239,247, + 285,369,299] + + v[0, :] = 1 + for i in range(1, dim_num): + j = poly[i] + m = int(1 + np.floor(np.log2(j))) + includ = np.array([int(x) for x in "{0:b}".format(j)], dtype = int) + for j in range(m - 1, SEEDBITS): + newv = v[i, j - m + 1] + l = 1 + for k in range(1, m): + l *= 2 + if includ[k]: + newv = np.bitwise_xor(newv, l * v[i, j - k]) + v[i, j] = newv + v = np.multiply(v, np.power(2, np.arange(SEEDBITS)[::-1])) + recipd = np.power(2., - SEEDBITS) + lastq = np.zeros(dim_num, dtype = int) + #### + if seed < 0: + seed = 0 + for seed_temp in range(seed): + l = low2bit(seed_temp) + lastq = np.bitwise_xor(lastq, v[:, l]) + r = np.empty((n, dim_num)) + for seed_temp in range(seed, seed + n): + l = low2bit(seed_temp) + r[seed_temp - seed, :] = lastq * recipd + lastq = np.bitwise_xor(lastq, v[:, l]) + return r + diff --git a/rrompy/utilities/types.py b/rrompy/utilities/types.py new file mode 100644 index 0000000..96524b4 --- /dev/null +++ b/rrompy/utilities/types.py @@ -0,0 +1,50 @@ +#!/usr/bin/python +# 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 typing import TypeVar, List, Tuple, Dict, Any + +__all__ = ['TupleAny','ListAny','DictAny','Np1D','Np2D','Np1DLst','N2FSExpr', + 'FenSpace','FenExpr','FenFunc','HFEng','HSEng','ROMEng','GenExpr', + 'strLst','BfSExpr'] + +# ANY +TupleAny = Tuple[Any] +ListAny = List[Any] +DictAny = Dict[Any, Any] + +# NUMPY +Np1D = TypeVar("NumPy 1D array") +Np2D = TypeVar("NumPy 2D array") +Np1DLst = TypeVar("NumPy 1D array or list of NumPy 1D array") +N2FSExpr = TypeVar("NumPy 2D array, float or str") + +# FENICS +FenSpace = TypeVar("FEniCS FESpace") +FenExpr = TypeVar("FEniCS expression") +FenFunc = TypeVar("FEniCS function") + +# ROM +HFEng = TypeVar("High fidelity engine") +HSEng = TypeVar("Hilbert space engine") +ROMEng = TypeVar("ROM engine") + +# OTHERS +GenExpr = TypeVar("Generic expression") +strLst = TypeVar("str or list of str") +BfSExpr = TypeVar("Boolean function or string") diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..ecef72c --- /dev/null +++ b/setup.py @@ -0,0 +1,47 @@ +#!/usr/lib/python +# Copyright (C) 2015-2018 by the RBniCS authors +# +# This file is part of RBniCS. +# +# RBniCS 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. +# +# RBniCS 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 RBniCS. If not, see . +# + +import os +from setuptools import find_packages, setup + +rrompy_directory = os.path.abspath(os.path.dirname(os.path.realpath(__file__))) +#rrompy_directory = os.path.join(rrompy_directory, 'rrompy') + +setup(name="RROMPy", + description="Rational reduced order modelling in Python", + long_description="Rational reduced order modelling in Python", + author="Davide Pradovera", + author_email="davide.pradovera@epfl.ch", + version="0.0.dev1", + license="GNU Library or Lesser General Public License (LGPL)", + classifiers=[ + "Development Status :: 1 - Planning" + "Intended Audience :: Developers", + "Intended Audience :: Science/Research", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.4", + "Programming Language :: Python :: 3.5", + "Programming Language :: Python :: 3.6", + "License :: OSI Approved :: GNU Library or Lesser General Public License (LGPL)", + "Topic :: Scientific/Engineering :: Mathematics", + "Topic :: Software Development :: Libraries :: Python Modules", + ], + packages=find_packages(rrompy_directory), + zip_safe=False + )