diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py index 4317de6..aa3c9d8 100644 --- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py +++ b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py @@ -1,185 +1,162 @@ # 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 numpy.polynomial.polynomial import polyfit as fit import fenics as fen from rrompy.utilities.base.types import (ScOp, List, Tuple, paramVal, Np1D, FenExpr) -from rrompy.solver.fenics import fenZERO, fenONE +from rrompy.solver.fenics import fenZERO from rrompy.hfengines.linear_problem.helmholtz_problem_engine import ( HelmholtzProblemEngine) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['HelmholtzSquareDomainProblemEngine'] class HelmholtzSquareDomainProblemEngine(HelmholtzProblemEngine): """ Solver for square Helmholtz problems with parametric laplacian. - \dxx u - mu_2**2 \dyy u - mu_1**2 * u = f(mu_2) in \Omega = [0,\pi]^2 u = 0 on \partial\Omega """ def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [12. ** .5, 1.], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.nAs, self.nbs = 6, 11 * 12 // 2 self.npar = 2 self.rescalingExp = [2., 1.] + self._theta = theta + self._kappa = kappa pi = np.pi mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi), 3 * n, 3 * n) self.V = fen.FunctionSpace(mesh, "P", 1) - - c, s = np.cos(theta), np.sin(theta) - x, y = fen.SpatialCoordinate(mesh)[:] - self.forcingTerm = [fen.cos(kappa * (c * x + s / self.mu0(0, 1) * y)), - fen.sin(kappa * (c * x + s / self.mu0(0, 1) * y))] - self.forcingTermDer = kappa * s * y - def getExtraFactorB(self, mu : paramVal = [], - derI : int = 0) -> Tuple[FenExpr, FenExpr]: - """Compute extra expression in RHS.""" + def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: + """Compute forcing term.""" mu = self.checkParameter(mu) - def getPowMinusj(x, power): - powR = x ** power - powI = fenZERO - if power % 2 == 1: - powR, powI = powI, powR - if power % 4 > 1: - powR, powI = - powR, - powI - return powR, powI - vbMng(self, "INIT", - "Assembling auxiliary expression for forcing term derivative.", - 25) - if derI == 0: return fenONE, fenZERO - coeffs = np.zeros(derI + 1) - coeffs[1] = - 2. - for j in range(2, derI + 1): - coeffs[1 :] = (-2. / j * coeffs[1 :] - - (3 - (1 + 2 * np.arange(derI)) / j) * coeffs[: -1]) - for j in range(derI): - powR, powI = getPowMinusj(self.forcingTermDer, derI - j) - mupBase = coeffs[j + 1] * mu(0, 1) ** (- 3 * derI + 2 * j) - mupR, mupI = np.real(mupBase), np.imag(mupBase) - if j == 0: - exprR = mupR * powR - mupI * powI - exprI = mupI * powR + mupR * powI - else: - exprR += mupR * powR - mupI * powI - exprI += mupI * powR + mupR * powI - vbMng(self, "DEL", "Done assembling auxiliary expression.", 25) - return exprR, exprI - + vbMng(self, "INIT", ("Assembling base expression for forcing term " + "at {}.").format(mu), 25) + mu = mu(0, 1) + c, s = np.cos(self._theta), np.sin(self._theta) + x, y = fen.SpatialCoordinate(self.V.mesh())[:] + forcingTerm = [fen.cos(self._kappa * (c * x + s / mu * y)), + fen.sin(self._kappa * (c * x + s / mu * y))] + vbMng(self, "DEL", "Done assembling base expression.", 25) + return forcingTerm + def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) self.autoSetDS() for j in range(2, 5): if derI <= j and self.As[j] is None: self.As[j] = self.checkAInBounds(-1) if derI <= 0 and self.As[0] is None: vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 1 and self.As[1] is None: vbMng(self, "INIT", "Assembling operator term A1.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) nRe, nIm = self.refractionIndex n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm parsRe = self.iterReduceQuadratureDegree(zip([n2Re], ["refractionIndexSquaredReal"])) parsIm = self.iterReduceQuadratureDegree(zip([n2Im], ["refractionIndexSquaredImag"])) a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0) + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 5 and self.As[5] is None: vbMng(self, "INIT", "Assembling operator term A5.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a5Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx self.As[5] = fenics2Sparse(a5Re, {}, DirichletBC0, 0) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" - raise Exception("Deprecated.") mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) + bDEIMCoeffs = None for j in range(derI, nbsTot): derH = hashI(j, self.npar) if bs[j] is None: if np.sum(derH) != derH[-1]: if homogeneized: self.bsH[j] = self.checkbInBounds(-1) + Ader = self.A(0, derH) + self.bsH[j] -= Ader.dot(self.liftedDirichletDatum) else: self.bs[j] = self.checkbInBounds(-1) continue vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) - if j == 0: - u0Re, u0Im = self.DirichletDatum - else: - u0Re, u0Im = fenZERO, fenZERO - if j < self.nbs: - fRe, fIm = self.forcingTerm - cRe, cIm = self.getExtraFactorB(self.mu0, derH[-1]) - cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe - else: - cfRe, cfIm = fenZERO, fenZERO - parsRe = self.iterReduceQuadratureDegree(zip([cfRe], - ["forcingTermDer{}Real".format(j)])) - parsIm = self.iterReduceQuadratureDegree(zip([cfIm], - ["forcingTermDer{}Imag".format(j)])) - L0Re = fen.dot(cfRe, self.v) * fen.dx - L0Im = fen.dot(cfIm, self.v) * fen.dx - DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary) - DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary) - b = (fenics2Vector(L0Re, parsRe, DBCR, 1) - + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1)) + if bDEIMCoeffs is None: + bDEIM = np.empty((self.nbs, self.spacedim()), + dtype = np.complex) + muDEIM = np.linspace(.5, 4., + np.sum(hashI(nbsTot - 1, self.npar))) + for jj, muD in enumerate(muDEIM): + fRe, fIm = self.getForcingTerm(muD) + parsRe = self.iterReduceQuadratureDegree(zip([fRe], + ["forcingTerm{}Real".format(jj)])) + parsIm = self.iterReduceQuadratureDegree(zip([fIm], + ["forcingTerm{}Imag".format(jj)])) + LR = fen.dot(fRe, self.v) * fen.dx + LI = fen.dot(fIm, self.v) * fen.dx + DBC0 = fen.DirichletBC(self.V, fenZERO, + self.DirichletBoundary) + bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1) + + 1.j * fenics2Vector(LI, parsIm, DBC0, 1)) + bDEIMCoeffs = fit(muDEIM, bDEIM, len(muDEIM) - 1) + b = bDEIMCoeffs[derH[-1]] if homogeneized: - Ader = self.A(self.mu0, derH) + Ader = self.A(0, derH) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) - return self._assembleb(mu, der, derI, homogeneized, self.mu0) + return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py b/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py index ff0c864..73cd317 100644 --- a/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py +++ b/rrompy/hfengines/linear_problem/bidimensional/laplace_disk_gaussian_2.py @@ -1,115 +1,115 @@ # 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.utilities.base.types import Np1D, Tuple, List, FenExpr, paramVal from rrompy.hfengines.linear_problem.laplace_disk_gaussian import ( LaplaceDiskGaussian) from rrompy.solver.fenics import fenZERO, fenONE from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyException from rrompy.solver.fenics import fenics2Vector __all__ = ['LaplaceDiskGaussian2'] class LaplaceDiskGaussian2(LaplaceDiskGaussian): """ Solver for disk Laplace problems with parametric forcing term center. - \Delta u = C exp(-.5 * ||\cdot - (mu1, mu2)||^2) in \Omega = B(0, 5) u = 0 on \partial\Omega. """ def __init__(self, n:int, mu0 : paramVal = [0., 0.], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(n = n, mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.nbs = 1 self.npar = 2 def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" mu = self.checkParameter(mu) if mu != self.forcingTermMu: vbMng(self, "INIT", "Assembling base expression for forcing term.", 25) x, y = fen.SpatialCoordinate(self.V.mesh())[:] C = np.exp(-.5 * (mu(0, 0) ** 2. + mu(0, 1) ** 2.)) CR, CI = np.real(C), np.imag(C) f0 = (2 * np.pi) ** -.5 * fen.exp(-.5 * (x ** 2. + y ** 2.)) muxR, muxI = np.real(mu(0, 0)), np.imag(mu(0, 0)) muyR, muyI = np.real(mu(0, 1)), np.imag(mu(0, 1)) f1R = fen.exp(muxR * x + muyR * y) * fen.cos(muxI * x + muyI * y) f1I = fen.exp(muxR * x + muyR * y) * fen.sin(muxI * x + muyI * y) self.forcingTerm = [f0 * (CR * f1R - CI * f1I), f0 * (CR * f1I + CI * f1R)] self.forcingTermMu = mu vbMng(self, "DEL", "Done assembling base expression.", 25) return self.forcingTerm def computebsFactors(self): pass def getExtraFactorB(self, mu : paramVal = [], derI : int = 0) -> Tuple[FenExpr, FenExpr]: if derI == 0: return [fenONE, fenZERO] raise RROMPyException("Not implemented.") def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) for j in range(derI, nbsTot): if bs[j] is None: vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) if j < self.nbs: fRe, fIm = self.getForcingTerm(mu) cRe, cIm = self.getExtraFactorB(mu, j) cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe else: cfRe, cfIm = fenZERO, fenZERO parsRe = self.iterReduceQuadratureDegree(zip([cfRe], ["forcingTermDer{}Real".format(j)])) parsIm = self.iterReduceQuadratureDegree(zip([cfIm], ["forcingTermDer{}Imag".format(j)])) L0Re = fen.dot(cfRe, self.v) * fen.dx L0Im = fen.dot(cfIm, self.v) * fen.dx DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1) + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1)) if homogeneized: - Ader = self.A(self.mu0, hashI(j, self.npar)) + Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py index 716069e..a50171c 100644 --- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py @@ -1,196 +1,162 @@ # 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 +from numpy.polynomial.polynomial import polyfit as fit import fenics as fen from rrompy.utilities.base.types import (Np1D, ScOp, Tuple, List, FenExpr, paramVal) from rrompy.solver.fenics import fenZERO from .helmholtz_problem_engine import HelmholtzProblemEngine from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['HelmholtzSquareBubbleDomainProblemEngine'] class HelmholtzSquareBubbleDomainProblemEngine(HelmholtzProblemEngine): """ Solver for square bubble Helmholtz problems with parametric domain heigth. - \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi] u = 0 on \Gamma_mu = \partial\Omega_mu with exact solution square bubble times plane wave. """ def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [1.], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) - self.nAs, self.nbs = 3, 20 + self.nAs, self.nbs = 3, 15 self.kappa = kappa self.theta = theta self.forcingTermMu = np.nan mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(np.pi, np.pi), 3 * n, 3 * n) self.V = fen.FunctionSpace(mesh, "P", 1) self.rescalingExp = [1.] def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" mu = self.checkParameter(mu) - if mu != self.forcingTermMu: - vbMng(self, "INIT", "Assembling base expression for forcing term.", - 25) - pi = np.pi - c, s = np.cos(self.theta), np.sin(self.theta) - x, y = fen.SpatialCoordinate(self.V.mesh())[:] - muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0)) - mu2R, mu2I = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.) - C = 16. / pi ** 4. - bR = C * (2 * (x * (pi - x) + y * (pi - y)) - + (self.kappa * s) ** 2. * (mu2R - 1.) - * x * (pi - x) * y * (pi - y)) - bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y * (pi - y) - + s * x * (pi - x) * (pi - 2 * y)) - + (self.kappa * s) ** 2. * mu2I - * x * (pi - x) * y * (pi - y)) - wR = (fen.cos(self.kappa * (c * x + s * muR * y)) - * fen.exp(self.kappa * s * muI * y)) - wI = (fen.sin(self.kappa * (c * x + s * muR * y)) - * fen.exp(self.kappa * s * muI * y)) - self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI] - self.forcingTermMu = mu - vbMng(self, "DEL", "Done assembling base expression.", 25) - return self.forcingTerm - - def getExtraFactorB(self, mu : paramVal = [], - derI : int = 0) -> Tuple[FenExpr, FenExpr]: - """Compute extra expression in RHS.""" - mu = self.checkParameter(mu) - def getPowMinusj(x, power): - powR = x ** power - powI = fenZERO - if power % 2 == 1: - powR, powI = powI, powR - if (power + 3) % 4 < 2: - powR, powI = - powR, - powI - return powR, powI - vbMng(self, "INIT", - "Assembling auxiliary expression for forcing term derivative.", - 25) - from scipy.special import factorial as fact - y = fen.SpatialCoordinate(self.V.mesh())[1] - powR, powI = [(self.kappa * np.sin(self.theta)) ** derI * k\ - for k in getPowMinusj(y, derI)] + vbMng(self, "INIT", ("Assembling base expression for forcing term " + "at {}.").format(mu), 25) + pi = np.pi + c, s = np.cos(self.theta), np.sin(self.theta) + x, y = fen.SpatialCoordinate(self.V.mesh())[:] + muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0)) mu2R, mu2I = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.) - exprR = mu2R * powR - mu2I * powI - exprI = mu2I * powR + mu2R * powI - if derI >= 1: - muR, muI = np.real(2. * mu(0, 0)), np.imag(2. * mu(0, 0)) - powR, powI = [(self.kappa * np.sin(self.theta)) ** (derI - 1) * k\ - * derI for k in getPowMinusj(y, derI - 1)] - exprR += muR * powR - muI * powI - exprI += muI * powR + muR * powI - if derI >= 2: - powR, powI = [(self.kappa * np.sin(self.theta)) ** (derI - 2) * k\ - * derI * (derI - 1) for k in getPowMinusj(y, derI - 2)] - exprR += powR - exprI += powI - fac = fact(derI) - vbMng(self, "DEL", "Done assembling auxiliary expression.", 25) - return [exprR / fac, exprI / fac] - + C = 16. / pi ** 4. + bR = C * (2 * (x * (pi - x) + y * (pi - y)) + + (self.kappa * s) ** 2. * (mu2R - 1.) + * x * (pi - x) * y * (pi - y)) + bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y * (pi - y) + + s * x * (pi - x) * (pi - 2 * y)) + + (self.kappa * s) ** 2. * mu2I + * x * (pi - x) * y * (pi - y)) + wR = (fen.cos(self.kappa * (c * x + s * muR * y)) + * fen.exp(self.kappa * s * muI * y)) + wI = (fen.sin(self.kappa * (c * x + s * muR * y)) + * fen.exp(self.kappa * s * muI * y)) + fRe, fIm = bR * wR + bI * wI, bI * wR - bR * wI + cRe, cIm = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.) + forcingTerm = [cRe * fRe - cIm * fIm + fenZERO, + cRe * fIm + cIm * fRe + fenZERO] + vbMng(self, "DEL", "Done assembling base expression.", 25) + return forcingTerm + def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) self.autoSetDS() if derI <= 0 and self.As[0] is None: vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 1 and self.As[1] is None: self.As[1] = self.checkAInBounds(-1) if derI <= 2 and self.As[2] is None: vbMng(self, "INIT", "Assembling operator term A2.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) nRe, nIm = self.refractionIndex n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2) k2n2Re = k2Re * n2Re - k2Im * n2Im k2n2Im = k2Re * n2Im + k2Im * n2Re parsRe = self.iterReduceQuadratureDegree(zip([k2n2Re], ["kappaSquaredRefractionIndexSquaredReal"])) parsIm = self.iterReduceQuadratureDegree(zip([k2n2Im], ["kappaSquaredRefractionIndexSquaredImag"])) a2Re = (fen.dot(self.u.dx(0), self.v.dx(0)) - k2n2Re * fen.dot(self.u, self.v)) * fen.dx a2Im = - k2n2Im * fen.dot(self.u, self.v) * fen.dx self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0) + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" - raise Exception("Deprecated.") mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) + bDEIMCoeffs = None for j in range(derI, nbsTot): if bs[j] is None: vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) - if j < self.nbs: - fRe, fIm = self.getForcingTerm(self.mu0) - cRe, cIm = self.getExtraFactorB(self.mu0, j) - cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe - else: - cfRe, cfIm = fenZERO, fenZERO - parsRe = self.iterReduceQuadratureDegree(zip([cfRe], - ["forcingTermDer{}Real".format(j)])) - parsIm = self.iterReduceQuadratureDegree(zip([cfIm], - ["forcingTermDer{}Imag".format(j)])) - L0Re = fen.dot(cfRe, self.v) * fen.dx - L0Im = fen.dot(cfIm, self.v) * fen.dx - DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + if bDEIMCoeffs is None: + bDEIM = np.empty((self.nbs, self.spacedim()), + dtype = np.complex) + muDEIM = np.linspace(.5, 4., self.nbs) + for jj, muD in enumerate(muDEIM): + fRe, fIm = self.getForcingTerm(muD) + parsRe = self.iterReduceQuadratureDegree(zip([fRe], + ["forcingTerm{}Real".format(jj)])) + parsIm = self.iterReduceQuadratureDegree(zip([fIm], + ["forcingTerm{}Imag".format(jj)])) + LR = fen.dot(fRe, self.v) * fen.dx + LI = fen.dot(fIm, self.v) * fen.dx + DBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) - b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1) - + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1)) + bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1) + + 1.j * fenics2Vector(LI, parsIm, DBC0, 1)) + bDEIMCoeffs = fit(muDEIM, bDEIM, self.nbs - 1) + b = bDEIMCoeffs[j] if homogeneized: - Ader = self.A(self.mu0, hashI(j, self.npar)) + Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) - return self._assembleb(mu, der, derI, homogeneized, self.mu0) + return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py index a5e4cfb..8899593 100644 --- a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py +++ b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py @@ -1,337 +1,337 @@ # 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.base.problem_engine_base import ProblemEngineBase from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal from rrompy.solver.fenics import (fenZERO, fenONE, L2InverseNormMatrix, H1NormMatrix, Hminus1NormMatrix) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.parameter import checkParameter from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['LaplaceBaseProblemEngine'] class LaplaceBaseProblemEngine(ProblemEngineBase): """ Solver for generic Laplace problems. - \nabla \cdot (a \nabla u) = f in \Omega u = u0 on \Gamma_D \partial_nu = g1 on \Gamma_N \partial_nu + h u = g2 on \Gamma_R Attributes: verbosity: Verbosity level. BCManager: Boundary condition manager. V: Real FE space. u: Generic trial functions for variational form evaluation. v: Generic test functions for variational form evaluation. As: Scipy sparse array representation (in CSC format) of As. bs: Numpy array representation of bs. bsH: Numpy array representation of homogeneized bs. energyNormMatrix: Scipy sparse matrix representing inner product over V. energyNormDualMatrix: Scipy sparse matrix representing inner product over V'. dualityMatrix: Scipy sparse matrix representing duality V-V'. energyNormPartialDualMatrix: Scipy sparse matrix representing dual inner product between Riesz representers V-V. liftedDirichletDatum: Dofs of Dirichlet datum lifting. mu0BC: Mu value of last Dirichlet datum lifting. degree_threshold: Threshold for ufl expression interpolation degree. omega: Value of omega. diffusivity: Value of a. forcingTerm: Value of f. DirichletDatum: Value of u0. NeumannDatum: Value of g1. RobinDatumG: Value of g2. RobinDatumH: Value of h. DirichletBoundary: Function handle to \Gamma_D. NeumannBoundary: Function handle to \Gamma_N. RobinBoundary: Function handle to \Gamma_R. ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). A0: Scipy sparse array representation (in CSC format) of A0. b0: Numpy array representation of b0. dsToBeSet: Whether ds needs to be set. """ _energyDualNormCompress = None def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.mu0 = checkParameter(mu0) self.npar = self.mu0.shape[1] self.omega = np.abs(self.mu0(0, 0)) if self.npar > 0 else 0. self.diffusivity = fenONE self.forcingTerm = fenZERO self.DirichletDatum = fenZERO self.NeumannDatum = fenZERO self.RobinDatumG = fenZERO self.RobinDatumH = fenZERO @property def V(self): """Value of V.""" return self._V @V.setter def V(self, V): ProblemEngineBase.V.fset(self, V) self.dsToBeSet = True @property def diffusivity(self): """Value of a.""" return self._diffusivity @diffusivity.setter def diffusivity(self, diffusivity): self.resetAs() if not isinstance(diffusivity, (list, tuple,)): diffusivity = [diffusivity, fenZERO] self._diffusivity = diffusivity @property def forcingTerm(self): """Value of f.""" return self._forcingTerm @forcingTerm.setter def forcingTerm(self, forcingTerm): self.resetbs() if not isinstance(forcingTerm, (list, tuple,)): forcingTerm = [forcingTerm, fenZERO] self._forcingTerm = forcingTerm @property def DirichletDatum(self): """Value of u0.""" return self._DirichletDatum @DirichletDatum.setter def DirichletDatum(self, DirichletDatum): self.resetbs() 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): self.resetbs() 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): self.resetbs() if not isinstance(RobinDatumG, (list, tuple,)): RobinDatumG = [RobinDatumG, fenZERO] self._RobinDatumG = RobinDatumG @property def RobinDatumH(self): """Value of h.""" return self._RobinDatumH @RobinDatumH.setter def RobinDatumH(self, RobinDatumH): self.resetAs() if not isinstance(RobinDatumH, (list, tuple,)): RobinDatumH = [RobinDatumH, fenZERO] self._RobinDatumH = RobinDatumH @property def DirichletBoundary(self): """Function handle to DirichletBoundary.""" return self.BCManager.DirichletBoundary @DirichletBoundary.setter def DirichletBoundary(self, DirichletBoundary): self.resetAs() self.resetbs() self.BCManager.DirichletBoundary = DirichletBoundary @property def NeumannBoundary(self): """Function handle to NeumannBoundary.""" return self.BCManager.NeumannBoundary @NeumannBoundary.setter def NeumannBoundary(self, NeumannBoundary): self.resetAs() self.resetbs() self.dsToBeSet = True self.BCManager.NeumannBoundary = NeumannBoundary @property def RobinBoundary(self): """Function handle to RobinBoundary.""" return self.BCManager.RobinBoundary @RobinBoundary.setter def RobinBoundary(self, RobinBoundary): self.resetAs() self.resetbs() self.dsToBeSet = True self.BCManager.RobinBoundary = RobinBoundary def autoSetDS(self): """Set FEniCS boundary measure based on boundary function handles.""" if self.dsToBeSet: vbMng(self, "INIT", "Initializing boundary measures.", 20) mesh = self.V.mesh() NB = self.NeumannBoundary RB = self.RobinBoundary boundary_markers = fen.MeshFunction("size_t", mesh, mesh.topology().dim() - 1) NB.mark(boundary_markers, 0) RB.mark(boundary_markers, 1) self.ds = fen.Measure("ds", domain = mesh, subdomain_data = boundary_markers) self.dsToBeSet = False vbMng(self, "DEL", "Done assembling boundary measures.", 20) def buildEnergyNormForm(self): """ Build sparse matrix (in CSR format) representative of scalar product. """ vbMng(self, "INIT", "Assembling energy matrix.", 20) self.energyNormMatrix = H1NormMatrix(self.V, np.abs(self.omega)**2) vbMng(self, "DEL", "Done assembling energy matrix.", 20) def buildEnergyNormDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product. """ vbMng(self, "INIT", "Assembling energy dual matrix.", 20) self.energyNormDualMatrix = Hminus1NormMatrix( self.V, np.abs(self.omega)**2, compressRank = self._energyDualNormCompress) vbMng(self, "DEL", "Done assembling energy dual matrix.", 20) def buildDualityPairingForm(self): """Build sparse matrix (in CSR format) representative of duality.""" vbMng(self, "INIT", "Assembling duality matrix.", 20) self.dualityMatrix = L2InverseNormMatrix( self.V, solverType = self._solver, solverArgs = self._solverArgs, compressRank = self._dualityCompress) vbMng(self, "DEL", "Done assembling duality matrix.", 20) def buildEnergyNormPartialDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product without duality. """ vbMng(self, "INIT", "Assembling energy partial dual matrix.", 20) self.energyNormPartialDualMatrix = Hminus1NormMatrix( self.V, np.abs(self.omega)**2, compressRank = self._energyDualNormCompress, duality = False) vbMng(self, "DEL", "Done assembling energy partial dual matrix.", 20) def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) if derI <= 0 and self.As[0] is None: self.autoSetDS() vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) aRe, aIm = self.diffusivity hRe, hIm = self.RobinDatumH termNames = ["diffusivity", "RobinDatumH"] parsRe = self.iterReduceQuadratureDegree(zip([aRe, hRe], [x + "Real" for x in termNames])) parsIm = self.iterReduceQuadratureDegree(zip([aIm, hIm], [x + "Imag" for x in termNames])) 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)) self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1) + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) for j in range(derI, nbsTot): if bs[j] is None: self.autoSetDS() vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) termNames, terms = [], [] if j == 0: u0Re, u0Im = self.DirichletDatum fRe, fIm = self.forcingTerm g1Re, g1Im = self.NeumannDatum g2Re, g2Im = self.RobinDatumG termNames += ["forcingTerm", "NeumannDatum", "RobinDatumG"] terms += [[fRe, fIm], [g1Re, g1Im], [g2Re, g2Im]] else: u0Re, u0Im = fenZERO, fenZERO fRe, fIm = fenZERO, fenZERO g1Re, g1Im = fenZERO, fenZERO g2Re, g2Im = fenZERO, fenZERO if len(termNames) > 0: parsRe = self.iterReduceQuadratureDegree(zip( [term[0] for term in terms], [x + "Real" for x in termNames])) parsIm = self.iterReduceQuadratureDegree(zip( [term[1] for term in terms], [x + "Imag" for x in termNames])) else: parsRe, parsIm = {}, {} 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)) DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary) b = (fenics2Vector(L0Re, parsRe, DBCR, 1) + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1)) if homogeneized: - Ader = self.A(self.mu0, hashI(j, self.npar)) + Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py index 30dab29..5980548 100644 --- a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py +++ b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py @@ -1,148 +1,112 @@ # 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 numpy.polynomial.polynomial import polyfit as fit import fenics as fen from rrompy.utilities.base.types import Np1D, Tuple, FenExpr, paramVal from .laplace_base_problem_engine import LaplaceBaseProblemEngine from rrompy.solver.fenics import fenZERO, fenONE from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.solver.fenics import fenics2Vector __all__ = ['LaplaceDiskGaussian'] class LaplaceDiskGaussian(LaplaceBaseProblemEngine): """ Solver for disk Laplace problems with parametric forcing term center. - \Delta u = C exp(-.5 * ||\cdot - (mu, 0)||^2) in \Omega = B(0, 5) u = 0 on \partial\Omega. """ - def __init__(self, n:int, mu0 : paramVal = [0.], degree_threshold : int = np.inf, - verbosity : int = 10, timestamp : bool = True): + def __init__(self, n:int, mu0 : paramVal = [0.], + degree_threshold : int = np.inf, verbosity : int = 10, + timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) - self.nbs = 20 - self.computebsFactors() + self.nbs = 19 self.forcingTermMu = np.nan import mshr mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), 3 * n) self.V = fen.FunctionSpace(mesh, "P", 1) def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" mu = self.checkParameter(mu) - if mu != self.forcingTermMu: - vbMng(self, "INIT", "Assembling base expression for forcing term.", - 25) - x, y = fen.SpatialCoordinate(self.V.mesh())[:] - C = np.exp(-.5 * mu(0, 0) ** 2.) - CR, CI = np.real(C), np.imag(C) - f0 = (2 * np.pi) ** -.5 * fen.exp(-.5 * (x ** 2. + y ** 2.)) - muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0)) - f1R = fen.exp(muR * x) * fen.cos(muI * x) - f1I = fen.exp(muR * x) * fen.sin(muI * x) - self.forcingTerm = [f0 * (CR * f1R - CI * f1I), - f0 * (CR * f1I + CI * f1R)] - self.forcingTermMu = mu - vbMng(self, "DEL", "Done assembling base expression.", 25) - return self.forcingTerm - - def computebsFactors(self): - self.bsFactors = np.zeros((self.nbs, self.nbs), dtype = float) - self.bsFactors[0, 0] = 1. - self.bsFactors[1, 1] = 1. - for j in range(2, self.nbs): - l = (j + 1) % 2 + 1 - J = np.arange(l, j + 1, 2) - self.bsFactors[j, J] = self.bsFactors[j - 1, J - 1] - if l == 2: - l = 0 - J = np.arange(l, j, 2) - self.bsFactors[j, J] += np.multiply(- 1 - J, - self.bsFactors[j - 1, J + 1]) - self.bsFactors[j, l : j + 2 : 2] /= j - - def getExtraFactorB(self, mu : paramVal = [], - derI : int = 0) -> Tuple[FenExpr, FenExpr]: - """Compute extra expression in RHS.""" - mu = self.checkParameter(mu) - vbMng(self, "INIT", - "Assembling auxiliary expression for forcing term derivative.", - 25) + vbMng(self, "INIT", ("Assembling base expression for forcing term " + "at {}.").format(mu), 25) + x, y = fen.SpatialCoordinate(self.V.mesh())[:] + C = np.exp(-.5 * mu(0, 0) ** 2.) + CR, CI = np.real(C), np.imag(C) + f0 = (2 * np.pi) ** -.5 * fen.exp(-.5 * (x ** 2. + y ** 2.)) muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0)) - x = fen.SpatialCoordinate(self.V.mesh())[0] - l = derI % 2 - if l == 0: - powR, powI = fenONE, fenZERO - else: - powR, powI = x - muR, fen.Constant(muI) - exprR, exprI = [self.bsFactors[derI, l] * k for k in [powR, powI]] - for j in range(l + 2, derI + 1, 2): - for _ in range(2): - powR, powI = (powR * (x - muR) - powI * muI, - powR * muI + powI * (x - muR)) - exprR += self.bsFactors[derI, j] * powR - exprI += self.bsFactors[derI, j] * powI - vbMng(self, "DEL", "Done assembling auxiliary expression.", 25) - return[exprR, exprI] - + f1R = fen.exp(muR * x) * fen.cos(muI * x) + f1I = fen.exp(muR * x) * fen.sin(muI * x) + forcingTerm = [f0 * (CR * f1R - CI * f1I) + fenZERO, + f0 * (CR * f1I + CI * f1R) + fenZERO] + vbMng(self, "DEL", "Done assembling base expression.", 25) + return forcingTerm + def b(self, mu : paramVal = [], der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" - raise Exception("Deprecated.") mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) + bDEIMCoeffs = None for j in range(derI, nbsTot): if bs[j] is None: vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) - if j < self.nbs: - fRe, fIm = self.getForcingTerm(self.mu0) - cRe, cIm = self.getExtraFactorB(self.mu0, j) - cfRe, cfIm = cRe * fRe - cIm * fIm, cRe * fIm + cIm * fRe - else: - cfRe, cfIm = fenZERO, fenZERO - parsRe = self.iterReduceQuadratureDegree(zip([cfRe], - ["forcingTermDer{}Real".format(j)])) - parsIm = self.iterReduceQuadratureDegree(zip([cfIm], - ["forcingTermDer{}Imag".format(j)])) - L0Re = fen.dot(cfRe, self.v) * fen.dx - L0Im = fen.dot(cfIm, self.v) * fen.dx - DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + if bDEIMCoeffs is None: + bDEIM = np.empty((self.nbs, self.spacedim()), + dtype = np.complex) + muDEIM = 3. * np.linspace(0., 1., self.nbs // 2 + 1) ** 2. + muDEIM = np.concatenate((-muDEIM[:0:-1], muDEIM)) + for jj, muD in enumerate(muDEIM): + fRe, fIm = self.getForcingTerm(muD) + parsRe = self.iterReduceQuadratureDegree(zip([fRe], + ["forcingTerm{}Real".format(jj)])) + parsIm = self.iterReduceQuadratureDegree(zip([fIm], + ["forcingTerm{}Imag".format(jj)])) + LR = fen.dot(fRe, self.v) * fen.dx + LI = fen.dot(fIm, self.v) * fen.dx + DBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) - b = (fenics2Vector(L0Re, parsRe, DirichletBC0, 1) - + 1.j * fenics2Vector(L0Im, parsIm, DirichletBC0, 1)) + bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1) + + 1.j * fenics2Vector(LI, parsIm, DBC0, 1)) + bDEIMCoeffs = (fit(muDEIM / 3., bDEIM, self.nbs - 1).T + * np.power(3., - np.arange(self.nbs))).T + b = bDEIMCoeffs[j] if homogeneized: - Ader = self.A(self.mu0, hashI(j, self.npar)) + Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py index 79568ae..25dcbe3 100644 --- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py +++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py @@ -1,363 +1,363 @@ # 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.base.vector_problem_engine_base import \ VectorProblemEngineBase from rrompy.utilities.base.types import Np1D, List, ScOp, paramVal from rrompy.solver.fenics import (fenZERO, fenZEROS, fenONE, L2InverseNormMatrix, elasticNormMatrix, elasticDualNormMatrix) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.parameter import checkParameter from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['LinearElasticityProblemEngine'] class LinearElasticityProblemEngine(VectorProblemEngineBase): """ Solver for generic linear elasticity problems. - div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = f in \Omega u = u0 on \Gamma_D \partial_nu = g1 on \Gamma_N \partial_nu + h u = g2 on \Gamma_R Attributes: verbosity: Verbosity level. BCManager: Boundary condition manager. V: Real vector FE space. u: Generic vector trial functions for variational form evaluation. v: Generic vector test functions for variational form evaluation. As: Scipy sparse array representation (in CSC format) of As. bs: Numpy array representation of bs. energyNormMatrix: Scipy sparse matrix representing inner product over V. energyNormDualMatrix: Scipy sparse matrix representing inner product over V'. dualityMatrix: Scipy sparse matrix representing duality V-V'. energyNormPartialDualMatrix: Scipy sparse matrix representing dual inner product between Riesz representers V-V. liftedDirichletDatum: Dofs of Dirichlet datum lifting. mu0BC: Mu value of last Dirichlet datum lifting. degree_threshold: Threshold for ufl expression interpolation degree. lambda_: Value of lambda_. mu_: Value of mu_. forcingTerm: Value of f. DirichletDatum: Value of u0. NeumannDatum: Value of g1. RobinDatumG: Value of g2. RobinDatumH: Value of h. DirichletBoundary: Function handle to \Gamma_D. NeumannBoundary: Function handle to \Gamma_N. RobinBoundary: Function handle to \Gamma_R. ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries). A0: Scipy sparse array representation (in CSC format) of A0. b0: Numpy array representation of b0. dsToBeSet: Whether ds needs to be set. """ _energyDualNormCompress = None def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.lambda_ = fenONE self.mu_ = fenONE self.mu0 = checkParameter(mu0) self.npar = self.mu0.shape[1] self.forcingTerm = fenZEROS(self.V.mesh().topology().dim()) self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim()) self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim()) self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim()) self.RobinDatumH = fenZERO @property def V(self): """Value of V.""" return self._V @V.setter def V(self, V): VectorProblemEngineBase.V.fset(self, V) self.forcingTerm = fenZEROS(self.V.mesh().topology().dim()) self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim()) self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim()) self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim()) self.dsToBeSet = True @property def lambda_(self): """Value of lambda_.""" return self._lambda_ @lambda_.setter def lambda_(self, lambda_): self.resetAs() if not isinstance(lambda_, (list, tuple,)): lambda_ = [lambda_, fenZERO] self._lambda_ = lambda_ @property def mu_(self): """Value of mu_.""" return self._mu_ @mu_.setter def mu_(self, mu_): self.resetAs() if not isinstance(mu_, (list, tuple,)): mu_ = [mu_, fenZERO] self._mu_ = mu_ @property def forcingTerm(self): """Value of f.""" return self._forcingTerm @forcingTerm.setter def forcingTerm(self, forcingTerm): self.resetbs() if not isinstance(forcingTerm, (list, tuple,)): forcingTerm = [forcingTerm, fenZEROS(self.V.mesh().topology().dim())] self._forcingTerm = forcingTerm @property def DirichletDatum(self): """Value of u0.""" return self._DirichletDatum @DirichletDatum.setter def DirichletDatum(self, DirichletDatum): self.resetbs() if not isinstance(DirichletDatum, (list, tuple,)): DirichletDatum = [DirichletDatum, fenZEROS(self.V.mesh().topology().dim())] self._DirichletDatum = DirichletDatum @property def NeumannDatum(self): """Value of g1.""" return self._NeumannDatum @NeumannDatum.setter def NeumannDatum(self, NeumannDatum): self.resetbs() if not isinstance(NeumannDatum, (list, tuple,)): NeumannDatum = [NeumannDatum, fenZEROS(self.V.mesh().topology().dim())] self._NeumannDatum = NeumannDatum @property def RobinDatumG(self): """Value of g2.""" return self._RobinDatumG @RobinDatumG.setter def RobinDatumG(self, RobinDatumG): self.resetbs() if not isinstance(RobinDatumG, (list, tuple,)): RobinDatumG = [RobinDatumG, fenZEROS(self.V.mesh().topology().dim())] self._RobinDatumG = RobinDatumG @property def RobinDatumH(self): """Value of h.""" return self._RobinDatumH @RobinDatumH.setter def RobinDatumH(self, RobinDatumH): self.resetAs() if not isinstance(RobinDatumH, (list, tuple,)): RobinDatumH = [RobinDatumH, fenZERO] self._RobinDatumH = RobinDatumH @property def DirichletBoundary(self): """Function handle to DirichletBoundary.""" return self.BCManager.DirichletBoundary @DirichletBoundary.setter def DirichletBoundary(self, DirichletBoundary): self.resetAs() self.resetbs() self.BCManager.DirichletBoundary = DirichletBoundary @property def NeumannBoundary(self): """Function handle to NeumannBoundary.""" return self.BCManager.NeumannBoundary @NeumannBoundary.setter def NeumannBoundary(self, NeumannBoundary): self.resetAs() self.resetbs() self.dsToBeSet = True self.BCManager.NeumannBoundary = NeumannBoundary @property def RobinBoundary(self): """Function handle to RobinBoundary.""" return self.BCManager.RobinBoundary @RobinBoundary.setter def RobinBoundary(self, RobinBoundary): self.resetAs() self.resetbs() self.dsToBeSet = True self.BCManager.RobinBoundary = RobinBoundary def autoSetDS(self): """Set FEniCS boundary measure based on boundary function handles.""" if self.dsToBeSet: vbMng(self, "INIT", "Initializing boundary measures.", 20) NB = self.NeumannBoundary RB = self.RobinBoundary boundary_markers = fen.MeshFunction("size_t", self.V.mesh(), self.V.mesh().topology().dim() - 1) NB.mark(boundary_markers, 0) RB.mark(boundary_markers, 1) self.ds = fen.Measure("ds", domain = self.V.mesh(), subdomain_data = boundary_markers) self.dsToBeSet = False vbMng(self, "DEL", "Done initializing boundary measures.", 20) def buildEnergyNormForm(self): """ Build sparse matrix (in CSR format) representative of scalar product. """ vbMng(self, "INIT", "Assembling energy matrix.", 20) self.energyNormMatrix = elasticNormMatrix(self.V, self.lambda_[0], self.mu_[0]) vbMng(self, "DEL", "Done assembling energy matrix.", 20) def buildEnergyNormDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product. """ vbMng(self, "INIT", "Assembling energy dual matrix.", 20) self.energyNormDualMatrix = elasticDualNormMatrix( self.V, self.lambda_[0], self.mu_[0], compressRank = self._energyDualNormCompress) vbMng(self, "DEL", "Done assembling energy dual matrix.", 20) def buildDualityPairingForm(self): """Build sparse matrix (in CSR format) representative of duality.""" vbMng(self, "INIT", "Assembling duality matrix.", 20) self.dualityMatrix = L2InverseNormMatrix( self.V, solverType = self._solver, solverArgs = self._solverArgs, compressRank = self._dualityCompress) vbMng(self, "DEL", "Done assembling duality matrix.", 20) def buildEnergyNormPartialDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product without duality. """ vbMng(self, "INIT", "Assembling energy partial dual matrix.", 20) self.energyNormPartialDualMatrix = elasticDualNormMatrix( self.V, self.lambda_[0], self.mu_[0], compressRank = self._energyDualNormCompress, duality = False) vbMng(self, "DEL", "Done assembling energy partial dual matrix.", 20) def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) self.autoSetDS() if derI <= 0 and self.As[0] is None: vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(self.V.mesh().topology().dim()), self.DirichletBoundary) lambda_Re, lambda_Im = self.lambda_ mu_Re, mu_Im = self.mu_ hRe, hIm = self.RobinDatumH termNames = ["lambda_", "mu_", "RobinDatumH"] parsRe = self.iterReduceQuadratureDegree(zip( [lambda_Re, mu_Re, hRe], [x + "Real" for x in termNames])) parsIm = self.iterReduceQuadratureDegree(zip( [lambda_Im, mu_Re, hIm], [x + "Imag" for x in termNames])) epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u)) sigma = lambda u, l_, m_: ( l_ * fen.div(u) * fen.Identity(u.geometric_dimension()) + 2. * m_ * epsilon(u)) a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re), epsilon(self.v)) * fen.dx + hRe * fen.inner(self.u, self.v) * self.ds(1)) a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im), epsilon(self.v)) * fen.dx + hIm * fen.inner(self.u, self.v) * self.ds(1)) self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1) + 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu != self.mu0BC: self.liftDirichletData(self.mu) fenZEROSEff = fenZEROS(self.V.mesh().topology().dim()) for j in range(derI, nbsTot): if bs[j] is None: self.autoSetDS() vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) if j == 0: u0Re, u0Im = self.DirichletDatum fRe, fIm = self.forcingTerm g1Re, g1Im = self.NeumannDatum g2Re, g2Im = self.RobinDatumG else: u0Re, u0Im = fenZEROSEff, fenZEROSEff fRe, fIm = fenZEROSEff, fenZEROSEff g1Re, g1Im = fenZEROSEff, fenZEROSEff g2Re, g2Im = fenZEROSEff, fenZEROSEff termNames = ["forcingTerm", "NeumannDatum", "RobinDatumG"] parsRe = self.iterReduceQuadratureDegree(zip( [fRe, g1Re, g2Re], [x + "Real" for x in termNames])) parsIm = self.iterReduceQuadratureDegree(zip( [fIm, g1Im, g2Im], [x + "Imag" for x in termNames])) L0Re = (fen.inner(fRe, self.v) * fen.dx + fen.inner(g1Re, self.v) * self.ds(0) + fen.inner(g2Re, self.v) * self.ds(1)) L0Im = (fen.inner(fIm, self.v) * fen.dx + fen.inner(g1Im, self.v) * self.ds(0) + fen.inner(g2Im, self.v) * self.ds(1)) DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary) b = (fenics2Vector(L0Re, parsRe, DBCR, 1) + 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1)) if homogeneized: - Ader = self.A(self.mu0, hashI(j, self.npar)) + Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/tests/hfengines/helmholtz_internal.py b/tests/hfengines/helmholtz_internal.py index cd3cc3d..d6fde8c 100644 --- a/tests/hfengines/helmholtz_internal.py +++ b/tests/hfengines/helmholtz_internal.py @@ -1,101 +1,97 @@ # 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 pytest import os, shutil import numpy as np from rrompy.hfengines.linear_problem import ( HelmholtzSquareBubbleDomainProblemEngine, HelmholtzSquareBubbleProblemEngine, HelmholtzSquareTransmissionProblemEngine) def test_helmholtz_square_io(): solver = HelmholtzSquareBubbleProblemEngine(kappa = 4, theta = 1., n = 20, verbosity = 0) mu = 5 uh = solver.solve(mu)[0] assert np.isclose(solver.norm(uh), 145.0115, rtol = 1e-3) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), 1.19934e-11, rtol = 1e-1) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), solver.norm(solver.residual(uh, mu, duality = False)[0], dual = True, duality = False), rtol = 1e-1) if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache") filesOut = [x for x in os.listdir("./.pytest_cache") if (x[-4:] == ".pvd" and x[:9] == "outSquare")] filesOutData = [x for x in os.listdir("./.pytest_cache") if (x[-4:] == ".vtu" and x[:9] == "outSquare")] for fileOut in filesOut: os.remove("./.pytest_cache/" + fileOut) for fileOut in filesOutData: os.remove("./.pytest_cache/" + fileOut) solver.outParaview(uh, what = ["MESH", "ABS"], filename = ".pytest_cache/outSquare", forceNewFile = False) filesOut = [x for x in os.listdir("./.pytest_cache") if (x[-4:] == ".pvd" and x[:9] == "outSquare")] filesOutData = [x for x in os.listdir("./.pytest_cache") if (x[-4:] == ".vtu" and x[:9] == "outSquare")] assert len(filesOut) == 1 assert len(filesOutData) == 1 os.remove("./.pytest_cache/" + filesOut[0]) os.remove("./.pytest_cache/" + filesOutData[0]) def test_helmholtz_transmission_io(): solver = HelmholtzSquareTransmissionProblemEngine(nT = 1, nB = 2, theta = np.pi * 40 / 180, kappa = 4., n = 20, verbosity = 0) mu = 5. uh = solver.solve(mu)[0] assert np.isclose(solver.norm(uh), 138.6609, rtol = 1e-2) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), 3.7288565e-12, rtol = 1e-1) if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache") solver.outParaviewTimeDomain(uh, omega = mu, filename = ".pytest_cache/outTrans", forceNewFile = False, folder = True) filesOut = [x for x in os.listdir("./.pytest_cache/outTrans") if (x[-4:] == ".pvd" and x[:8] == "outTrans")] filesOutData = [x for x in os.listdir("./.pytest_cache/outTrans") if (x[-4:] == ".vtu" and x[:8] == "outTrans")] assert len(filesOut) == 1 assert len(filesOutData) == 20 shutil.rmtree("./.pytest_cache/outTrans") -@pytest.mark.xfail(reason = ("no graphical interface and deprecated RHS " - "derivative computation")) def test_helmholtz_domain_io(): solver = HelmholtzSquareBubbleDomainProblemEngine(kappa = 4, theta = 1., n = 10, mu0 = 1.5, verbosity = 0) mu = 1.5 uh = solver.solve(mu)[0] if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache") solver.plot(uh, save = "./.pytest_cache/outDomain", show = False) filesOut = [x for x in os.listdir("./.pytest_cache") if (x[-4:] == ".eps" and x[:9] == "outDomain")] assert len(filesOut) == 1 os.remove("./.pytest_cache/" + filesOut[0]) assert np.isclose(solver.norm(uh), 10.07843, rtol = 1e-2) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), 6.14454989e-13, rtol = 1e-1) - diff --git a/tests/hfengines/laplace.py b/tests/hfengines/laplace.py index 86f0c92..cae0297 100644 --- a/tests/hfengines/laplace.py +++ b/tests/hfengines/laplace.py @@ -1,44 +1,42 @@ # 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 pytest import numpy as np from rrompy.hfengines.linear_problem import LaplaceDiskGaussian from rrompy.hfengines.linear_problem.bidimensional import LaplaceDiskGaussian2 -@pytest.mark.xfail(reason = "deprecated RHS derivative computation") def test_laplace_disk(): solver = LaplaceDiskGaussian(n = 20, verbosity = 0) mu = 1.5 solver.setSolver("BICG", {"tol" : 1e-15}) uh = solver.solve(mu)[0] assert np.isclose(solver.norm(uh), 1.053403077447029, rtol = 1e-2) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), - 3.4591353e-08, rtol = 1e-1) + 4.66728e-10, rtol = 1e-1) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), solver.norm(solver.residual(uh, mu, duality = False)[0], dual = True, duality = False), rtol = 1e-1) def test_laplace_disk_2(): solver = LaplaceDiskGaussian2(n = 20, verbosity = 0) mu = [[0., 1.5]] mu = [0., 1.5] uh = solver.solve(mu)[0] assert np.isclose(solver.norm(uh), 1.0534030774205372, rtol = 1e-2) assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True), 2.40043363e-08, rtol = 1e-1)