diff --git a/examples/1_symmetric_disk/symmetric_disk.py b/examples/1_symmetric_disk/symmetric_disk.py index 9511194..ca9e8a7 100644 --- a/examples/1_symmetric_disk/symmetric_disk.py +++ b/examples/1_symmetric_disk/symmetric_disk.py @@ -1,65 +1,65 @@ import numpy as np from symmetric_disk_engine import SymmetricDiskEngine as engine from rrompy.reduction_methods.standard import (RationalInterpolant as RI, ReducedBasis as RB) from rrompy.reduction_methods.greedy import (RationalInterpolantGreedy as RIG, ReducedBasisGreedy as RBG) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS ks = [10., 20.] k0, n = np.mean(np.power(ks, 2.)) ** .5, 150 solver = engine(k0, n) k = 12. -method = "RI" -method = "RB" -method = "RI_GREEDY" -method = "RB_GREEDY" - -if method == "RI": - params = {'N':39, 'M':39, 'S':40, 'POD':True, 'polybasis':"CHEBYSHEV", - 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} - algo = RI -if method == "RB": - params = {'R':40, 'S':40, 'POD':True, - 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} - algo = RB -if method == "RI_GREEDY": - params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, - 'sampler':QS(ks, "UNIFORM", scalingExp = 2.), - 'errorEstimatorKind':"DISCREPANCY"} - algo = RIG -if method == "RB_GREEDY": - params = {'S':10, 'POD':True, 'greedyTol':1e-2, - 'sampler':QS(ks, "UNIFORM", scalingExp = 2.)} - algo = RBG +for method in ["RI", "RB", "RI_GREEDY", "RB_GREEDY"]: + print("Testing {} method".format(method)) -approx = algo(solver, mu0 = k0, approx_state = True, - approxParameters = params, verbosity = 20) -if len(method) == 2: - approx.setupApprox() -else: - approx.setupApprox("ALL") - -approx.plotApprox(k, name = 'u_app') -approx.plotHF(k, name = 'u_HF') -approx.plotErr(k, name = 'err') -approx.plotRes(k, name = 'res') -normErr = approx.normErr(k)[0] -normSol = approx.normHF(k)[0] -normRes = approx.normRes(k)[0] -normRHS = approx.normRHS(k)[0] -print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + if method == "RI": + params = {'N':39, 'M':39, 'S':40, 'POD':True, 'polybasis':"CHEBYSHEV", + 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} + algo = RI + if method == "RB": + params = {'R':40, 'S':40, 'POD':True, + 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} + algo = RB + if method == "RI_GREEDY": + params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, + 'sampler':QS(ks, "UNIFORM", scalingExp = 2.), + 'errorEstimatorKind':"DISCREPANCY"} + algo = RIG + if method == "RB_GREEDY": + params = {'S':10, 'POD':True, 'greedyTol':1e-2, + 'sampler':QS(ks, "UNIFORM", scalingExp = 2.)} + algo = RBG + + approx = algo(solver, mu0 = k0, approx_state = True, + approxParameters = params, verbosity = 20) + if len(method) == 2: + approx.setupApprox() + else: + approx.setupApprox("LAST") + + approx.plotApprox(k, name = 'u_app') + approx.plotHF(k, name = 'u_HF') + approx.plotErr(k, name = 'err') + approx.plotRes(k, name = 'res') + normErr = approx.normErr(k)[0] + normSol = approx.normHF(k)[0] + normRes = approx.normRes(k)[0] + normRHS = approx.normRHS(k)[0] + print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( normSol, normErr, normErr / normSol)) -print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + + if method[:2] == "RI": + poles, residues = approx.getResidues() + if method[:2] == "RB": + poles = approx.getPoles() + print("Poles:\n{}".format(poles)) + if method[:2] == "RI": + for pol, res in zip(poles, residues.T): + solver.plot(res) + print("pole = {:.5e}".format(pol)) -if method[:2] == "RI": - poles, residues = approx.getResidues() -if method[:2] == "RB": - poles = approx.getPoles() -print("Poles:\n{}".format(poles)) -if method[:2] == "RI": - for pol, res in zip(poles, residues.T): - solver.plot(res) - print("pole = {:.5e}".format(pol)) + print("\n") diff --git a/examples/2_double_slit/double_slit.py b/examples/2_double_slit/double_slit.py index 4551bde..e379bc3 100644 --- a/examples/2_double_slit/double_slit.py +++ b/examples/2_double_slit/double_slit.py @@ -1,59 +1,62 @@ import numpy as np from double_slit_engine import DoubleSlitEngine as engine from rrompy.reduction_methods.standard import RationalInterpolant as RI from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG from rrompy.parameter.parameter_sampling import QuadratureSampler as QS from rrompy.solver.fenics import interp_project ks = [10., 15.] k0, n = np.mean(ks), 150 solver = engine(k0, n) k = 11. -method = "RI" -#method = "RI_GREEDY" +for method in ["RI", "RI_GREEDY"]: + print("Testing {} method".format(method)) -if method == "RI": - params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV", - 'sampler':QS(ks, "CHEBYSHEV")} - algo = RI -if method == "RI_GREEDY": - params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, - 'sampler':QS(ks, "UNIFORM"), 'errorEstimatorKind':"LOOK_AHEAD", - 'trainSetGenerator':QS(ks, "CHEBYSHEV")} - algo = RIG + if method == "RI": + params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV", + 'sampler':QS(ks, "CHEBYSHEV")} + algo = RI + if method == "RI_GREEDY": + params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, + 'sampler':QS(ks, "UNIFORM"), + 'errorEstimatorKind':"LOOK_AHEAD", + 'trainSetGenerator':QS(ks, "CHEBYSHEV")} + algo = RIG + + approx = algo(solver, mu0 = k0, approx_state = True, + approxParameters = params, verbosity = 20) + if len(method) == 2: + approx.setupApprox() + else: + approx.setupApprox("LAST") -approx = algo(solver, mu0 = k0, approx_state = True, - approxParameters = params, verbosity = 20) -if len(method) == 2: - approx.setupApprox() -else: - approx.setupApprox("ALL") - -approx.plotApprox(k, name = 'u_app') -approx.plotHF(k, name = 'u_HF') -approx.plotErr(k, name = 'err') -approx.plotRes(k, name = 'res') -normErr = approx.normErr(k)[0] -normSol = approx.normHF(k)[0] -normRes = approx.normRes(k)[0] -normRHS = approx.normRHS(k)[0] -print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + approx.plotApprox(k, name = 'u_app') + approx.plotHF(k, name = 'u_HF') + approx.plotErr(k, name = 'err') + approx.plotRes(k, name = 'res') + normErr = approx.normErr(k)[0] + normSol = approx.normHF(k)[0] + normRes = approx.normRes(k)[0] + normRHS = approx.normRHS(k)[0] + print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( normSol, normErr, normErr / normSol)) -print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + + uIncR, uIncI = solver.getDirichletValues(k) + uIncR = interp_project(uIncR, solver.V) + uIncI = interp_project(uIncI, solver.V) + uInc = np.array(uIncR.vector()) + 1.j * np.array(uIncI.vector()) + uEx = approx.getHF(k)[0] - uInc + uApp = approx.getApprox(k)[0] - uInc + solver.plot(uEx, name = 'uex_tot') + solver.plot(uApp, name = 'u_app_tot') + + poles, residues = approx.getResidues() + print("Poles:\n{}".format(poles)) + for pol, res in zip(poles, residues.T): + solver.plot(res) + print("pole = {:.5e}".format(pol)) -uIncR, uIncI = solver.getDirichletValues(k) -uIncR = interp_project(uIncR, solver.V) -uIncI = interp_project(uIncI, solver.V) -uInc = np.array(uIncR.vector()) + 1.j * np.array(uIncI.vector()) -uEx = approx.getHF(k)[0] - uInc -uApp = approx.getApprox(k)[0] - uInc -solver.plot(uEx, name = 'uex_tot') -solver.plot(uApp, name = 'u_app_tot') - -poles, residues = approx.getResidues() -print("Poles:\n{}".format(poles)) -for pol, res in zip(poles, residues.T): - solver.plot(res) - print("pole = {:.5e}".format(pol)) + print("\n") diff --git a/examples/2_double_slit/double_slit_engine.py b/examples/2_double_slit/double_slit_engine.py index 7a6bf67..ce049ea 100644 --- a/examples/2_double_slit/double_slit_engine.py +++ b/examples/2_double_slit/double_slit_engine.py @@ -1,53 +1,55 @@ import numpy as np import ufl import fenics as fen import mshr -from rrompy.hfengines.base import LinearNonAffineEngine +from rrompy.utilities.base import affine_construct, nonaffine_construct from rrompy.hfengines.linear_problem import ScatteringProblemEngine from rrompy.utilities.numerical import hashDerivativeToIdx as hashD from rrompy.solver.fenics import fenZERO, fenics2Vector -class DoubleSlitEngine(LinearNonAffineEngine, ScatteringProblemEngine): +class DoubleSlitEngine(ScatteringProblemEngine): def __init__(self, k0:float, n:int): super().__init__(mu0 = [k0]) self._affinePoly = False delta, eps = .1, .01 mesh = mshr.generate_mesh( mshr.Circle(fen.Point(0., 0.), 3.) - mshr.Rectangle(fen.Point(-3., -delta), fen.Point(-.75, delta)) - mshr.Rectangle(fen.Point(-.5, -delta), fen.Point(.5, delta)) - mshr.Rectangle(fen.Point(.75, -delta), fen.Point(3., delta)), n) self.V = fen.FunctionSpace(mesh, "P", 1) self.DirichletBoundary = lambda x, on_boundary: (on_boundary and np.abs(x[1]) <= delta and np.abs(x[1]) > delta - eps) self.NeumannBoundary = lambda x, on_boundary: (on_boundary and np.abs(x[1]) <= delta - eps) self.RobinBoundary = "REST" def getDirichletValues(self, mu = []): mu = self.checkParameter(mu) x, y = fen.SpatialCoordinate(self.V.mesh())[:] c, s = .5, - .5 * 3. ** .5 muR, muI = np.real(mu[0])[0], np.imag(mu[0])[0] mod = - muI * (c * x + s * y) angle = muR * (c * x + s * y) DR = fen.exp(mod) * fen.cos(angle) DI = fen.exp(mod) * fen.sin(angle) DR = ufl.conditional(ufl.ge(y, 0), DR, fenZERO) DI = ufl.conditional(ufl.ge(y, 0), DI, fenZERO) return DR, DI + @affine_construct def A(self, mu = [], der = 0): return ScatteringProblemEngine.A(self, mu, der) + @nonaffine_construct def b(self, mu = [], der = 0): derI = hashD(der) if hasattr(der, "__len__") else der if derI < 0: return self.baselineb() if derI > 0: raise Exception("Derivatives not implemented.") fen0 = fen.inner(fenZERO, self.v) * fen.dx DR, DI = self.getDirichletValues(mu) DBCR = fen.DirichletBC(self.V, DR, self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, DI, self.DirichletBoundary) return (fenics2Vector(fen0, {}, DBCR, 1) + 1.j * fenics2Vector(fen0, {}, DBCI, 1)) diff --git a/examples/3_sector_angle/sector_angle.py b/examples/3_sector_angle/sector_angle.py index 5c53c1c..d693ea1 100644 --- a/examples/3_sector_angle/sector_angle.py +++ b/examples/3_sector_angle/sector_angle.py @@ -1,80 +1,90 @@ import numpy as np import matplotlib.pyplot as plt from sector_angle_engine import SectorAngleEngine as engine from rrompy.reduction_methods.pivoted import ( RationalInterpolantPivoted as RIP, RationalInterpolantGreedyPivoted as RIGP) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS -ks, ts = [10., 15.], [.35, .65] +ks, ts = [10., 15.], [.4, .6] k0, t0, n = np.mean(np.power(ks, 2.)) ** .5, np.mean(ts), 50 solver = engine(k0, t0, n) murange = [[ks[0], ts[0]], [ks[-1], ts[-1]]] mu = [12., .535] -method = "RI" -#method = "RI_GREEDY" +for method in ["RI", "RI_GREEDY"]: + print("Testing {} method".format(method)) -if method == "RI": - params = {'N':19, 'M':19, 'S':20, 'MMarginal':2, 'SMarginal':7, 'POD':True, - 'polybasis':"CHEBYSHEV", 'polybasisMarginal':"MONOMIAL", - 'matchingWeight':1., 'cutOffTolerance': 1., - 'samplerPivot':QS(ks, "CHEBYSHEV", 2.), - 'samplerMarginal':QS(ts, "UNIFORM")} - algo = RIP -if method == "RI_GREEDY": - params = {'S':10, 'MMarginal':2, 'SMarginal':7, 'POD':True, - 'polybasis':"LEGENDRE", 'polybasisMarginal':"MONOMIAL", - 'matchingWeight':1., 'cutOffTolerance': 1., - 'samplerPivot':QS(ks, "UNIFORM", 2.), - 'greedyTol':1e-2, 'errorEstimatorKind':"INTERPOLATORY", - 'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.), - 'samplerMarginal':QS(ts, "UNIFORM")} - algo = RIGP + if method == "RI": + params = {'N':19, 'M':19, 'S':20, 'MMarginal':6, 'SMarginal':7, + 'POD':True, 'polybasis':"CHEBYSHEV", + 'polybasisMarginal':"MONOMIAL", 'matchingWeight':1., + 'cutOffTolerance': 1., + 'samplerPivot':QS(ks, "CHEBYSHEV", 2.), + 'samplerMarginal':QS(ts, "UNIFORM")} + algo = RIP + if method == "RI_GREEDY": + params = {'S':10, 'MMarginal':6, 'SMarginal':7, 'POD':True, + 'polybasis':"LEGENDRE", 'polybasisMarginal':"MONOMIAL", + 'matchingWeight':1., 'cutOffTolerance': 1., + 'samplerPivot':QS(ks, "UNIFORM", 2.), + 'greedyTol':1e-3, 'errorEstimatorKind':"INTERPOLATORY", + 'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.), + 'samplerMarginal':QS(ts, "UNIFORM")} + algo = RIGP + + approx = algo([0], solver, mu0 = [k0, t0], approx_state = True, + approxParameters = params, verbosity = 20) + if len(method) == 2: + approx.setupApprox() + else: + approx.setupApprox("LAST") -approx = algo([0], solver, mu0 = [k0, t0], approx_state = True, - approxParameters = params, verbosity = 20) -if len(method) == 2: - approx.setupApprox() -else: - approx.setupApprox("ALL") - -approx.plotApprox(mu, name = 'u_app') -approx.plotHF(mu, name = 'u_HF') -approx.plotErr(mu, name = 'err') -approx.plotRes(mu, name = 'res') -normErr = approx.normErr(mu)[0] -normSol = approx.normHF(mu)[0] -normRes = approx.normRes(mu)[0] -normRHS = approx.normRHS(mu)[0] -print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + approx.plotApprox(mu, name = 'u_app') + approx.plotHF(mu, name = 'u_HF') + approx.plotErr(mu, name = 'err') + approx.plotRes(mu, name = 'res') + normErr = approx.normErr(mu)[0] + normSol = approx.normHF(mu)[0] + normRes = approx.normRes(mu)[0] + normRHS = approx.normRHS(mu)[0] + print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( normSol, normErr, normErr / normSol)) -print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + + verb = approx.verbosity + approx.verbosity = 0 + tspace = np.linspace(ts[0], ts[-1], 100) + for j, t in enumerate(tspace): + pls = approx.getPoles([None, t]) + pls[np.abs(np.imag(pls ** 2.)) > 1e-5] = np.nan + if j == 0: poles = np.empty((len(tspace), len(pls))) + poles[j] = np.real(pls) + approx.verbosity = verb + plt.figure(figsize = (12, 5)) + plt.subplot(121) + plt.plot(poles, tspace) + plt.ylim(ts) + plt.xlabel('mu_1') + plt.ylabel('mu_2') + plt.grid() + plt.subplot(122) + plt.plot(poles, tspace) + plt.xlim(ks) + plt.ylim(ts) + plt.xlabel('mu_1') + plt.ylabel('mu_2') + plt.grid() + plt.show() + plt.close() + + poles, residues = approx.getResidues([None, mu[1]]) + for pol, res in zip(poles, residues.T): + if (np.abs(np.imag(pol ** 2.)) < 1e-5 and np.real(pol) >= ks[0] + and np.real(pol) <= ks[-1]): + solver.plot(res) + print("pole = {:.5e}".format(pol)) -verb = approx.verbosity -approx.verbosity = 0 -tspace = np.linspace(ts[0], ts[-1], 100) -for j, t in enumerate(tspace): - pls = approx.getPoles([None, t]) - pls[np.abs(np.imag(pls ** 2.)) > 1e-5] = np.nan - if j == 0: poles = np.empty((len(tspace), len(pls))) - poles[j] = np.real(pls) -approx.verbosity = verb -plt.figure() -plt.plot(poles, tspace) -plt.xlim(ks) -plt.ylim(ts) -plt.xlabel('mu_1') -plt.ylabel('mu_2') -plt.grid() -plt.show() -plt.close() - -poles, residues = approx.getResidues([None, mu[1]]) -for pol, res in zip(poles, residues.T): - if (np.abs(np.imag(pol ** 2.)) < 1e-5 and np.real(pol) >= ks[0] - and np.real(pol) <= ks[-1]): - solver.plot(res) - print("pole = {:.5e}".format(pol)) + print("\n") diff --git a/examples/3_sector_angle/sector_angle_engine.py b/examples/3_sector_angle/sector_angle_engine.py index 94f79b9..9356f52 100644 --- a/examples/3_sector_angle/sector_angle_engine.py +++ b/examples/3_sector_angle/sector_angle_engine.py @@ -1,41 +1,43 @@ import numpy as np import fenics as fen import mshr -from rrompy.hfengines.base import LinearNonAffineEngine +from rrompy.utilities.base import nonaffine_construct from rrompy.hfengines.linear_problem import HelmholtzProblemEngine -class SectorAngleEngine(LinearNonAffineEngine, HelmholtzProblemEngine): +class SectorAngleEngine(HelmholtzProblemEngine): def __init__(self, k0:float, t0:float, n:int): super().__init__(mu0 = [k0, t0]) self._affinePoly = False self.npar = 2 self.rescalingExp = [2., 1.] mesh = mshr.generate_mesh( mshr.Circle(fen.Point(0., 0.), 1.) - mshr.Rectangle(fen.Point(-1., -1.), fen.Point(0., 1.)) - mshr.Rectangle(fen.Point(-1., -1.), fen.Point(1., 0.)), n) self.V = fen.FunctionSpace(mesh, "P", 1) x, y = fen.SpatialCoordinate(mesh)[:] self.forcingTerm = [fen.exp(x + y) * (1. - x ** 2. - y ** 2.), fen.exp(x - y) * (1. - x ** 2. - y ** 2.)] self._tBoundary = np.nan def setBoundary(self, t:float): while hasattr(t, "__len__"): t = t[0] if not np.isclose(t, self._tBoundary): eps = 1e-2 self._tBoundary = t self.DirichletBoundary = lambda x, on_boundary: ( on_boundary and x[0] >= eps and x[1] <= eps + np.sin(t * np.pi / 2.)) self.NeumannBoundary = "REST" + @nonaffine_construct def A(self, mu = [], der = 0): mu = self.checkParameter(mu) self.setBoundary(mu(1)) return HelmholtzProblemEngine.A(self, mu, der) + @nonaffine_construct def b(self, mu = [], der = 0): mu = self.checkParameter(mu) self.setBoundary(mu(1)) return HelmholtzProblemEngine.b(self, mu, der)