diff --git a/examples/base/solver.py b/examples/base/solver.py index b34ceae..3426d8f 100644 --- a/examples/base/solver.py +++ b/examples/base/solver.py @@ -1,72 +1,91 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import ( - HelmholtzSquareTransmissionProblemEngine as HSTPE) -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE -from rrompy.hfengines.scipy import ( - HelmholtzCavityScatteringProblemEngine as HCSPE) +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzCavityScatteringProblemEngine as HCSPE +from rrompy.hfengines.linear_mixed_problem import \ + MixedHelmholtzProblemEngine as MHPE -testNo = 3 +testNo = 5 verb = 0 if testNo == 1: solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, verbosity = verb) mu = 12.**.5 uh = solver.solve(mu) solver.plotmesh() print(solver.norm(uh)) solver.plot(uh) solver.plot(solver.residual(uh, mu), 'res') ########### -elif testNo == [2, -2]: +elif testNo in [2, -2]: solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4., n = 50, verbosity = verb) mu = 4. uref = solver.liftDirichletData(mu) if testNo > 0: uh = solver.solve(mu) utot = uh - uref else: utot = solver.solve(mu, homogeneized = True) uh = utot + uref print(solver.norm(uh)) print(solver.norm(uref)) solver.plot(uh) solver.plot(uref, name = 'u_Dir') solver.plot(utot, name = 'u_tot') solver.plot(solver.residual(uh, mu), 'res') solver.plot(solver.residual(utot, mu, homogeneized = True), 'res_tot') ########### elif testNo in [3, -3]: solver = HBSPE(R = 5, kappa = 12**.5, theta = - np.pi * 60 / 180, n = 30, verbosity = verb) mu = 12**.5 uref = solver.liftDirichletData(mu) if testNo > 0: uh = solver.solve(mu) utot = uh - uref else: utot = solver.solve(mu, homogeneized = True) uh = utot + uref solver.plotmesh() print(solver.norm(uh)) print(solver.norm(utot)) solver.plot(uh) solver.plot(utot, name = 'u_tot') solver.plot(solver.residual(uh, mu), 'res') solver.plot(solver.residual(utot, mu, homogeneized = True), 'res_tot') ########### elif testNo == 4: solver = HCSPE(kappa = 5, n = 30, verbosity = verb) mu = 10 uh = solver.solve(mu) solver.plotmesh() print(solver.norm(uh)) solver.plot(uh) solver.plot(solver.residual(uh, mu), 'res') + +########### +elif testNo == 5: + solver = MHPE(verbosity = verb) + solver.NeumannBoundary = "ALL" + import fenics as fen + solver.forcingTerm = fen.Expression(("10*exp(-(pow(x[0] - 0.5, 2) + " + "pow(x[1] - 0.5, 2)) / 0.02)"), + degree = 2) + solver.NeumannDatum = fen.Expression("-sin(5*x[0])", degree = 2) + uh = solver.solve(np.pi) + solver.plotmesh() + print(solver.norm(uh)) + solver.plot(uh) + solver.plot(solver.residual(uh, np.pi), 'res') diff --git a/examples/greedy/squareScatteringHomog.py b/examples/from_papers/greedy_internalBox.py similarity index 70% copy from examples/greedy/squareScatteringHomog.py copy to examples/from_papers/greedy_internalBox.py index d037701..2caff05 100644 --- a/examples/greedy/squareScatteringHomog.py +++ b/examples/from_papers/greedy_internalBox.py @@ -1,97 +1,110 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as \ - HCSPE +import fenics as fen +from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB +dim = 3 + verb = 2 timed = False algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" -errorEstimatorKind = "SIMPLIFIED" -#errorEstimatorKind = "EXACT" -k0s = np.linspace(5, 12, 100) -k0 = np.mean(k0s) +k0s = np.power(np.linspace(500 ** 2., 2250 ** 2., 200), .5) +k0 = np.mean(np.power(k0s, 2.)) ** .5 kl, kr = min(k0s), max(k0s) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, - 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis, - 'errorEstimatorKind':errorEstimatorKind} + 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis} -if timed: - verb = 0 +if dim == 2: + mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(.1, .15), 10, 15) + x, y = fen.SpatialCoordinate(mesh)[:] + f = fen.exp(- 1e2 * (x + y)) +else:#if dim == 3: + mesh = fen.BoxMesh(fen.Point(0., 0., 0.), fen.Point(.1, .15, .25), 4, 6,10) + x, y, z = fen.SpatialCoordinate(mesh)[:] + f = fen.exp(- 1e2 * (x + y + z)) -solver = HCSPE(kappa = 5, n = 10, verbosity = verb) +solver = HPE(verbosity = verb) solver.omega = np.real(k0) +solver.V = fen.FunctionSpace(mesh, "P", 3) +solver.refractionIndex = fen.Constant(1. / 54.6) +solver.forcingTerm = f +solver.NeumannBoundary = "ALL" + +######### + if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 +kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j])) / approx.estNormer.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus), - 1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') + 1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/greedy/squareScatteringAirfoil.py b/examples/from_papers/greedy_scatteringAirfoil.py similarity index 97% rename from examples/greedy/squareScatteringAirfoil.py rename to examples/from_papers/greedy_scatteringAirfoil.py index 58b9024..b90b954 100644 --- a/examples/greedy/squareScatteringAirfoil.py +++ b/examples/from_papers/greedy_scatteringAirfoil.py @@ -1,145 +1,146 @@ import numpy as np import fenics as fen import ufl -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HSP +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HSP from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB from rrompy.utilities.base.fenics import fenONE verb = 2 timed = False algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" homog = True homog = False k0s = np.linspace(5, 20, 25) k0 = np.mean(k0s) kl, kr = min(k0s), max(k0s) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis} ######### PI = np.pi R = 2 def Dboundary(x, on_boundary): return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R kappa = 10 theta = PI * - 45 / 180. mu = 1.1 epsilon = .1 mesh = fen.Mesh('../data/mesh/airfoil.xml') c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] u0R = - fen.cos(kappa * (c * x + s * y)) u0I = - fen.sin(kappa * (c * x + s * y)) checkReal = x**2-x+y**2 rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25 phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2 phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2 kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/ ((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.) )**.5 - mu kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/ ((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.) )**.5 - mu Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1 Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1 cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE) c_F = fen.Constant(.1) cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE) c_F = fen.Constant(.1) cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F) cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F) a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF) solver = HSP(R, np.real(k0), theta, n = 1, verbosity = verb, degree_threshold = 8) solver.omega = np.real(k0) solver.V = fen.FunctionSpace(mesh, "P", 3) solver.diffusivity = a solver.DirichletBoundary = Dboundary solver.RobinBoundary = "REST" solver.DirichletDatum = [u0R, u0I] ######### if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j], homogeneized=homog)) / approx.estNormer.norm(approx.getRHS(k0s[j], homogeneized=homog))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus), 1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/pod/membraneTaylor.py b/examples/from_papers/pod_membraneTaylor.py similarity index 97% rename from examples/pod/membraneTaylor.py rename to examples/from_papers/pod_membraneTaylor.py index dee9594..d6d03c9 100644 --- a/examples/pod/membraneTaylor.py +++ b/examples/from_papers/pod_membraneTaylor.py @@ -1,96 +1,96 @@ import fenics as fen import numpy as np -from rrompy.hfengines.base import HelmholtzProblemEngine as HPE +from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper verb = 0 test = "poles" test = "error" k0 = 10 ktars = np.linspace(78**.5, 122**.5, 50) def boundaryNeumann(x, on_boundary): return on_boundary and x[1] > .25 and x[0] > 0.995 and x[0] < 1.005 meshname = '../data/mesh/crack_coarse.xml' #meshname = '../data/mesh/crack_fine.xml' mesh = fen.Mesh(meshname) x, y = fen.SpatialCoordinate(mesh)[:] x0, y0 = .5, .5 Rr, Ri = .1, .1 forcingTerm = fen.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Rr**2) solver = HPE(verbosity = verb) solver.omega = np.real(k0) solver.V = fen.FunctionSpace(mesh, "P", 3) solver.forcingTerm = forcingTerm solver.NeumannBoundary = boundaryNeumann solver.DirichletBoundary = 'rest' if test == "poles": appPoles = {} Emax = 13 params = {'N':6, 'M':0, 'E':6, 'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} approxPade = TP(solver, mu0 = k0, approxParameters = params, verbosity = verb) for E in range(6, Emax + 1): approxPade.E = E appPoles[E] = np.sort(approxPade.getPoles()) a = fen.dot(fen.grad(solver.u), fen.grad(solver.v)) * fen.dx A = fen.assemble(a) fen.DirichletBC(solver.V, fen.Constant(0.), solver.DirichletBoundary).apply(A) AMat = fen.as_backend_type(A).mat() Ar, Ac, Av = AMat.getValuesCSR() import scipy.sparse as scsp A = scsp.csr_matrix((Av, Ac, Ar), shape = AMat.size) m = fen.dot(solver.u, solver.v) * fen.dx M = fen.assemble(m) fen.DirichletBC(solver.V, fen.Constant(0.), solver.DirichletBoundary).apply(M) MMat = fen.as_backend_type(M).mat() Mr, Mc, Mv = MMat.getValuesCSR() import scipy.sparse as scsp M = scsp.csr_matrix((Mv, Mc, Mr), shape = MMat.size) poles = scsp.linalg.eigs(A, k = 7, M = M, sigma = 100., return_eigenvectors = False) II = np.argsort(np.abs(poles - k0)) poles = poles[II] print('Exact', end = ': ') [print('{},{}'.format(np.real(x), np.imag(x)), end = ',') for x in poles] print() for E in range(6, Emax + 1): print(E, end = ': ') [print('{},{}'.format(np.real(x), np.imag(x)), end = ',')\ for x in np.sort(appPoles[E])] print() elif test == "error": M0 = 5 Emax = 8 params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} paramsSetsPade = [None] * (Emax - M0 + 1) for M in range(M0, Emax + 1): paramsSetsPade[M - M0] = {'N':M0 + 1, 'M':M, 'E':max(M, M0 + 1)} approxPade = TP(solver, mu0 = k0, approxParameters = params, verbosity = verb) sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep('../data/output/membrane_error.dat', outputs = 'ALL') sweeper.plot(filenamePade, ['muRe'], ['normHF', 'normApp'], ['M'], onePlot = True) sweeper.plot(filenamePade, ['muRe'], ['normErr'], ['M']) diff --git a/examples/pod/airfoil.py b/examples/from_papers/pod_scatteringAirfoil.py similarity index 98% rename from examples/pod/airfoil.py rename to examples/from_papers/pod_scatteringAirfoil.py index 2b6169e..f489dc2 100644 --- a/examples/pod/airfoil.py +++ b/examples/from_papers/pod_scatteringAirfoil.py @@ -1,212 +1,213 @@ from copy import deepcopy as copy import numpy as np import fenics as fen import ufl -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HSP +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HSP from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from rrompy.utilities.base.fenics import fenONE from operator import itemgetter def subdict(d, ks): return dict(zip(ks, itemgetter(*ks)(d))) verb = 0 #################### homog = True #homog = False #################### test = "solve" test = "Taylor" test = "Lagrange" test = "TaylorSweep" test = "LagrangeSweep" plotSamples = True k0 = 10 kLeft, kRight = 8 + 0.j, 12 + 0.j ktar = 11 ktars = np.linspace(8, 12, 21) + 0.j PI = np.pi R = 2 def Dboundary(x, on_boundary): return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R kappa = 10 theta = PI * - 45 / 180. mu = 1.1 epsilon = .1 mesh = fen.Mesh('../data/mesh/airfoil.xml') c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] u0R = - fen.cos(kappa * (c * x + s * y)) u0I = - fen.sin(kappa * (c * x + s * y)) checkReal = x**2-x+y**2 rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25 phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2 phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2 kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/ ((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.) )**.5 - mu kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/ ((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.) )**.5 - mu Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1 Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1 cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE) c_F = fen.Constant(.1) cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE) c_F = fen.Constant(.1) cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F) cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F) a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF) ### solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb, degree_threshold = 8) solver.V = fen.FunctionSpace(mesh, "P", 3) solver.diffusivity = a solver.DirichletBoundary = Dboundary solver.RobinBoundary = "REST" solver.DirichletDatum = [u0R, u0I] ### if test == "solve": uinc = solver.liftDirichletData(k0) if homog: uhtot = solver.solve(k0, homogeneized = homog) uh = uhtot + uinc else: uh = solver.solve(k0, homogeneized = homog) uhtot = uh - uinc print(solver.norm(uh)) print(solver.norm(uhtot)) solver.plot(fen.project(a, solver.V).vector(), what = 'Real', name = 'a') solver.plot(uinc, what = 'Real', name = 'u_inc') solver.plot(uh, what = 'ABS') solver.plot(uhtot, what = 'ABS', name = 'u + u_inc') elif test in ["Taylor", "Lagrange"]: if test == "Taylor": params = {'N':8, 'M':8, 'R':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD']) parRB = subdict(params, ['R', 'E', 'sampleType', 'POD']) approxPade = TP(solver, mu0 = k0, approxParameters = parPade, verbosity = verb, homogeneized = homog) approxRB = TRB(solver, mu0 = k0, approxParameters = parRB, verbosity = verb, homogeneized = homog) else: params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"CHEBYSHEV", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler', 'basis']) parRB = subdict(params, ['R', 'S', 'POD', 'sampler']) approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parPade, verbosity = verb, homogeneized = homog) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parRB, verbosity = verb, homogeneized = homog) approxPade.setupApprox() approxRB.setupApprox() if plotSamples: approxPade.plotSamples() approxPade.plotHF(ktar, name = 'u_HF') approxPade.plotApp(ktar, name = 'u_Pade''') approxPade.plotErr(ktar, name = 'err_Pade''') approxPade.plotRes(ktar, name = 'res_Pade''') approxRB.plotApp(ktar, name = 'u_RB') approxRB.plotErr(ktar, name = 'err_RB') approxRB.plotRes(ktar, name = 'res_RB') HFNorm, RHSNorm = approxPade.normHF(ktar), approxPade.normRHS(ktar) PadeRes, PadeErr = approxPade.normRes(ktar), approxPade.normErr(ktar) RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar) print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm)) print('PadeRes:\t{}\nPadeErr:\t{}'.format(PadeRes, PadeErr)) print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr)) print('\nPoles Pade'':') print(approxPade.getPoles()) elif test in ["TaylorSweep", "LagrangeSweep"]: if test == "TaylorSweep": shift = 10 nsets = 2 stride = 3 Emax = stride * (nsets - 1) + shift + 1 params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift + 1, 'M':stride*i+shift+1, 'E':stride*i+shift + 1} paramsSetsRB[i] = {'E':stride*i+shift + 1,'R':stride*i+shift + 2} approxPade = TP(solver, mu0 = kappa,approxParameters = params, verbosity = verb, homogeneized = homog) approxRB = TRB(solver, mu0 = kappa, approxParameters = params, verbosity = verb, homogeneized = homog) else: shift = 10 nsets = 2 stride = 3 Smax = stride * (nsets - 1) + shift + 2 paramsPade = {'S':Smax, 'POD':True, 'basis':"CHEBYSHEV", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} paramsRB = copy(paramsPade) paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift+1, 'S': stride*i+shift+2} paramsSetsRB[i] = {'R': stride*i+shift+2, 'S': stride*i+shift+2} approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsPade, verbosity = verb, homogeneized = homog) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsRB, verbosity = verb, homogeneized = homog) sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') homogMSG = "" if not homog: homogMSG = "Non" filenamebase = '../data/output/airfoil' + test[:-5] + homogMSG + "Homog" sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase +'Pade.dat') sweeper.ROMEngine = approxRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase +'RB.dat') if test == "TaylorSweep": constr = ['E'] else: constr = ['S'] sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normHF', 'normApp'], constr, onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB"]) sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normResRel'], constr, save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB"]) sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normErrRel'], constr, save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB"]) diff --git a/examples/greedy/squareScatteringHomog.py b/examples/greedy/parametricDomain.py similarity index 80% copy from examples/greedy/squareScatteringHomog.py copy to examples/greedy/parametricDomain.py index d037701..516a4fb 100644 --- a/examples/greedy/squareScatteringHomog.py +++ b/examples/greedy/parametricDomain.py @@ -1,97 +1,99 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as \ - HCSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleDomainProblemEngine as HSBDPE from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB verb = 2 -timed = False +timed = True algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" errorEstimatorKind = "SIMPLIFIED" -#errorEstimatorKind = "EXACT" +errorEstimatorKind = "EXACT" -k0s = np.linspace(5, 12, 100) -k0 = np.mean(k0s) +k0s = np.power(np.linspace(9, 19, 100), .5) +k0 = np.mean(np.power(k0s, 2.)) ** .5 kl, kr = min(k0s), max(k0s) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis, 'errorEstimatorKind':errorEstimatorKind} if timed: verb = 0 -solver = HCSPE(kappa = 5, n = 10, verbosity = verb) +solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, mu0 = k0, + degree_threshold = 15, verbosity = verb) solver.omega = np.real(k0) + if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j])) / approx.estNormer.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() -plt.plot(k0s, norm) -plt.plot(k0s, normApp, '--') -plt.plot(np.real(approx.mus), - 1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') +plt.semilogy(k0s, norm) +plt.semilogy(k0s, normApp, '--') +plt.semilogy(np.real(approx.mus), + 1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/greedy/squareBubbleHomog.py b/examples/greedy/squareBubbleHomog.py index be15aa4..ef8a3b1 100644 --- a/examples/greedy/squareBubbleHomog.py +++ b/examples/greedy/squareBubbleHomog.py @@ -1,108 +1,109 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB from rrompy.utilities.base import squareResonances verb = 2 timed = True algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" errorEstimatorKind = "SIMPLIFIED" errorEstimatorKind = "EXACT" -k0s = np.power(np.linspace(95, 149, 200), .5) +k0s = np.power(np.linspace(95, 149, 500), .5) k0s = np.power(np.linspace(95, 129, 20), .5) k0 = np.mean(np.power(k0s, 2.)) ** .5 kl, kr = min(k0s), max(k0s) polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5)) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis, 'errorEstimatorKind':errorEstimatorKind} if timed: verb = 0 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 30, verbosity = verb) solver.omega = np.real(k0) if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 from matplotlib import pyplot as plt normApp = np.zeros_like(k0s) norm = np.zeros_like(k0s) res = np.zeros_like(k0s) err = np.zeros_like(k0s) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j])) / approx.estNormer.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.semilogy(k0s, norm) plt.semilogy(k0s, normApp, '--') plt.semilogy(polesexact, 2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.') plt.semilogy(np.real(approx.mus), 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(polesexact, 2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.semilogy(polesexact, 2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.plot(np.real(polesexact), np.imag(polesexact), 'm.') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/greedy/squareScatteringHomog.py b/examples/greedy/squareScatteringHomog.py index d037701..5996be4 100644 --- a/examples/greedy/squareScatteringHomog.py +++ b/examples/greedy/squareScatteringHomog.py @@ -1,97 +1,97 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as \ - HCSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzCavityScatteringProblemEngine as HCSPE from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB verb = 2 -timed = False +timed = True algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" errorEstimatorKind = "SIMPLIFIED" -#errorEstimatorKind = "EXACT" +errorEstimatorKind = "EXACT" k0s = np.linspace(5, 12, 100) k0 = np.mean(k0s) kl, kr = min(k0s), max(k0s) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, 'greedyTol':1e-2, 'nTestPoints':2, 'basis':polyBasis, 'errorEstimatorKind':errorEstimatorKind} if timed: verb = 0 solver = HCSPE(kappa = 5, n = 10, verbosity = verb) solver.omega = np.real(k0) if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j])) / approx.estNormer.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus), 1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/greedy/squareTransmissionNonHomog.py b/examples/greedy/squareTransmissionNonHomog.py index 5803a52..e0674d6 100644 --- a/examples/greedy/squareTransmissionNonHomog.py +++ b/examples/greedy/squareTransmissionNonHomog.py @@ -1,108 +1,108 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine \ - as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangeRBGreedy as RB from rrompy.reduction_methods.lagrange_greedy import \ ApproximantLagrangePadeGreedy as Pade timed = False verb = 2 algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" homog = True #homog = False errorEstimatorKind = "SIMPLIFIED" errorEstimatorKind = "EXACT" k0s = np.power(np.linspace(4, 15, 100), .5) k0 = np.mean(np.power(k0s, 2.)) ** .5 kl, kr = min(k0s), max(k0s) rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0, 'greedyTol':1e-2, 'nTestPoints':5, 'basis':polyBasis, 'errorEstimatorKind':errorEstimatorKind} solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4., n = 20, verbosity = verb) solver.omega = np.real(k0) if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) else: approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApp(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estNormer.norm(approx.getRes(k0s[j], homogeneized=homog)) / approx.estNormer.norm(approx.getRHS(k0s[j], homogeneized=homog))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) polesApp = approx.getPoles() polesApp = polesApp[np.abs(np.imag(polesApp)) < 1e-3] plt.figure() plt.semilogy(k0s, norm) plt.semilogy(k0s, normApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') plt.semilogy(np.real(polesApp), 2.*np.max(norm)*np.ones_like(polesApp, dtype = float), 'k.') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, res) plt.semilogy(k0s, resApp, '--') plt.semilogy(np.real(approx.mus), 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') plt.semilogy(np.real(polesApp), 2.*np.max(resApp)*np.ones_like(polesApp, dtype = float), 'k.') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.semilogy(np.real(polesApp), 2.*np.max(err)*np.ones_like(polesApp, dtype = float), 'k.') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/pod/LagrangePoles.py b/examples/pod/LagrangePoles.py index bbda022..b3af094 100644 --- a/examples/pod/LagrangePoles.py +++ b/examples/pod/LagrangePoles.py @@ -1,46 +1,47 @@ from matplotlib import pyplot as plt import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from rrompy.utilities.base import squareResonances verb = 0 ks = [1, 46 ** .5] solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, verbosity = verb) k0 = np.mean(np.power(ks, 2.)) ** .5 k0 = 3.46104724 solver.omega = np.real(k0) rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) sampler = QS(ks, "UNIFORM", rescaling, rescalingInv) nsets = 15 paramsPade = {'S':2, 'POD':True, 'basis':"LEGENDRE", 'sampler':sampler} approx = Pade(solver, mu0 = k0, approxParameters = paramsPade, verbosity = verb) poles = [None] * (nsets - 1) polesexact = np.unique(np.power(squareResonances(ks[0]**2., ks[1]**2., False), .5)) for i in range(1, nsets): print("N = {}".format(4 * i)) approx.approxParameters = {'N': 4 * i, 'M': 4 * i, 'S': 4 * i + 1} approx.setupApprox() poles[i - 1] = approx.getPoles() for i in range(1, nsets): plt.figure() plt.plot(np.real(poles[i - 1]), np.imag(poles[i - 1]), 'kx') plt.plot(polesexact, np.zeros_like(polesexact), 'm.') plt.plot(k0, 0, 'r*') plt.xlim(ks) plt.ylim((ks[0] - ks[1]) / 2., (ks[1] - ks[0]) / 2.) plt.title("N = {}, Neff = {}".format(4 * i, len(poles[i - 1]))) plt.grid() plt.show() plt.close() \ No newline at end of file diff --git a/examples/pod/LagrangeSweep.py b/examples/pod/LagrangeSweep.py index c291cd0..9025425 100644 --- a/examples/pod/LagrangeSweep.py +++ b/examples/pod/LagrangeSweep.py @@ -1,108 +1,108 @@ from copy import copy import numpy as np -#from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as \ -# HBSPE -from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as \ - HBSPE +#from rrompy.hfengines.linear_problem import \ +# HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzCavityScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from rrompy.utilities.parameter_sampling import WarpingFunction as WF from rrompy.utilities.parameter_sampling import WarpingSampler as WS verb = 0 npoints = 100 homog = True homog = False LSratio = 2. / 3. sampling = "Uniform" #sampling = "Cheby" #sampling = "SinC" assert LSratio <= 1. + np.finfo(float).eps ks = [10 + 0.j, 14 + 0.j] solver = HBSPE(kappa = 3, n = 15) solver.omega = np.real(np.mean(ks)) mutars = np.linspace(9, 15, npoints) homogMSG = "Homog" if not homog: homogMSG = "Non" + homogMSG filenamebase = '../data/output/ScatteringSquareLSweep' #filenamebase = '../data/plots/LagrangeScatteringSquare1p5' #filenamebase = filenamebase + sampling + "/HelmholtzBoxLSweep" + homogMSG k0 = np.mean(ks) shift = 7 shift = np.int(8 / LSratio - 1) nsets = 3 stride = np.int(8 / LSratio) Smax = stride * (nsets - 1) + shift + 2 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) if sampling == "Uniform": polyBasis = "MONOMIAL" sampler = QS(ks, "UNIFORM", rescaling, rescalingInv) if sampling == "Cheby": polyBasis = "CHEBYSHEV" sampler = QS(ks, "CHEBYSHEV", rescaling, rescalingInv) if sampling == "SinC": polyBasis = "MONOMIAL" warping = WF(call = lambda x: (x - 2. * (1. - LSratio) / np.pi * np.sin(np.pi * x)), repr = "x-{}*sin(pi*x)".format(2. * (1. - LSratio) / np.pi)) sampler = WS(ks, warping, rescaling, rescalingInv) paramsPade = {'S':Smax, 'POD':True, 'basis':polyBasis, 'sampler':sampler} paramsRB = copy(paramsPade) paramsPoly = copy(paramsPade) paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': np.int(LSratio * (stride * i + shift + 1)), 'M': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} paramsSetsRB[i] = {'R': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} paramsSetsPoly[i] = {'N': 0, 'M': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} appPade = Pade(solver, mu0 = k0, approxParameters = paramsPade, verbosity = verb, homogeneized = homog) appRB = RB(solver, mu0 = k0, approxParameters = paramsRB, verbosity = verb, homogeneized = homog) appPoly = Pade(solver, mu0 = k0, approxParameters = paramsPoly, verbosity = verb, homogeneized = homog) sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx') sweeper.ROMEngine = appPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = appRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = appPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['S'], onePlot = True, save = filenamebase + 'Norm', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['S'], save = filenamebase + 'Res', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['S'], save = filenamebase + 'Err', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) diff --git a/examples/pod/LagrangeSweepUnstable.py b/examples/pod/LagrangeSweepUnstable.py index 27836aa..e05297e 100644 --- a/examples/pod/LagrangeSweepUnstable.py +++ b/examples/pod/LagrangeSweepUnstable.py @@ -1,75 +1,76 @@ from copy import copy import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import ManualSampler as MS npoints = 50 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 30, verbosity = 0) mutars = np.linspace(6**.5, 16**.5, npoints) filenamebase = '../data/output/HelmholtzBubbleLagrange' k0 = np.mean(mutars) rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) paramsChoices = { 'Stable': {'S':4, 'POD':True, 'basis':"MONOMIAL", 'sampler':MS([6**.5, 16**.5], [[14**.5], [12**.5], [9**.5], [11**.5]])}, 'Unstable': {'S':4, 'POD':True, 'basis':"MONOMIAL", 'sampler':MS([6**.5, 16**.5], [[8**.5], [10**.5], [13**.5], [11**.5]])} } j = 0 filenamePade = [None] * len(paramsChoices.keys()) filenameRB = [None] * len(paramsChoices.keys()) for typeP in paramsChoices.keys(): paramsPade = paramsChoices[typeP] paramsRB = copy(paramsPade) paramsSetsPade = [{'N': 3, 'M': 3}] paramsSetsRB = [{'R': 4}] appPade = Pade(solver, mu0 = k0, approxParameters = paramsPade, verbosity = 0) appRB = RB(solver, mu0 = k0, approxParameters = paramsRB, verbosity = 0) sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx') sweeper.ROMEngine = appPade sweeper.params = paramsSetsPade filenamePade[j] = sweeper.sweep(filenamebase + 'Pade{}.dat'.format(typeP), verbose = 0) sweeper.ROMEngine = appRB sweeper.params = paramsSetsRB filenameRB[j] = sweeper.sweep(filenamebase + 'RB{}.dat'.format(typeP), verbose = 0) j += 1 print("Pade'") sweeper.plotCompare(filenamePade, ['muRe'], ['normHF', 'normApp'], ['S'], onePlot = True, save = filenamebase + 'PadeNorm', saveFormat = "png", labels = list(paramsChoices.keys())) sweeper.plotCompare(filenamePade, ['muRe'], ['normResRel'], ['S'], save = filenamebase + 'PadeRes', saveFormat = "png", labels = list(paramsChoices.keys())) sweeper.plotCompare(filenamePade, ['muRe'], ['normErrRel'], ['S'], save = filenamebase + 'PadeErr', saveFormat = "png", labels = list(paramsChoices.keys())) print("RB") sweeper.plotCompare(filenameRB, ['muRe'], ['normHF', 'normApp'], ['S'], onePlot = True, save = filenamebase + 'RBNorm', saveFormat = "png", labels = list(paramsChoices.keys())) sweeper.plotCompare(filenameRB, ['muRe'], ['normResRel'], ['S'], save = filenamebase + 'RBRes', saveFormat = "png", labels = list(paramsChoices.keys())) sweeper.plotCompare(filenameRB, ['muRe'], ['normErrRel'], ['S'], save = filenamebase + 'RBErr', saveFormat = "png", labels = list(paramsChoices.keys())) diff --git a/examples/pod/PadeLagrange.py b/examples/pod/PadeLagrange.py index f4eb9c2..3ea0ca6 100644 --- a/examples/pod/PadeLagrange.py +++ b/examples/pod/PadeLagrange.py @@ -1,111 +1,113 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine \ - as HSTPE -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.utilities.parameter_sampling import QuadratureSampler as QS testNo = 1 verb = 5 polyBasis = "CHEBYSHEV" polyBasis = "LEGENDRE" #polyBasis = "MONOMIAL" homog = True homog = False if testNo == 1: k0s = np.power([10 + 0.j, 14 + 0.j], .5) k0 = np.mean(k0s) ktar = (11 + 0.j) ** .5 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'N':4, 'M':3, 'S':5, 'POD':True, 'basis':polyBasis, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 2: k0s = [3.85 + 0.j, 4.15 + 0.j] k0 = np.mean(k0s) ktar = 4 + 0.j rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'N':8, 'M':9, 'S':10, 'POD':True, 'basis':polyBasis, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 3: k0s = [2, 5] k0 = np.mean(k0s) ktar = 4.5 - .1j params = {'N':10, 'M':10, 'S':11, 'POD':True, 'basis':polyBasis, 'sampler':QS(k0s, "CHEBYSHEV")} solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) diff --git a/examples/pod/PadeTaylor.py b/examples/pod/PadeTaylor.py index ed114a5..16c8187 100644 --- a/examples/pod/PadeTaylor.py +++ b/examples/pod/PadeTaylor.py @@ -1,97 +1,99 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import ( - HelmholtzSquareTransmissionProblemEngine as HSTPE) -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade testNo = 1 verb = 5 homog = True #homog = False if testNo == 1: params = {'N':4, 'M':3, 'E':4, 'sampleType':'Arnoldi', 'POD':True} k0 = 12 ** .5 ktar = 10.5 ** .5 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 2: params = {'N':6, 'M':7, 'E':7, 'sampleType':'Arnoldi', 'POD':True} k0 = 16 ** .5 ktar = 15 ** .5 solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo in [3, 4]: if testNo == 3: params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True} else: params = {'N':7, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} k0 = 3 ktar = 4.+0.j solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) diff --git a/examples/pod/RBLagrange.py b/examples/pod/RBLagrange.py index ac5474f..6677973 100644 --- a/examples/pod/RBLagrange.py +++ b/examples/pod/RBLagrange.py @@ -1,99 +1,102 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB from rrompy.utilities.parameter_sampling import QuadratureSampler as QS testNo = 3 verb = 0 homog = True homog = False if testNo == 1: k0s = np.power([10 + 0.j, 14 + 0.j], .5) k0 = np.mean(k0s) ktar = (11 + .5j) ** .5 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'S':5, 'R':4, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 2: k0s = [3.85 + 0.j, 4.15 + 0.j] k0 = np.mean(k0s) ktar = 4 + .15j rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'S':10, 'R':9, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 3: k0s = [2, 5] k0 = np.mean(k0s) ktar = 4.5 - 0.j params = {'S':15, 'R':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")} solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) diff --git a/examples/pod/RBTaylor.py b/examples/pod/RBTaylor.py index 00b729a..c1df6f1 100644 --- a/examples/pod/RBTaylor.py +++ b/examples/pod/RBTaylor.py @@ -1,90 +1,93 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB testNo = 4 verb = 0 homog = True #homog = False if testNo == 1: params = {'E':4, 'R':4, 'sampleType':'Arnoldi', 'POD':True} k0 = 12 ** .5 ktar = 10.5 ** .5 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 2: params = {'E':7, 'R':7, 'sampleType':'Arnoldi', 'POD':True} k0 = 16**.5 ktar = 15**.5 solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 3., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo in [3, 4]: if testNo == 3: params = {'E':8, 'sampleType':'Krylov', 'POD':True} else: params = {'E':8, 'sampleType':'Arnoldi', 'POD':True} k0 = 3 ktar = 4.25+.5j solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) diff --git a/examples/pod/TaylorPoles.py b/examples/pod/TaylorPoles.py index d59d092..d2c5c71 100644 --- a/examples/pod/TaylorPoles.py +++ b/examples/pod/TaylorPoles.py @@ -1,68 +1,69 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB from rrompy.utilities.base import squareResonances verb = 0 k0 = (12+0.j) ** .5 Nmin, Nmax = 2, 10 Nvals = np.arange(Nmin, Nmax + 1, 2) params = {'N':Nmin, 'M':0, 'Emax':Nmax, 'POD':True, 'sampleType':'Arnoldi'} #, 'robustTol':1e-14} #boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x) # - np.real(z0)) #cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 25, verbosity = verb) solver.omega = np.real(k0) approxP = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb)#, # equilibration = True, cleanupParameters = cleanupParameters) approxR = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) rP, rE = [None] * len(Nvals), [None] * len(Nvals) verbose = 1 for j, N in enumerate(Nvals): if verbose > 0: print('N = E = {}'.format(N)) approxP.approxParameters = {'N':N, 'E':N} approxR.approxParameters = {'R':N, 'E':N} if verbose > 1: print(approxP.approxParameters) print(approxR.approxParameters) rP[j] = approxP.getPoles() rE[j] = approxR.getPoles() if verbose > 2: print(rP) print(rE) from matplotlib import pyplot as plt plotRows = int(np.ceil(len(Nvals) / 3)) fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows)) for j, N in enumerate(Nvals): i1, i2 = int(np.floor(j / 3)), j % 3 axes[i1, i2].set_title('N = E = {}'.format(N)) axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb', label="Pade'", markersize = 8) axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), 'Pr', label="RB", markersize = 8) axes[i1, i2].axhline(linewidth=1, color='k') xmin, xmax = axes[i1, i2].get_xlim() height = (xmax - xmin) / 2. res = np.power(squareResonances(xmin**2., xmax**2., False), .5) axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4) axes[i1, i2].plot(np.real(k0), np.imag(k0), 'om', markersize = 5) axes[i1, i2].plot(np.real(k0) * np.ones(2), 1.5 * height * np.arange(-1, 3, 2), '--m') axes[i1, i2].grid() axes[i1, i2].set_xlim(xmin, xmax) axes[i1, i2].set_ylim(- height, height) p = axes[i1, i2].legend() plt.tight_layout() for j in range((len(Nvals) - 1) % 3 + 1, 3): axes[plotRows - 1, j].axis('off') diff --git a/examples/pod/TaylorSweep.py b/examples/pod/TaylorSweep.py index 3cf37e6..32878d3 100644 --- a/examples/pod/TaylorSweep.py +++ b/examples/pod/TaylorSweep.py @@ -1,77 +1,79 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE -from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.hfengines.linear_problem import \ + HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper testNo = 1 verb = 0 homog = True homog = False k0 = (12 + 0.j) ** .5 npoints = 100 shift = 5 nsets = 2 stride = 2 Emax = stride * (nsets - 1) + shift + 1 if testNo == 1: solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 15, verbosity = verb) solver.omega = np.real(k0) params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True} ktars = np.linspace(7**.5, 16**.5, npoints) filenamebase = '../data/output/HelmholtzBubbleTSweep' homog = False elif testNo == 2: solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10, verbosity = verb) solver.omega = np.real(k0) params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} ktars = np.linspace(7**.5, 17**.5, npoints) homogMSG = "" if not homog: homogMSG = "Non" filenamebase = '../data/output/HelmholtzBoxTSweep' + homogMSG + "Homog" paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift+1, 'E':stride*i+shift+1} paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+2} paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift+1, 'E':stride*i+shift+1} appPade = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) appRB = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) appPoly = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb, homogeneized = homog) sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') sweeper.ROMEngine = appPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = appRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = appPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['E'], onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['E'], save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['E'], save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) diff --git a/examples/pod/laplaceGaussianTaylor.py b/examples/pod/laplaceGaussianTaylor.py index bd16914..8b0c564 100644 --- a/examples/pod/laplaceGaussianTaylor.py +++ b/examples/pod/laplaceGaussianTaylor.py @@ -1,100 +1,100 @@ import numpy as np -from rrompy.hfengines.scipy import LaplaceDiskGaussian as LDG +from rrompy.hfengines.linear_problem import LaplaceDiskGaussian as LDG from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from operator import itemgetter def subdict(d, ks): return dict(zip(ks, itemgetter(*ks)(d))) testNo = 3 verb = 0 if testNo == 1: mu = 4. solver = LDG(n = 40, verbosity = verb) uh = solver.solve(mu) solver.plotmesh() print(solver.norm(uh)) solver.plot(uh) ############ if testNo == 2: params = {'N':8, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} # params = {'N':8, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True} mu0 = 0. solver = LDG(n = 20, degree_threshold = 15, verbosity = verb) approxP = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) paramsRB = subdict(params, ['E', 'sampleType', 'POD']) approxR = RB(solver, mu0 = mu0, approxParameters = paramsRB, verbosity = verb) approxP.setupApprox() approxR.setupApprox() # approxP.plotSamples() mutar = 3.25 approxP.plotHF(mutar, name = 'u_HF') approxP.plotApp(mutar, name = 'u_Pade''') approxR.plotApp(mutar, name = 'u_RB') approxP.plotErr(mutar, name = 'err_Pade''') approxR.plotErr(mutar, name = 'err_RB') solNorm = approxP.normHF(mutar) appPErr = approxP.normErr(mutar) appRErr = approxR.normErr(mutar) print(('SolNorm:\t{}\nErrRelP:\t{}\nErrRelR:\t{}').format(solNorm, appPErr / solNorm, appRErr / solNorm)) ############ elif testNo == 3: mu0 = 0. mutars = np.linspace(0., 4., 60) solver = LDG(n = 10, degree_threshold = 15, verbosity = verb) shift = 2 nsets = 6 stride = 2 Emax = stride * (nsets - 1) + shift params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} params = {'Emax':Emax, 'sampleType':'Krylov', 'POD':False} paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift, 'M':stride*i+shift, 'E':stride*i+shift} paramsSetsRB[i] = {'E':stride*i+shift} paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift, 'E':stride*i+shift} approxPade = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) approxRB = RB(solver, mu0 = mu0, approxParameters = params, verbosity = verb) approxPoly = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) filenamebase = '../data/output/LapGTaylor' sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx') sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = approxRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = approxPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['E'], onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['E'], save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['E'], save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) diff --git a/examples/pod/parametricDomain.py b/examples/pod/parametricDomain.py index f5c5dc4..ce2df61 100644 --- a/examples/pod/parametricDomain.py +++ b/examples/pod/parametricDomain.py @@ -1,103 +1,104 @@ import numpy as np -from rrompy.hfengines.scipy import HelmholtzSquareBubbleDomainProblemEngine as HSBDPE +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleDomainProblemEngine as HSBDPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from operator import itemgetter def subdict(d, ks): return dict(zip(ks, itemgetter(*ks)(d))) testNo = 3 verb = 0 if testNo == 1: mu = 7 ** .5 solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, mu0 = mu, degree_threshold = 15, verbosity = verb) uh = solver.solve(mu) solver.plotmesh() print(solver.norm(uh)) solver.plot(uh) ############ if testNo == 2: params = {'N':8, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} # params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True} mu0 = 7 ** .5 mutar = (7. + .1j) ** .5 solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, mu0 = mu0, degree_threshold = 15, verbosity = verb) approxP = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) paramsRB = subdict(params, ['E', 'sampleType', 'POD']) approxR = RB(solver, mu0 = mu0, approxParameters = paramsRB, verbosity = verb) approxP.setupApprox() approxR.setupApprox() # approxP.plotSamples() approxP.plotHF(mutar, name = 'u_HF') approxP.plotApp(mutar, name = 'u_Pade''') approxR.plotApp(mutar, name = 'u_RB') approxP.plotErr(mutar, name = 'err_Pade''') approxR.plotErr(mutar, name = 'err_RB') solNorm = approxP.normHF(mutar) appPErr = approxP.normErr(mutar) appRErr = approxR.normErr(mutar) print(('SolNorm:\t{}\nErrRelP:\t{}\nErrRelR:\t{}').format(solNorm, appPErr / solNorm, appRErr / solNorm)) print('\nPoles Pade'':') print(approxP.getPoles()) ############ elif testNo == 3: mu0 = 14 ** .5 mutars = np.linspace(9**.5, 19**.5, 100) solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, mu0 = mu0, degree_threshold = 15, verbosity = verb) shift = 10 nsets = 1 stride = 2 Emax = stride * (nsets - 1) + shift params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift, 'M':stride*i+shift, 'E':stride*i+shift} paramsSetsRB[i] = {'E':stride*i+shift} paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift, 'E':stride*i+shift} approxPade = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) approxRB = RB(solver, mu0 = mu0, approxParameters = params, verbosity = verb) approxPoly = Pade(solver, mu0 = mu0, approxParameters = params, verbosity = verb) filenamebase = '../data/output/domainTaylor' sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx') sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = approxRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = approxPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['E'], onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['E'], save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['E'], save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) diff --git a/examples/pod/scatteringSquare.py b/examples/pod/scatteringSquare.py index 5832d25..657df94 100644 --- a/examples/pod/scatteringSquare.py +++ b/examples/pod/scatteringSquare.py @@ -1,172 +1,172 @@ from copy import copy import numpy as np -from rrompy.hfengines.scipy import ( - HelmholtzCavityScatteringProblemEngine as CSPE) +from rrompy.hfengines.linear_problem import \ + HelmholtzCavityScatteringProblemEngine as CSPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from operator import itemgetter def subdict(d, ks): return dict(zip(ks, itemgetter(*ks)(d))) verb = 0 #################### test = "solve" test = "Taylor" test = "Lagrange" test = "TaylorSweep" #test = "LagrangeSweep" plotSamples = True k0 = 10 kLeft, kRight = 9, 11 ktar = 9.5 ktars = np.linspace(8.5, 11.5, 125) #ktars = np.array([k0]) kappa = 5 n = 50 solver = CSPE(kappa = kappa, n = n, verbosity = verb) solver.omega = k0 if test == "solve": uh = solver.solve(k0) print(solver.norm(uh)) solver.plot(uh, what = ['ABS', 'REAL']) elif test in ["Taylor", "Lagrange"]: if test == "Taylor": params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Krylov', 'POD':True} params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD']) parRB = subdict(params, ['R', 'E', 'sampleType', 'POD']) approxPade = TP(solver, mu0 = k0, approxParameters = parPade, verbosity = verb) approxRB = TRB(solver, mu0 = k0, approxParameters = parRB, verbosity = verb) else: params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"MONOMIAL", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"CHEBYSHEV", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler', 'basis']) parRB = subdict(params, ['R', 'S', 'POD', 'sampler']) approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parPade, verbosity = verb) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parRB, verbosity = verb) approxPade.setupApprox() approxRB.setupApprox() if plotSamples: approxPade.plotSamples() PadeErr, solNorm = approxPade.normErr(ktar), approxPade.normHF(ktar) RBErr = approxRB.normErr(ktar) print(('SolNorm:\t{}\nErrPade:\t{}\nErrRelPade:\t{}\nErrRB:\t\t{}' '\nErrRelRB:\t{}').format(solNorm, PadeErr, np.divide(PadeErr, solNorm), RBErr, np.divide(RBErr, solNorm))) print('\nPoles Pade'':') print(approxPade.getPoles()) print('\nPoles RB:') print(approxRB.getPoles()) approxPade.plotHF(ktar, name = 'u_ex') approxPade.plotApp(ktar, name = 'u_Pade''') approxRB.plotApp(ktar, name = 'u_RB') approxPade.plotErr(ktar, name = 'errPade''') approxRB.plotErr(ktar, name = 'errRB') elif test in ["TaylorSweep", "LagrangeSweep"]: if test == "TaylorSweep": shift = 5 nsets = 4 stride = 3 Emax = stride * (nsets - 1) + shift + 1 params = {'Emax':Emax, 'sampleType':'Krylov', 'POD':True} params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift,#+1, 'E': stride*i+shift+1} paramsSetsRB[i] = {'R': stride*i+shift+1,#+1, 'E': stride*i+shift+1} approxPade = TP(solver, mu0 = k0,approxParameters = params, verbosity = verb) approxRB = TRB(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: shift = 7 nsets = 4 stride = 3 Smax = stride * (nsets - 1) + shift + 2 paramsPade = {'S':Smax, 'POD':True, 'basis':"MONOMIAL", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} paramsPade = {'S':Smax, 'POD':True, 'basis':"CHEBYSHEV", 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} paramsRB = copy(paramsPade) paramsPoly = copy(paramsPade) paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift+1, 'S': stride*i+shift+2} paramsSetsRB[i] = {'R': stride*i+shift+2, 'S': stride*i+shift+2} paramsSetsPoly[i] = {'N': 0, 'M': stride*i+shift+1, 'S': stride*i+shift+2} approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsPade, verbosity = verb) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsRB, verbosity = verb) approxPoly = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsPoly, verbosity = verb) filenamebase = '../data/output/scatSquare' + test[:-5] sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = approxRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = approxPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') if test == "TaylorSweep": constr = ['E'] else: constr = ['S'] sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], constr, onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], constr, save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], constr, save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) diff --git a/rrompy/hfengines/base/__init__.py b/rrompy/hfengines/base/__init__.py index de93c79..00f3c2e 100644 --- a/rrompy/hfengines/base/__init__.py +++ b/rrompy/hfengines/base/__init__.py @@ -1,34 +1,30 @@ # 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 .problem_engine_base import ProblemEngineBase +from .mixed_problem_engine_base import MixedProblemEngineBase from .boundary_conditions import BoundaryConditions -from .laplace_base_problem_engine import LaplaceBaseProblemEngine -from .helmholtz_problem_engine import HelmholtzProblemEngine -from .scattering_problem_engine import ScatteringProblemEngine __all__ = [ 'ProblemEngineBase', - 'BoundaryConditions', - 'LaplaceBaseProblemEngine', - 'HelmholtzProblemEngine', - 'ScatteringProblemEngine' + 'MixedProblemEngineBase', + 'BoundaryConditions' ] diff --git a/rrompy/hfengines/base/mixed_problem_engine_base.py b/rrompy/hfengines/base/mixed_problem_engine_base.py new file mode 100644 index 0000000..e09384c --- /dev/null +++ b/rrompy/hfengines/base/mixed_problem_engine_base.py @@ -0,0 +1,154 @@ +# 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 scipy.sparse import csr_matrix +from matplotlib import pyplot as plt +from rrompy.utilities.base.types import Np1D, strLst +from .problem_engine_base import ProblemEngineBase +from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth + +__all__ = ['MixedProblemEngineBase'] + +class MixedProblemEngineBase(ProblemEngineBase): + """ + Generic solver for parametric mixed problems. + + Attributes: + verbosity: Verbosity level. + BCManager: Boundary condition manager. + V: Real mixed FE space. + u: Generic mixed trial functions for variational form evaluation. + v: Generic mixed 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. + bsmu: Mu value of last bs evaluation. + liftDirichletDatamu: Mu value of last Dirichlet datum evaluation. + liftedDirichletDatum: Dofs of Dirichlet datum lifting. + mu0BC: Mu value of last Dirichlet datum lifting. + degree_threshold: Threshold for ufl expression interpolation degree. + """ + + nAs, nbs = 1, 1 + functional = lambda self, u: 0. + + def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10): + super().__init__(degree_threshold = degree_threshold, + verbosity = verbosity) + V1 = fen.FiniteElement("P", fen.triangle, 1) + V2 = fen.FiniteElement("R", fen.triangle, 0) + element = fen.MixedElement([V1, V2]) + self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), element) + self.primalIndices = [0] + + def buildEnergyNormForm(self): # L2 primal + """ + Build sparse matrix (in CSR format) representative of scalar product. + """ + if self.verbosity >= 20: + verbosityDepth("INIT", "Assembling energy matrix.", end = "") + a = 0. + for j in self.primalIndices: + a = a + fen.dot(self.u[j], self.v[j]) + normMatFen = fen.assemble(a * fen.dx) + normMat = fen.as_backend_type(normMatFen).mat() + self.energyNormMatrix = csr_matrix(normMat.getValuesCSR()[::-1], + shape = normMat.size) + if self.verbosity >= 20: + verbosityDepth("DEL", " Done.", inline = True) + + def plot(self, u:Np1D, name : strLst = "u", save : str = None, + what : strLst = 'all', saveFormat : str = "eps", + saveDPI : int = 100, **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'. + save(optional): Where to save plot(s). Defaults to None, i.e. no + saving. + what(optional): Which plots to do. If list, can contain 'ABS', + 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. + Defaults to 'ALL'. + saveFormat(optional): Format for saved plot(s). Defaults to "eps". + saveDPI(optional): DPI for saved plot(s). Defaults to 100. + 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 = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], + listname = self.name() + ".what", baselevel = 1) + if len(what) == 0: return + if 'figsize' not in figspecs.keys(): + figspecs['figsize'] = (13. * len(what) / 4, 3) + if isinstance(name, (str,)): + name = ["{}_sub{}".format(name, j + 1)\ + for j in range(self.V.num_sub_spaces())] + + for j in range(self.V.num_sub_spaces()): + subplotcode = 100 + len(what) * 10 + Vj = self.V.sub(j) + dofs = Vj.dofmap().dofs() + uj = u[dofs] + plt.figure(**figspecs) + Vj = Vj.collapse() + plt.jet() + if 'ABS' in what: + uAb = fen.Function(Vj) + uAb.vector()[:] = np.array(np.abs(uj), dtype = float) + subplotcode = subplotcode + 1 + plt.subplot(subplotcode) + p = fen.plot(uAb, title = "|{}|".format(name[j])) + plt.colorbar(p) + if 'PHASE' in what: + uPh = fen.Function(Vj) + uPh.vector()[:] = np.array(np.angle(uj), dtype = float) + subplotcode = subplotcode + 1 + plt.subplot(subplotcode) + p = fen.plot(uPh, title = "phase({})".format(name[j])) + plt.colorbar(p) + if 'REAL' in what: + uRe = fen.Function(Vj) + uRe.vector()[:] = np.array(np.real(uj), dtype = float) + subplotcode = subplotcode + 1 + plt.subplot(subplotcode) + p = fen.plot(uRe, title = "Re({})".format(name[j])) + plt.colorbar(p) + if 'IMAG' in what: + uIm = fen.Function(Vj) + uIm.vector()[:] = np.array(np.imag(uj), dtype = float) + subplotcode = subplotcode + 1 + plt.subplot(subplotcode) + p = fen.plot(uIm, title = "Im({})".format(name[j])) + plt.colorbar(p) + if save is not None: + save = save.strip() + plt.savefig(getNewFilename("{}_sub{}_fig_".format(save, j + 1), + saveFormat), + format = saveFormat, dpi = saveDPI) + plt.show() + plt.close() diff --git a/rrompy/hfengines/base/problem_engine_base.py b/rrompy/hfengines/base/problem_engine_base.py index 5d01f7d..452d221 100644 --- a/rrompy/hfengines/base/problem_engine_base.py +++ b/rrompy/hfengines/base/problem_engine_base.py @@ -1,366 +1,368 @@ # 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 fenics as fen import numpy as np from scipy.sparse import csr_matrix import scipy.sparse as scsp import scipy.sparse.linalg as scspla from matplotlib import pyplot as plt -from rrompy.utilities.base.types import (Np1D, ScOp, strLst, FenFunc, +from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, strLst, FenFunc, Tuple, List) from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth from .boundary_conditions import BoundaryConditions __all__ = ['ProblemEngineBase'] class ProblemEngineBase: """ Generic solver for parametric problems. 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. energyNormMatrix: Scipy sparse matrix representing inner product over V. bsmu: Mu value of last bs evaluation. liftDirichletDatamu: Mu value of last Dirichlet datum evaluation. liftedDirichletDatum: Dofs of Dirichlet datum lifting. mu0BC: Mu value of last Dirichlet datum lifting. degree_threshold: Threshold for ufl expression interpolation degree. """ nAs, nbs = 1, 1 functional = lambda self, u: 0. def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10): self.BCManager = BoundaryConditions("Dirichlet") self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1) self.verbosity = verbosity self.resetAs() self.resetbs() self.bsmu = np.nan self.liftDirichletDatamu = np.nan self.mu0BC = np.nan self.degree_threshold = degree_threshold def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] @property def V(self): """Value of V.""" return self._V @V.setter def V(self, V): self.resetAs() self.resetbs() 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) - def innerProduct(self, u:Np1D, v:Np1D) -> float: + def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D: """Hilbert space scalar product.""" if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() + if onlyDiag: + return np.sum(self.energyNormMatrix.dot(u) * v.conj(), axis = 0) return v.conj().T.dot(self.energyNormMatrix.dot(u)) - def buildEnergyNormForm(self): + def buildEnergyNormForm(self): # L2 """ Build sparse matrix (in CSR format) representative of scalar product. """ if self.verbosity >= 20: verbosityDepth("INIT", "Assembling energy matrix.", end = "") normMatFen = fen.assemble(fen.dot(self.u, self.v) * fen.dx) normMat = fen.as_backend_type(normMatFen).mat() self.energyNormMatrix = csr_matrix(normMat.getValuesCSR()[::-1], shape = normMat.size) if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) - def norm(self, u:Np1D) -> float: - return np.abs(self.innerProduct(u, u)) ** .5 + def norm(self, u:Np2D) -> Np1D: + return np.abs(self.innerProduct(u, u, onlyDiag = True)) ** .5 def rescaling(self, x:Np1D) -> Np1D: """Rescaling in parameter dependence.""" return x def rescalingInv(self, x:Np1D) -> Np1D: """Inverse rescaling in parameter dependence.""" return x def checkAInBounds(self, der : int = 0): """Check if derivative index is oob for operator of linear system.""" if der < 0 or der >= self.nAs: d = self.V.dim() return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)), shape = (d, d), dtype = np.complex) def checkbInBounds(self, der : int = 0, homogeneized : bool = False): """Check if derivative index is oob for RHS of linear system.""" if der < 0 or der >= max(self.nbs, self.nAs * homogeneized): return np.zeros(self.V.dim(), dtype = np.complex) def setDirichletDatum(self, mu:complex): """Set Dirichlet datum if parametric.""" if hasattr(self, "liftedDirichletDatum"): self.liftDirichletDatamu = mu def liftDirichletData(self, mu:complex) -> Np1D: """Lift Dirichlet datum.""" self.setDirichletDatum(mu) if not np.isclose(self.liftDirichletDatamu, mu): try: liftRe = fen.interpolate(self.DirichletDatum[0], self.V) except: liftRe = fen.project(self.DirichletDatum[0], self.V) try: liftIm = fen.interpolate(self.DirichletDatum[1], self.V) except: liftIm = fen.project(self.DirichletDatum[1], self.V) self.liftedDirichletDatum = (np.array(liftRe.vector()) + 1.j * np.array(liftIm.vector())) return self.liftedDirichletDatum def resetAs(self): """Reset (derivatives of) operator of linear system.""" self.As = [None] * self.nAs def resetbs(self): """Reset (derivatives of) RHS of linear system.""" self.bs = {True: [None] * max(self.nbs, self.nAs), False: [None] * self.nbs} def reduceQuadratureDegree(self, fun:FenFunc, name:str): """Check whether to reduce compiler parameters to degree threshold.""" if not np.isinf(self.degree_threshold): from ufl.algorithms.estimate_degrees import ( estimate_total_polynomial_degree as ETPD) try: deg = ETPD(fun) except: return False if deg > self.degree_threshold: if self.verbosity >= 15: verbosityDepth("MAIN", ("Reducing quadrature degree from " "{} to {} for {}.").format( deg, self.degree_threshold, name)) return True return False def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]): """ Iterate reduceQuadratureDegree over list and define reduce compiler parameters. """ if funsNames is not None: for fun, name in funsNames: if self.reduceQuadratureDegree(fun, name): return {"quadrature_degree" : self.degree_threshold} return {} @abstractmethod def A(self, mu:complex, der : int = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" Anull = self.checkAInBounds(der) if Anull is not None: return Anull if self.As[der] is None: self.As[der] = 0. return self.As[der] @abstractmethod def b(self, mu:complex, der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" bnull = self.checkbInBounds(der, homogeneized) if bnull is not None: return bnull if self.bs[homogeneized][der] is None: self.bs[homogeneized][der] = 0. return self.bs[homogeneized][der] def affineBlocksA(self, mu : complex = 0.) -> Tuple[List[Np1D], callable]: """Assemble affine blocks of operator of linear system.""" def lambdas(x, j): if j == 0: return np.ones(np.size(x)) if j in range(1, self.nAs): return np.power(self.rescaling(x) - self.rescaling(mu), j) raise Exception("Wrong j value.") As = [None] * self.nAs for j in range(self.nAs): As[j] = self.A(mu, j) return As, lambdas def affineBlocksb(self, mu : complex = 0., homogeneized : bool = False)\ -> Tuple[List[Np1D], callable]: """Assemble affine blocks of RHS of linear system.""" def lambdas(x, j): if j == 0: return np.ones(np.size(x)) if j in range(1, self.nbsEff): return np.power(self.rescaling(x) - self.rescaling(mu), j) raise Exception("Wrong j value.") if homogeneized: self.nbsEff = max(self.nAs, self.nbs) else: self.nbsEff = self.nbs bs = [None] * self.nbsEff for j in range(self.nbsEff): bs[j] = self.b(mu, j, homogeneized) return bs, lambdas def solve(self, mu:complex, RHS : Np1D = None, homogeneized : bool = False) -> Np1D: """ Find solution of linear system. Args: mu: parameter value. RHS: RHS of linear system. If None, defaults to that of parametric system. Defaults to None. """ A = self.A(mu) if RHS is None: RHS = self.b(mu, 0, homogeneized) return scspla.spsolve(A, RHS) def residual(self, u:Np1D, mu:complex, homogeneized : bool = False) -> Np1D: """ Find residual of linear system for given approximate solution. Args: u: numpy complex array with function dofs. If None, set to 0. mu: parameter value. """ A = self.A(mu) RHS = self.b(mu, 0, homogeneized) if u is None: return RHS return RHS - A.dot(u) def plot(self, u:Np1D, name : str = "u", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, **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'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. 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 = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], listname = self.name() + ".what", baselevel = 1) if len(what) == 0: return if 'figsize' not in figspecs.keys(): figspecs['figsize'] = (13. * len(what) / 4, 3) subplotcode = 100 + len(what) * 10 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)) 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)) plt.colorbar(p) if save is not None: save = save.strip() plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat), format = saveFormat, dpi = saveDPI) plt.show() plt.close() def plotmesh(self, name : str = "Mesh", save : str = None, saveFormat : str = "eps", saveDPI : int = 100, **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): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ plt.figure(**figspecs) fen.plot(self.V.mesh()) if save is not None: save = save.strip() plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat), format = saveFormat, dpi = saveDPI) plt.show() plt.close() diff --git a/rrompy/utilities/base/fenics.py b/rrompy/hfengines/linear_mixed_problem/__init__.py similarity index 74% copy from rrompy/utilities/base/fenics.py copy to rrompy/hfengines/linear_mixed_problem/__init__.py index 955b774..da9a255 100644 --- a/rrompy/utilities/base/fenics.py +++ b/rrompy/hfengines/linear_mixed_problem/__init__.py @@ -1,30 +1,28 @@ # 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 +from .mixed_laplace_base_problem_engine import MixedLaplaceBaseProblemEngine +from .mixed_helmholtz_problem_engine import MixedHelmholtzProblemEngine -__all__ = ['fenZERO','fenONE','bdrTrue','bdrFalse'] +__all__ = [ + 'MixedLaplaceBaseProblemEngine', + 'MixedHelmholtzProblemEngine' + ] -# 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/hfengines/base/helmholtz_problem_engine.py b/rrompy/hfengines/linear_mixed_problem/mixed_helmholtz_problem_engine.py similarity index 77% copy from rrompy/hfengines/base/helmholtz_problem_engine.py copy to rrompy/hfengines/linear_mixed_problem/mixed_helmholtz_problem_engine.py index 9775a24..52ef35c 100644 --- a/rrompy/hfengines/base/helmholtz_problem_engine.py +++ b/rrompy/hfengines/linear_mixed_problem/mixed_helmholtz_problem_engine.py @@ -1,164 +1,167 @@ # 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 scipy.sparse import csr_matrix import fenics as fen -from .laplace_base_problem_engine import LaplaceBaseProblemEngine +from .mixed_laplace_base_problem_engine import MixedLaplaceBaseProblemEngine from rrompy.utilities.base.types import Np1D, ScOp -from rrompy.utilities.base.fenics import fenZERO, fenONE +from rrompy.utilities.base.fenics import fenZERO, fenZEROS, fenONE from rrompy.utilities.base import verbosityDepth -__all__ = ['HelmholtzProblemEngine'] +__all__ = ['MixedHelmholtzProblemEngine'] -class HelmholtzProblemEngine(LaplaceBaseProblemEngine): +class MixedHelmholtzProblemEngine(MixedLaplaceBaseProblemEngine): """ - Solver for generic Helmholtz problems with parametric wavenumber. + Solver for generic mixed Helmholtz problems with parametric wavenumber. - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega + \int_Omega u = 0 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. + V: Real mixed FE space. + u: Generic mixed trial functions for variational form evaluation. + v: Generic mixed 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. bsmu: Mu value of last bs evaluation. liftDirichletDatamu: Mu value of last Dirichlet datum evaluation. 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. refractionIndex: Value of n. 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. A1: Scipy sparse array representation (in CSC format) of A1. b0: Numpy array representation of b0. dsToBeSet: Whether ds needs to be set. """ nAs = 2 def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = 1. self.refractionIndex = fenONE @property def refractionIndex(self): """Value of n.""" return self._refractionIndex @refractionIndex.setter def refractionIndex(self, refractionIndex): self.resetAs() if not isinstance(refractionIndex, (list, tuple,)): refractionIndex = [refractionIndex, fenZERO] self._refractionIndex = refractionIndex def rescaling(self, x:Np1D) -> Np1D: """Rescaling in parameter dependence.""" return np.power(x, 2.) def rescalingInv(self, x:Np1D) -> Np1D: """Inverse rescaling in parameter dependence.""" return np.power(x, .5) def A(self, mu:complex, der : int = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" Anull = self.checkAInBounds(der) if Anull is not None: return Anull self.autoSetDS() if der <= 0 and self.As[0] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A0.") - DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2), 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)) + a0Re = ((aRe * fen.dot(fen.grad(self.u[0]), fen.grad(self.v[0])) + + self.u[0] * self.v[1] + self.u[1] * self.v[0]) * fen.dx + + hRe * fen.dot(self.u[0], self.v[0]) * self.ds(1)) + a0Im = (aIm * fen.dot(fen.grad(self.u[0]), fen.grad(self.v[0])) \ + * fen.dx + + hIm * fen.dot(self.u[0], self.v[0]) * self.ds(1)) A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe) A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm) 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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), - shape = A0ReMat.size) - + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), - shape = A0ImMat.size)) + self.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer), + shape = A0ReMat.size) + + 1.j * csr_matrix((A0Imv, A0Imc, A0Imr), + shape = A0ImMat.size)) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") if der <= 1 and self.As[1] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A1.") - DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2), 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 + a1Re = - n2Re * fen.dot(self.u[0], self.v[0]) * fen.dx + a1Im = - n2Im * fen.dot(self.u[0], self.v[0]) * fen.dx A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe) A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm) 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.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer), - shape = A1ReMat.size) - + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr), - shape = A1ImMat.size)) + self.As[1] = (csr_matrix((A1Rev, A1Rec, A1Rer), + shape = A1ReMat.size) + + 1.j * csr_matrix((A1Imv, A1Imc, A1Imr), + shape = A1ImMat.size)) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") if der == 0: return self.As[0] + mu**2 * self.As[1] return self.As[1] diff --git a/rrompy/hfengines/base/laplace_base_problem_engine.py b/rrompy/hfengines/linear_mixed_problem/mixed_laplace_base_problem_engine.py similarity index 82% copy from rrompy/hfengines/base/laplace_base_problem_engine.py copy to rrompy/hfengines/linear_mixed_problem/mixed_laplace_base_problem_engine.py index dab145b..f9b5b33 100644 --- a/rrompy/hfengines/base/laplace_base_problem_engine.py +++ b/rrompy/hfengines/linear_mixed_problem/mixed_laplace_base_problem_engine.py @@ -1,320 +1,327 @@ # 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 scipy.sparse import csr_matrix import fenics as fen -from .problem_engine_base import ProblemEngineBase +from rrompy.hfengines.base.mixed_problem_engine_base import \ + MixedProblemEngineBase from rrompy.utilities.base.types import Np1D, ScOp -from rrompy.utilities.base.fenics import fenZERO, fenONE +from rrompy.utilities.base.fenics import fenZERO, fenZEROS, fenONE from rrompy.utilities.base import verbosityDepth -__all__ = ['LaplaceBaseProblemEngine'] +__all__ = ['MixedLaplaceBaseProblemEngine'] -class LaplaceBaseProblemEngine(ProblemEngineBase): +class MixedLaplaceBaseProblemEngine(MixedProblemEngineBase): """ - Solver for generic Laplace problems. + Solver for generic mixed Laplace problems. - \nabla \cdot (a \nabla u) = f in \Omega + \int_Omega u = 0 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. + V: Real mixed FE space. + u: Generic mixed trial functions for variational form evaluation. + v: Generic mixed 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. bsmu: Mu value of last bs evaluation. liftDirichletDatamu: Mu value of last Dirichlet datum evaluation. 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. """ def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = 0. self.diffusivity = fenONE self.forcingTerm = fenZERO - self.DirichletDatum = fenZERO + self.DirichletDatum = fenZEROS(2) 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) + MixedProblemEngineBase.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] + DirichletDatum = [DirichletDatum, fenZEROS(2)] 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: if self.verbosity >= 20: verbosityDepth("INIT", "Initializing boundary measures.", end = "") 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 if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) - def buildEnergyNormForm(self): + def buildEnergyNormForm(self): # H1_omega """ Build sparse matrix (in CSR format) representative of scalar product. """ if self.verbosity >= 20: verbosityDepth("INIT", "Assembling energy matrix.", end = "") - 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) + a = (fen.dot(fen.grad(self.u[0]), fen.grad(self.v[0])) + + np.abs(self.omega)**2 * fen.dot(self.u[0], self.v[0])) + normMatFen = fen.assemble(a * fen.dx) normMat = fen.as_backend_type(normMatFen).mat() - self.energyNormMatrix = scsp.csr_matrix(normMat.getValuesCSR()[::-1], - shape = normMat.size) + self.energyNormMatrix = csr_matrix(normMat.getValuesCSR()[::-1], + shape = normMat.size) if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) def A(self, mu:complex, der : int = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" Anull = self.checkAInBounds(der) if Anull is not None: return Anull self.autoSetDS() if self.As[0] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A0.") - DirichletBC0 = fen.DirichletBC(self.V, fenZERO, + DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2), 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)) + a0Re = ((aRe * fen.dot(fen.grad(self.u[0]), fen.grad(self.v[0])) + + self.u[0] * self.v[1] + self.u[1] * self.v[0]) * fen.dx + + hRe * fen.dot(self.u[0], self.v[0]) * self.ds(1)) + a0Im = (aIm * fen.dot(fen.grad(self.u[0]), fen.grad(self.v[0])) \ + * fen.dx + + hIm * fen.dot(self.u[0], self.v[0]) * self.ds(1)) A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe) A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm) 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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), - shape = A0ReMat.size) - + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), - shape = A0ImMat.size)) + self.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer), + shape = A0ReMat.size) + + 1.j * csr_matrix((A0Imv, A0Imc, A0Imr), + shape = A0ImMat.size)) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") return self.As[0] def b(self, mu:complex, der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" bnull = self.checkbInBounds(der, homogeneized) if bnull is not None: return bnull if homogeneized and not np.isclose(self.mu0BC, mu): self.u0BC = self.liftDirichletData(mu) if (max(self.nbs, self.nAs * homogeneized) > 1 and not np.isclose(self.bsmu, mu)): self.bsmu = mu self.resetbs() if self.bs[homogeneized][der] is None: self.autoSetDS() if self.verbosity >= 20: verbosityDepth("INIT", "Assembling forcing term b{}.".format( der)) if der < self.nbs: fRe, fIm = self.forcingTerm g1Re, g1Im = self.NeumannDatum g2Re, g2Im = self.RobinDatumG else: fRe, fIm = fenZERO, fenZERO g1Re, g1Im = fenZERO, fenZERO g2Re, g2Im = fenZERO, fenZERO 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.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)) + L0Re = (fen.dot(fRe, self.v[0]) * fen.dx + + fen.dot(g1Re, self.v[0]) * self.ds(0) + + fen.dot(g2Re, self.v[0]) * self.ds(1)) + L0Im = (fen.dot(fIm, self.v[0]) * fen.dx + + fen.dot(g1Im, self.v[0]) * self.ds(0) + + fen.dot(g2Im, self.v[0]) * self.ds(1)) b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe) b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm) if homogeneized: Ader = self.A(mu, der) b0Re[:] -= np.real(Ader.dot(self.u0BC)) b0Im[:] -= np.imag(Ader.dot(self.u0BC)) - DBCR = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) - DBCI = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) + DBCR = fen.DirichletBC(self.V, fenZEROS(2), + self.DirichletBoundary) + DBCI = fen.DirichletBC(self.V, fenZEROS(2), + self.DirichletBoundary) else: DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0], self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1], self.DirichletBoundary) DBCR.apply(b0Re) DBCI.apply(b0Im) self.bs[homogeneized][der] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling forcing term.") return self.bs[homogeneized][der] diff --git a/rrompy/hfengines/scipy/__init__.py b/rrompy/hfengines/linear_problem/__init__.py similarity index 85% rename from rrompy/hfengines/scipy/__init__.py rename to rrompy/hfengines/linear_problem/__init__.py index 054ddc1..e5a035c 100644 --- a/rrompy/hfengines/scipy/__init__.py +++ b/rrompy/hfengines/linear_problem/__init__.py @@ -1,41 +1,47 @@ # 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 .laplace_base_problem_engine import LaplaceBaseProblemEngine +from .helmholtz_problem_engine import HelmholtzProblemEngine +from .scattering_problem_engine import ScatteringProblemEngine from .helmholtz_box_scattering_problem_engine import \ HelmholtzBoxScatteringProblemEngine from .helmholtz_cavity_scattering_problem_engine import \ HelmholtzCavityScatteringProblemEngine from .helmholtz_square_bubble_problem_engine import \ HelmholtzSquareBubbleProblemEngine from .helmholtz_square_bubble_domain_problem_engine import \ HelmholtzSquareBubbleDomainProblemEngine from .helmholtz_square_transmission_problem_engine import \ HelmholtzSquareTransmissionProblemEngine from .laplace_disk_gaussian import LaplaceDiskGaussian __all__ = [ + 'LaplaceBaseProblemEngine', + 'HelmholtzProblemEngine', + 'ScatteringProblemEngine', 'HelmholtzBoxScatteringProblemEngine', 'HelmholtzCavityScatteringProblemEngine', 'HelmholtzSquareBubbleProblemEngine', 'HelmholtzSquareBubbleDomainProblemEngine', 'HelmholtzSquareTransmissionProblemEngine', 'LaplaceDiskGaussian' ] diff --git a/rrompy/hfengines/scipy/helmholtz_box_scattering_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py similarity index 94% rename from rrompy/hfengines/scipy/helmholtz_box_scattering_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py index 3c249c8..eedd5a4 100644 --- a/rrompy/hfengines/scipy/helmholtz_box_scattering_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_box_scattering_problem_engine.py @@ -1,58 +1,57 @@ # 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.scattering_problem_engine import \ - ScatteringProblemEngine +from .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, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = kappa 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 = "REST" c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] u0R = - fen.cos(kappa * (c * x + s * y)) u0I = - fen.sin(kappa * (c * x + s * y)) self.DirichletDatum = [u0R, u0I] diff --git a/rrompy/hfengines/scipy/helmholtz_cavity_scattering_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py similarity index 94% rename from rrompy/hfengines/scipy/helmholtz_cavity_scattering_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py index d3d11c5..ba217ab 100644 --- a/rrompy/hfengines/scipy/helmholtz_cavity_scattering_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_cavity_scattering_problem_engine.py @@ -1,60 +1,59 @@ # 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.scattering_problem_engine import \ - ScatteringProblemEngine +from .scattering_problem_engine import ScatteringProblemEngine __all__ = ['HelmholtzCavityScatteringProblemEngine'] class HelmholtzCavityScatteringProblemEngine(ScatteringProblemEngine): """ Solver for scattering problem inside a cavity 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, kappa:float, n:int, gamma : float = 0., signR : int = -1, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.signR = signR self.omega = kappa pi = np.pi mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi), n, n) self.V = fen.FunctionSpace(mesh, "P", 3) self.RobinBoundary = (lambda x, on_boundary: on_boundary and np.isclose(x[1], np.pi)) self.DirichletBoundary = "REST" x, y = fen.SpatialCoordinate(mesh)[:] C = 4. / pi ** 4. bR = C * (2 * (x * (pi - x) + y * (2 * pi - y)) + (kappa * gamma) ** 2. * x * (pi - x) * y * (2 * pi - y)) bI = C * signR * 2 * kappa * (gamma * (pi - 2 * x) * y * (pi - y) + 2 * x * (pi - x) * (pi - y)) wR = fen.cos(kappa * signR * (gamma * x + y)) wI = fen.sin(kappa * signR * (gamma * x + y)) self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI] diff --git a/rrompy/hfengines/base/helmholtz_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py similarity index 100% rename from rrompy/hfengines/base/helmholtz_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_problem_engine.py diff --git a/rrompy/hfengines/scipy/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py similarity index 98% rename from rrompy/hfengines/scipy/helmholtz_square_bubble_domain_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py index d5dd759..c63ac70 100644 --- a/rrompy/hfengines/scipy/helmholtz_square_bubble_domain_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py @@ -1,243 +1,242 @@ # 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.utilities.base.types import Np1D, ScOp, Tuple, FenExpr from rrompy.utilities.base.fenics import fenZERO -from rrompy.hfengines.base.helmholtz_problem_engine import \ - HelmholtzProblemEngine +from .helmholtz_problem_engine import HelmholtzProblemEngine from rrompy.utilities.base import verbosityDepth __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. """ nAs, nbs = 3, 20 def __init__(self, kappa:float, theta:float, n:int, mu0 : np.complex = 1., degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = kappa self.kappa = kappa self.theta = theta self.mu0 = mu0 self.forcingTermMu = np.nan mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n) self.V = fen.FunctionSpace(mesh, "P", 3) - def buildEnergyNormForm(self): + def buildEnergyNormForm(self): # H1 """ Build sparse matrix (in CSR format) representative of scalar product. """ if self.verbosity >= 20: verbosityDepth("INIT", "Assembling energy matrix.", end = "") mudx = np.abs(self.mu0) * fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx muM = np.abs(self.mu0) * fen.dot(self.u, self.v) * fen.dx imudy = 1. / np.abs(self.mu0) * (fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx) normMatFen = fen.assemble(mudx + imudy + muM) normMat = fen.as_backend_type(normMatFen).mat() self.energyNormMatrix = scsp.csr_matrix(normMat.getValuesCSR()[::-1], shape = normMat.size) if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) def getForcingTerm(self, mu:complex) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" if not np.isclose(mu, self.forcingTermMu): if self.verbosity >= 25: verbosityDepth("INIT", "Assembling base expression for forcing term.", end = "") 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), np.imag(mu) mu2R, mu2I = np.real(mu ** 2.), np.imag(mu ** 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 if self.verbosity >= 25: verbosityDepth("DEL", " Done.", inline = True) return self.forcingTerm def getExtraFactorB(self, mu:complex, der:int) -> Tuple[FenExpr, FenExpr]: """Compute extra expression in RHS.""" 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 if self.verbosity >= 25: verbosityDepth("INIT", ("Assembling auxiliary expression for " "forcing term derivative."), end = "") from math import factorial as fact y = fen.SpatialCoordinate(self.V.mesh())[1] powR, powI = [(self.kappa * np.sin(self.theta)) ** der * k\ for k in getPowMinusj(y, der)] mu2R, mu2I = np.real(mu ** 2.), np.imag(mu ** 2.) exprR = mu2R * powR - mu2I * powI exprI = mu2I * powR + mu2R * powI if der >= 1: muR, muI = np.real(2. * mu), np.imag(2. * mu) powR, powI = [(self.kappa * np.sin(self.theta)) ** (der - 1) * k\ * der for k in getPowMinusj(y, der - 1)] exprR += muR * powR - muI * powI exprI += muI * powR + muR * powI if der >= 2: powR, powI = [(self.kappa * np.sin(self.theta)) ** (der - 2) * k\ * der * (der - 1) for k in getPowMinusj(y, der - 2)] exprR += powR exprI += powI fac = fact(der) if self.verbosity >= 25: verbosityDepth("DEL", " Done.", inline = True) return [exprR / fac, exprI / fac] def rescaling(self, x:Np1D) -> Np1D: """Rescaling in parameter dependence.""" return x def rescalingInv(self, x:Np1D) -> Np1D: """Inverse rescaling in parameter dependence.""" return x def A(self, mu:complex, der : int = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" Anull = self.checkAInBounds(der) if Anull is not None: return Anull self.autoSetDS() if der <= 0 and self.As[0] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A0.") DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx A0Re = fen.assemble(a0Re) DirichletBC0.apply(A0Re) A0ReMat = fen.as_backend_type(A0Re).mat() A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR() self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer), shape = A0ReMat.size, dtype = np.complex) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") if der <= 2 and self.As[2] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A2.") 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 A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe) A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm) 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.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer), shape = A2ReMat.size) + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr), shape = A2ImMat.size)) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") if der == 0: return self.As[0] + mu ** 2 * self.As[2] if der == 1: return 2. * mu * self.As[2] return self.As[2] def b(self, mu:complex, der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" bnull = self.checkbInBounds(der, homogeneized) if bnull is not None: return bnull if homogeneized and not np.isclose(self.mu0BC, mu): self.u0BC = self.liftDirichletData(mu) if not np.isclose(self.bsmu, mu): self.bsmu = mu self.resetbs() if self.bs[homogeneized][der] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling forcing term b{}.".format(der)) if der < self.nbs: fRe, fIm = self.getForcingTerm(mu) cRe, cIm = self.getExtraFactorB(mu, der) cfRe = cRe * fRe - cIm * fIm cfIm = cRe * fIm + cIm * fRe else: cfRe, cfIm = fenZERO, fenZERO parsRe = self.iterReduceQuadratureDegree(zip([cfRe], ["forcingTermDer{}Real".format(der)])) parsIm = self.iterReduceQuadratureDegree(zip([cfIm], ["forcingTermDer{}Imag".format(der)])) L0Re = fen.dot(cfRe, self.v) * fen.dx L0Im = fen.dot(cfIm, self.v) * fen.dx b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe) b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm) if homogeneized: Ader = self.A(mu, der) b0Re[:] -= np.real(Ader.dot(self.u0BC)) b0Im[:] -= np.imag(Ader.dot(self.u0BC)) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) DirichletBC0.apply(b0Re) DirichletBC0.apply(b0Im) self.bs[homogeneized][der] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling forcing term.") return self.bs[homogeneized][der] diff --git a/rrompy/hfengines/scipy/helmholtz_square_bubble_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py similarity index 93% rename from rrompy/hfengines/scipy/helmholtz_square_bubble_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py index 1b8a05c..aab7986 100644 --- a/rrompy/hfengines/scipy/helmholtz_square_bubble_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_problem_engine.py @@ -1,53 +1,52 @@ # 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.helmholtz_problem_engine import \ - HelmholtzProblemEngine +from .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, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = kappa pi = np.pi mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi), n, n) self.V = fen.FunctionSpace(mesh, "P", 3) c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] C = 16. / pi ** 4. bR = C * 2 * (x * (pi - x) + y * (pi - y)) bI = C * 2 * kappa * (c * (pi - 2 * x) * y * (pi - y) + s * x * (pi - x) * (pi - 2 * y)) wR = fen.cos(kappa * (c * x + s * y)) wI = fen.sin(kappa * (c * x + s * y)) self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI] diff --git a/rrompy/hfengines/scipy/helmholtz_square_transmission_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py similarity index 95% rename from rrompy/hfengines/scipy/helmholtz_square_transmission_problem_engine.py rename to rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py index 277c9a2..1d2d9b7 100644 --- a/rrompy/hfengines/scipy/helmholtz_square_transmission_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_square_transmission_problem_engine.py @@ -1,78 +1,77 @@ # 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 import ufl -from rrompy.hfengines.base.helmholtz_problem_engine import \ - HelmholtzProblemEngine +from .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, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = kappa 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) 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 = fen.SpatialCoordinate(mesh)[:] TR, TI = np.real(T), np.imag(T) if np.isclose(np.imag(Ky), 0.): u0RT = (TR * fen.cos(Kx * x + np.real(Ky) * y) - TI * fen.sin(Kx * x + np.real(Ky) * y)) u0IT = (TR * fen.sin(Kx * x + np.real(Ky) * y) + TI * fen.cos(Kx * x + np.real(Ky) * y)) else: u0RT = fen.exp(- np.imag(Ky) * y) * (TR * fen.cos(Kx * x) - TI * fen.sin(Kx * x)) u0IT = fen.exp(- np.imag(Ky) * y) * (TR * fen.sin(Kx * x) + TI * fen.cos(Kx * x)) u0RB = (fen.cos(kappa * nB * (dx * x + dy * y)) + (TR - 1) * fen.cos(kappa * nB * (dx*x - dy*y)) - TI * fen.sin(kappa * nB * (dx*x - dy*y))) u0IB = (fen.sin(kappa * nB * (dx * x + dy * y)) + (TR - 1) * fen.sin(kappa * nB * (dx*x - dy*y)) + TI * fen.cos(kappa * nB * (dx*x - dy*y))) u0R = ufl.conditional(ufl.ge(y, 0.), u0RT, u0RB) u0I = ufl.conditional(ufl.ge(y, 0.), u0IT, u0IB) self.refractionIndex = ufl.conditional(ufl.ge(y, 0.), fen.Constant(nT), fen.Constant(nB)) self.DirichletDatum = [u0R, u0I] diff --git a/rrompy/hfengines/base/laplace_base_problem_engine.py b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py similarity index 99% rename from rrompy/hfengines/base/laplace_base_problem_engine.py rename to rrompy/hfengines/linear_problem/laplace_base_problem_engine.py index dab145b..95f39fa 100644 --- a/rrompy/hfengines/base/laplace_base_problem_engine.py +++ b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py @@ -1,320 +1,320 @@ # 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 .problem_engine_base import ProblemEngineBase +from rrompy.hfengines.base.problem_engine_base import ProblemEngineBase from rrompy.utilities.base.types import Np1D, ScOp from rrompy.utilities.base.fenics import fenZERO, fenONE from rrompy.utilities.base import verbosityDepth __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. energyNormMatrix: Scipy sparse matrix representing inner product over V. bsmu: Mu value of last bs evaluation. liftDirichletDatamu: Mu value of last Dirichlet datum evaluation. 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. """ def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.omega = 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: if self.verbosity >= 20: verbosityDepth("INIT", "Initializing boundary measures.", end = "") 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 if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) - def buildEnergyNormForm(self): + def buildEnergyNormForm(self): # H1_omega """ Build sparse matrix (in CSR format) representative of scalar product. """ if self.verbosity >= 20: verbosityDepth("INIT", "Assembling energy matrix.", end = "") 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() self.energyNormMatrix = scsp.csr_matrix(normMat.getValuesCSR()[::-1], shape = normMat.size) if self.verbosity >= 20: verbosityDepth("DEL", " Done.", inline = True) def A(self, mu:complex, der : int = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" Anull = self.checkAInBounds(der) if Anull is not None: return Anull self.autoSetDS() if self.As[0] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling operator term A0.") 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)) A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe) A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm) 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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer), shape = A0ReMat.size) + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr), shape = A0ImMat.size)) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling operator term.") return self.As[0] def b(self, mu:complex, der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" bnull = self.checkbInBounds(der, homogeneized) if bnull is not None: return bnull if homogeneized and not np.isclose(self.mu0BC, mu): self.u0BC = self.liftDirichletData(mu) if (max(self.nbs, self.nAs * homogeneized) > 1 and not np.isclose(self.bsmu, mu)): self.bsmu = mu self.resetbs() if self.bs[homogeneized][der] is None: self.autoSetDS() if self.verbosity >= 20: verbosityDepth("INIT", "Assembling forcing term b{}.".format( der)) if der < self.nbs: fRe, fIm = self.forcingTerm g1Re, g1Im = self.NeumannDatum g2Re, g2Im = self.RobinDatumG else: fRe, fIm = fenZERO, fenZERO g1Re, g1Im = fenZERO, fenZERO g2Re, g2Im = fenZERO, fenZERO 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.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, form_compiler_parameters = parsRe) b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm) if homogeneized: Ader = self.A(mu, der) b0Re[:] -= np.real(Ader.dot(self.u0BC)) b0Im[:] -= np.imag(Ader.dot(self.u0BC)) DBCR = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) else: DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0], self.DirichletBoundary) DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1], self.DirichletBoundary) DBCR.apply(b0Re) DBCI.apply(b0Im) self.bs[homogeneized][der] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling forcing term.") return self.bs[homogeneized][der] diff --git a/rrompy/hfengines/scipy/laplace_disk_gaussian.py b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py similarity index 97% rename from rrompy/hfengines/scipy/laplace_disk_gaussian.py rename to rrompy/hfengines/linear_problem/laplace_disk_gaussian.py index 245eab1..b6076cb 100644 --- a/rrompy/hfengines/scipy/laplace_disk_gaussian.py +++ b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py @@ -1,150 +1,149 @@ # 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, FenExpr -from rrompy.hfengines.base.laplace_base_problem_engine import ( - LaplaceBaseProblemEngine) +from .laplace_base_problem_engine import LaplaceBaseProblemEngine from rrompy.utilities.base.fenics import fenZERO, fenONE from rrompy.utilities.base import verbosityDepth __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. """ nbs = 20 def __init__(self, n:int, degree_threshold : int = np.inf, verbosity : int = 10): super().__init__(degree_threshold = degree_threshold, verbosity = verbosity) self.computebsFactors() self.forcingTermMu = np.nan import mshr mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), n) self.V = fen.FunctionSpace(mesh, "P", 3) def getForcingTerm(self, mu:complex) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" if not np.isclose(mu, self.forcingTermMu): if self.verbosity >= 25: verbosityDepth("INIT", "Assembling base expression for forcing term.", end = "") x, y = fen.SpatialCoordinate(self.V.mesh())[:] C = np.exp(-.5 * mu ** 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), np.imag(mu) 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 if self.verbosity >= 25: verbosityDepth("DEL", " Done.", inline = True) 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:complex, der:int) -> Tuple[FenExpr, FenExpr]: """Compute extra expression in RHS.""" if self.verbosity >= 25: verbosityDepth("INIT", ("Assembling auxiliary expression for " "forcing term derivative."), end = "") muR, muI = np.real(mu), np.imag(mu) x = fen.SpatialCoordinate(self.V.mesh())[0] l = der % 2 if l == 0: powR, powI = fenONE, fenZERO else: powR, powI = x - muR, fen.Constant(muI) exprR, exprI = [self.bsFactors[der, l] * k for k in [powR, powI]] for j in range(l + 2, der + 1, 2): for _ in range(2): powR, powI = (powR * (x - muR) - powI * muI, powR * muI + powI * (x - muR)) exprR += self.bsFactors[der, j] * powR exprI += self.bsFactors[der, j] * powI if self.verbosity >= 25: verbosityDepth("DEL", " Done.", inline = True) return[exprR, exprI] def b(self, mu:complex, der : int = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" bnull = self.checkbInBounds(der, homogeneized) if bnull is not None: return bnull if homogeneized and not np.isclose(self.mu0BC, mu): self.u0BC = self.liftDirichletData(mu) if not np.isclose(self.bsmu, mu): self.bsmu = mu self.resetbs() if self.bs[homogeneized][der] is None: if self.verbosity >= 20: verbosityDepth("INIT", "Assembling forcing term b{}.".format(der)) if der < self.nbs: fRe, fIm = self.getForcingTerm(mu) cRe, cIm = self.getExtraFactorB(mu, der) cfRe = cRe * fRe - cIm * fIm cfIm = cRe * fIm + cIm * fRe else: cfRe, cfIm = fenZERO, fenZERO parsRe = self.iterReduceQuadratureDegree(zip([cfRe], ["forcingTermDer{}Real".format(der)])) parsIm = self.iterReduceQuadratureDegree(zip([cfIm], ["forcingTermDer{}Imag".format(der)])) L0Re = fen.dot(cfRe, self.v) * fen.dx L0Im = fen.dot(cfIm, self.v) * fen.dx b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe) b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm) if homogeneized: Ader = self.A(mu, der) b0Re[:] -= np.real(Ader.dot(self.u0BC)) b0Im[:] -= np.imag(Ader.dot(self.u0BC)) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) DirichletBC0.apply(b0Re) DirichletBC0.apply(b0Im) self.bs[homogeneized][der] = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex) if self.verbosity >= 20: verbosityDepth("DEL", "Done assembling forcing term.") return self.bs[homogeneized][der] diff --git a/rrompy/hfengines/base/scattering_problem_engine.py b/rrompy/hfengines/linear_problem/scattering_problem_engine.py similarity index 100% rename from rrompy/hfengines/base/scattering_problem_engine.py rename to rrompy/hfengines/linear_problem/scattering_problem_engine.py diff --git a/rrompy/reduction_methods/base/fit_utils.py b/rrompy/reduction_methods/base/fit_utils.py index 9594a40..b593035 100644 --- a/rrompy/reduction_methods/base/fit_utils.py +++ b/rrompy/reduction_methods/base/fit_utils.py @@ -1,52 +1,55 @@ # 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 numpy import pi, polynomial as po from rrompy.utilities.base.types import DictAny from rrompy.utilities.warning_manager import warn __all__ = ['setupFitCallables'] def setupFitCallables(kind:str) -> DictAny: """Setup main callables and name for data fitting.""" kind = kind.upper() if kind == "CHEBYSHEV": val = po.chebyshev.chebval + valder = lambda x, c: po.chebyshev.chebval(x, po.chebyshev.chebder(c)) vander = po.chebyshev.chebvander fitname = "chebfit" roots = po.chebyshev.chebroots domcoeff = lambda n: 2. ** (n - 1) if n > 0 else 1. elif kind == "LEGENDRE": val = po.legendre.legval + valder = lambda x, c: po.legendre.legval(x, po.legendre.legder(c)) vander = po.legendre.legvander fitname = "legfit" roots = po.legendre.legroots from scipy.special import binom domcoeff = lambda n: (2. ** n * (pi * n) ** -.5 if n > 10 else .5 ** n * binom(2 * n, n)) else: if kind != "MONOMIAL": warn("Fitting basis not recognized. Overriding to 'MONOMIAL'.") val = po.polynomial.polyval + valder = lambda x, c: po.polynomial.polyval(x,po.polynomial.polyder(c)) vander = po.polynomial.polyvander fitname = "polyfit" roots = po.polynomial.polyroots domcoeff = lambda n: 1. - return {"val":val, "vander":vander, "fitname":fitname, + return {"val":val, "valder":valder, "vander":vander, "fitname":fitname, "roots":roots, "domcoeff":domcoeff} diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py index 096a0f6..d179887 100644 --- a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py +++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py @@ -1,479 +1,491 @@ # 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 numpy as np from rrompy.reduction_methods.base import (checkRobustTolerance, setupFitCallables) from .generic_approximant_lagrange import GenericApproximantLagrange from rrompy.utilities.base.types import Np1D, List, HFEng, DictAny from rrompy.utilities.base import verbosityDepth, purgeDict, customFit from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantLagrangePade'] class ApproximantLagrangePade(GenericApproximantLagrange): """ ROM Lagrange Pade' interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; defaults to 2; - 'sampler': sample point generator; defaults to uniform sampler on [0, 1]; - 'basis': type of basis for interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'E': coefficient of interpolant to be minimized; defaults to min(S, M + 1); - 'M': degree of Pade' interpolant numerator; defaults to 0; - 'N': degree of Pade' interpolant denominator; defaults to 0; - 'interpRcond': tolerance for interpolation via numpy.polyfit; defaults to None; - 'robustTol': tolerance for robust Pade' denominator management; defaults to 0. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. ws: Array of snapshot weigths. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'basis': type of basis for interpolation; - 'E': coefficient of interpolant to be minimized; - 'M': degree of Pade' interpolant numerator; - 'N': degree of Pade' interpolant denominator; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust Pade' denominator management. extraApproxParameters: List of approxParameters keys in addition to mother class's. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. M: Numerator degree of approximant. N: Denominator degree of approximant. POD: Whether to compute POD of snapshots. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust Pade' denominator management. samplingEngine: Sampling engine. 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. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["basis", "E", "M", "N", "interpRcond", "robustTol"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """ Value of approximant parameters. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["basis", "E", "M", "N", "interpRcond", "robustTol"], True, True, baselevel = 1) if hasattr(self, "M"): Mold = self.M self._M = 0 if hasattr(self, "N"): Nold = self.N self._N = 0 if hasattr(self, "E"): self._E = 0 GenericApproximantLagrange.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "basis" in keyList or not hasattr(self, "_val"): if "basis" in keyList: kind = approxParameters["basis"] else: kind = "MONOMIAL" setupFit = setupFitCallables(kind) for x in setupFit: super().__setattr__("_" + x, setupFit[x]) if "interpRcond" in keyList: self.interpRcond = approxParameters["interpRcond"] elif hasattr(self, "interpRcond"): self.interpRcond = self.interpRcond else: self.interpRcond = None if "robustTol" in keyList: self.robustTol = approxParameters["robustTol"] elif hasattr(self, "robustTol"): self.robustTol = self.robustTol else: self.robustTol = 0 if "M" in keyList: self.M = approxParameters["M"] elif hasattr(self, "M"): self.M = Mold else: self.M = 0 if "N" in keyList: self.N = approxParameters["N"] elif hasattr(self, "N"): self.N = Nold else: self.N = 0 if "E" in keyList: self.E = approxParameters["E"] else: self.E = min(self.S - 1, self.M + 1) @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: warn("Prescribed S is too small. Updating S to M + 1.") 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: warn("Prescribed S is too small. Updating S to N + 1.") self.S = self.N + 1 @property def E(self): """Value of E. Its assignment may change S.""" 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, "S") and self.S < self.E + 1: warn("Prescribed S is too small. Updating S to E + 1.") self.S = self.E + 1 @property def robustTol(self): """Value of tolerance for robust Pade' denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: warn("Overriding prescribed negative robustness tolerance to 0.") robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @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] * 3, {0:"M", 1:"N", 2:"E"} if hasattr(self, "M"): vals[0] = self.M if hasattr(self, "N"): vals[1] = self.N if hasattr(self, "E"): vals[2] = self.E idxmax = np.argmax(vals) if vals[idxmax] + 1 > S: warn("Prescribed S is too small. Updating S to {} + 1."\ .format(label[idxmax])) self.S = vals[idxmax] + 1 else: self._S = S self._approxParameters["S"] = self.S if Sold != self.S: self.resetSamples() def setupApprox(self): """ Compute Pade' interpolant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeScaleFactor() self.computeSnapshots() if self.N > 0: if self.verbosity >= 7: verbosityDepth("INIT", ("Starting computation of " "denominator.")) TE = self._vander(self.radiusPade(self.mus), self.E) while self.N > 0: RHS = np.zeros(self.E + 1) RHS[-1] = 1. fitOut = customFit(TE.T, RHS, w = self.ws, full = True, rcond = self.interpRcond) if self.verbosity >= 5: verbosityDepth("MAIN", ("Fitting {} samples with " "degree {} through {}... " "Conditioning of LS system: " "{:.4e}.").format( self.S, self.E, self._fitname, fitOut[1][2][0] / fitOut[1][2][-1])) if fitOut[1][1] < self.E + 1: Enew = fitOut[1][1] - 1 Nnew = min(self.N, Enew) Mnew = min(self.M, Enew) if Nnew == self.N: strN = "" else: strN = "N from {} to {} and ".format(self.N, Nnew) if Mnew == self.M: strM = "" else: strM = "M from {} to {} and ".format(self.M, Mnew) warn(("Polyfit is poorly conditioned.\nReducing {}{}E " "from {} to {}.").format(strN, strM, self.E, Enew)) newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew} self.approxParameters = newParameters continue G = (TE[:, : self.N + 1].T * fitOut[0]).T if self.POD: if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square " "root of gramian matrix.")) G = self.samplingEngine.RPOD.dot(G) _, s, eV = np.linalg.svd(G, full_matrices = False) ev = s[::-1] eV = eV[::-1, :].conj().T if self.verbosity >= 5: try: condev = s[0] / s[-1] except: condev = np.inf verbosityDepth("MAIN", ("Solved svd problem of " "size {} x {} with " "condition number " "{:.4e}.").format( self.S, self.N + 1, condev)) else: if self.verbosity >= 10: verbosityDepth("INIT", "Building gramian matrix.", end = "") G = self.samplingEngine.samples.dot(G) G2 = self.HFEngine.innerProduct(G, G) if self.verbosity >= 10: verbosityDepth("DEL", "Done building gramian.", inline = True) if self.verbosity >= 7: verbosityDepth("INIT", ("Solving eigenvalue " "problem for gramian " "matrix.")) ev, eV = np.linalg.eigh(G2) if self.verbosity >= 5: try: condev = ev[-1] / ev[0] except: condev = np.inf verbosityDepth("MAIN", ("Solved eigenvalue " "problem of size {} with " "condition number " "{:.4e}.").format( self.N + 1, condev)) if self.verbosity >= 7: verbosityDepth("DEL", ("Done solving eigenvalue " "problem.")) newParameters = checkRobustTolerance(ev, self.E, self.robustTol) if not newParameters: break self.approxParameters = newParameters if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) self.Q = eV[:, 0] if self.verbosity >= 7: verbosityDepth("DEL", "Done computing denominator.") else: self.Q = np.ones(1, dtype = np.complex) if self.verbosity >= 7: verbosityDepth("INIT", "Starting computation of numerator.") self.lastApproxParameters = copy(self.approxParameters) Qevaldiag = np.diag(self.getQVal(self.mus)) while self.M >= 0: fitVander = self._vander(self.radiusPade(self.mus), self.M) fitOut = customFit(fitVander, Qevaldiag, w = self.ws, full = True, rcond = self.interpRcond) if self.verbosity >= 5: verbosityDepth("MAIN", ("Fitting {} samples with degree " "{} through {}... Conditioning of " "LS system: {:.4e}.").format( self.S, self.M, self._fitname, fitOut[1][2][0] / fitOut[1][2][-1])) if fitOut[1][1] == self.M + 1: P = fitOut[0].T break warn(("Polyfit is poorly conditioned. Reducing M from {} to " "{}. Exact snapshot interpolation not guaranteed.")\ .format(self.M, fitOut[1][1] - 1)) self.M = fitOut[1][1] - 1 if self.M <= 0: raise Exception(("Instability in computation of numerator. " "Aborting.")) self.P = np.atleast_2d(P) if self.POD: self.P = self.samplingEngine.RPOD.dot(self.P) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing numerator.") self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") 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 ((self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)) / self.scaleFactor) def getPVal(self, mu:List[complex]): """ Evaluate Pade' numerator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating numerator at mu = {}.".format(mu), end = "") try: len(mu) except: mu = [mu] p = self._val(self.radiusPade(mu), self.P.T) if len(mu) == 1: p = p.flatten() if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return p def getQVal(self, mu:List[complex]): """ Evaluate Pade' denominator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating denominator at mu = {}.".format(mu), end = "") q = self._val(self.radiusPade(mu), self.Q) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return q def evalApproxReduced(self, mu:complex): """ Evaluate Pade' approximant at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Evaluating approximant at mu = {}.".format(mu)) self.uAppReduced = self.getPVal(mu) / self.getQVal(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done evaluating approximant.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.samplingEngine.samples.dot(self.uAppReduced) + def getResidues(self) -> Np1D: + """ + Obtain norm of approximant residues. + + Returns: + Numpy vector of norms of residues. + """ + poles = self.getPoles() + Pvals = self.samplingEngine.samples.dot(self.getPVal(poles)) + Qder = self._valder(self.radiusPade(poles), self.Q) + return Pvals / Qder + def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() return self.HFEngine.rescalingInv( self.scaleFactor * self._roots(self.Q) + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py index 7b783e0..0adc2f2 100644 --- a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py +++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py @@ -1,209 +1,209 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.base.types import DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth __all__ = ['GenericApproximantLagrange'] class GenericApproximantLagrange(GenericApproximant): """ ROM Lagrange interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; defaults to uniform sampler on [0, 1]. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. ws: Array of snapshot weigths. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'S': total number of snapshots current approximant relies upon; - 'sampler': sample point generator. extraApproxParameters: List of approxParameters keys in addition to mother class's. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. POD: Whether to compute POD of snapshots. samplingEngine: Sampling engine. 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. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["S", "sampler"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() def setupSampling(self): """Setup sampling engine.""" if not hasattr(self, "POD"): return if self.POD: - from rrompy.sampling.scipy.sampling_engine_lagrange_pod import ( - SamplingEngineLagrangePOD) + from rrompy.sampling.linear_problem.sampling_engine_lagrange_pod \ + import SamplingEngineLagrangePOD super().setupSampling(SamplingEngineLagrangePOD) else: - from rrompy.sampling.scipy.sampling_engine_lagrange import ( - SamplingEngineLagrange) + from rrompy.sampling.linear_problem.sampling_engine_lagrange \ + import SamplingEngineLagrange super().setupSampling(SamplingEngineLagrange) @property def mus(self): """Value of mus. Its assignment may reset snapshots.""" return self._mus @mus.setter def mus(self, mus): musOld = self.mus if hasattr(self, 'mus') else None self._mus = np.array(mus) _, musCounts = np.unique(self._mus, return_counts = True) if len(np.where(musCounts > 1)[0]) > 0: raise Exception("Repeated sample points not allowed.") if (musOld is None or len(self.mus) != len(musOld) or not np.allclose(self.mus, musOld, 1e-14)): self.resetSamples() self.autoNode = None @property def approxParameters(self): """Value of approximant parameters. Its assignment may change S.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["S", "sampler"], True, True, baselevel = 1) GenericApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "S" in keyList: self.S = approxParameters["S"] elif hasattr(self, "S"): self.S = self.S else: self.S = 2 if "sampler" in keyList: self.sampler = approxParameters["sampler"] elif not hasattr(self, "S"): from rrompy.utilities.parameter_sampling import QuadratureSampler self.sampler = QuadratureSampler([0., 1.], "UNIFORM") del QuadratureSampler @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): if S <= 0: raise ArithmeticError("S must be positive.") if hasattr(self, "S"): Sold = self.S else: Sold = -1 self._S = S self._approxParameters["S"] = self.S if Sold != self.S: self.resetSamples() @property def sampler(self): """Value of sampler.""" return self._sampler @sampler.setter def sampler(self, sampler): if 'generatePoints' not in dir(sampler): raise Exception("Sampler type not recognized.") if hasattr(self, '_sampler'): samplerOld = self.sampler self._sampler = sampler self._approxParameters["sampler"] = self.sampler if not 'samplerOld' in locals() or samplerOld != self.sampler: self.resetSamples() def computeSnapshots(self): """Compute snapshots of solution map.""" if self.samplingEngine.samples is None: if self.verbosity >= 5: verbosityDepth("INIT", "Starting computation of snapshots.") self.mus, self.ws = self.sampler.generatePoints(self.S) self.mus = np.array([x[0] for x in self.mus]) self.samplingEngine.iterSample(self.mus, homogeneized = self.homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Done computing snapshots.") 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.samplingEngine.samples is not None and super().checkComputedApprox()) def normApp(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ if not self.POD or self.homogeneized != homogeneized: return super().normApp(mu, homogeneized) return np.linalg.norm(self.getAppReduced(mu)) def computeScaleFactor(self): """Compute parameter rescaling factor.""" self.scaleFactor = .5 * np.abs(self.HFEngine.rescaling( self.sampler.lims[0][0]) - self.HFEngine.rescaling( self.sampler.lims[1][0])) diff --git a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py index 1596ba8..a155677 100644 --- a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py +++ b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py @@ -1,592 +1,593 @@ # 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 numpy as np from rrompy.reduction_methods.base import (checkRobustTolerance, setupFitCallables) from .generic_approximant_lagrange_greedy import ( GenericApproximantLagrangeGreedy) from rrompy.reduction_methods.lagrange import ApproximantLagrangePade from rrompy.utilities.base.types import DictAny, List, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth, customFit from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantLagrangePadeGreedy'] class ApproximantLagrangePadeGreedy(GenericApproximantLagrangeGreedy, ApproximantLagrangePade): """ ROM greedy Pade' interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'muBounds': list of bounds for parameter values; defaults to [[0], [1]]; - 'basis': type of basis for interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'Delta': difference between M and N in rational approximant; defaults to 0; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'errorEstimatorKind': kind of error estimator; available values include 'EXACT' and 'SIMPLIFIED'; defaults to 'SIMPLIFIED'; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of training points to be exhausted before training set refinement; defaults to 0.2; - 'nTrainingPoints': number of training points; defaults to maxIter / refinementRatio; - 'nTestPoints': number of starting test points; defaults to 1; - 'trainingSetGenerator': training sample points generator; defaults to uniform sampler within muBounds; - 'interpRcond': tolerance for interpolation via numpy.polyfit; defaults to None; - 'robustTol': tolerance for robust Pade' denominator management; defaults to 0. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'muBounds': list of bounds for parameter values; - 'basis': type of basis for interpolation; - 'Delta': difference between M and N in rational approximant; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'errorEstimatorKind': kind of error estimator; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of training points to be exhausted before training set refinement; - 'nTrainingPoints': number of training points; - 'nTestPoints': number of starting test points; - 'trainingSetGenerator': training sample points generator; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust Pade' denominator management. extraApproxParameters: List of approxParameters keys in addition to mother class's. POD: whether to compute POD of snapshots. muBounds: list of bounds for parameter values. greedyTol: uniform error tolerance for greedy algorithm. errorEstimatorKind: kind of error estimator. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTrainingPoints: number of training points. nTestPoints: number of starting test points. trainingSetGenerator: training sample points generator. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust Pade' denominator management. samplingEngine: Sampling engine. 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. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["basis", "Delta", "errorEstimatorKind", "interpRcond", "robustTol"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change robustTol. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["basis", "Delta", "errorEstimatorKind", "interpRcond", "robustTol"], True, True, baselevel = 1) if "Delta" in list(approxParameters.keys()): self._Delta = approxParameters["Delta"] elif hasattr(self, "Delta"): self._Delta = self.Delta else: self._Delta = 0 GenericApproximantLagrangeGreedy.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) self.Delta = self.Delta if "basis" in keyList or not hasattr(self, "_val"): if "basis" in keyList: kind = approxParameters["basis"] else: kind = "MONOMIAL" setupFit = setupFitCallables(kind) for x in setupFit: super().__setattr__("_" + x, setupFit[x]) if "errorEstimatorKind" in keyList: self.errorEstimatorKind = approxParameters["errorEstimatorKind"] elif hasattr(self, "errorEstimatorKind"): self.errorEstimatorKind = self.errorEstimatorKind else: self.errorEstimatorKind = "SIMPLIFIED" if "interpRcond" in keyList: self.interpRcond = approxParameters["interpRcond"] elif hasattr(self, "interpRcond"): self.interpRcond = self.interpRcond else: self.interpRcond = None if "robustTol" in keyList: self.robustTol = approxParameters["robustTol"] elif hasattr(self, "robustTol"): self.robustTol = self.robustTol else: self.robustTol = 0 @property def Delta(self): """Value of Delta.""" return self._Delta @Delta.setter def Delta(self, Delta): if not np.isclose(Delta, np.floor(Delta)): raise ArithmeticError("Delta must be an integer.") if Delta < 0: warn(("Error estimator unreliable for Delta < 0. Overloading of " "errorEstimator is suggested.")) else: Deltamin = (max(self.HFEngine.nbs, self.HFEngine.nAs * self.homogeneized) - 1 - 1 * (self.HFEngine.nAs > 1)) if Delta < Deltamin: warn(("Method may be unreliable for selected Delta. Suggested " "minimal value of Delta: {}.").format(Deltamin)) self._Delta = Delta self._approxParameters["Delta"] = self.Delta @property def errorEstimatorKind(self): """Value of errorEstimatorKind.""" return self._errorEstimatorKind @errorEstimatorKind.setter def errorEstimatorKind(self, errorEstimatorKind): errorEstimatorKind = errorEstimatorKind.upper() if errorEstimatorKind not in ["EXACT", "SIMPLIFIED"]: warn(("Error estimator kind not recognized. Overriding to " "'SIMPLIFIED'.")) errorEstimatorKind = "SIMPLIFIED" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind @property def nTestPoints(self): """Value of nTestPoints.""" return self._nTestPoints @nTestPoints.setter def nTestPoints(self, nTestPoints): if nTestPoints <= np.abs(self.Delta): warn(("nTestPoints must be at least abs(Delta) + 1. Increasing " "value to abs(Delta) + 1.")) nTestPoints = np.abs(self.Delta) + 1 if not np.isclose(nTestPoints, np.int(nTestPoints)): raise ArithmeticError("nTestPoints must be an integer.") nTestPoints = np.int(nTestPoints) if hasattr(self, "nTestPoints"): nTestPointsold = self.nTestPoints else: nTestPointsold = -1 self._nTestPoints = nTestPoints self._approxParameters["nTestPoints"] = self.nTestPoints if nTestPointsold != self.nTestPoints: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() self.resbb = None self.resAb = None self.resAA = None self.As = None self.bs = None def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]: """Standard residual-based error estimator.""" self.setupApprox() self.initEstNormer() PM = self.P[:, -1] if np.any(np.isnan(PM)) or np.any(np.isinf(PM)): err = np.empty(len(mus)) err[:] = np.inf return err nAs = self.HFEngine.nAs - 1 nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized) radiusmus = self.radiusPade(mus) radiussmus = self.radiusPade(self.mus) musTile = np.tile(radiusmus.reshape(-1, 1), [1, self.S]) smusCol = radiussmus.reshape(1, -1) num = np.prod(musTile[:, : self.S] - smusCol, axis = 1) den = self.getQVal(mus) self.assembleReducedResidualBlocks() vanderBase = np.polynomial.polynomial.polyvander(radiusmus, max(nAs, nbs)).T radiusb0 = vanderBase[: nbs + 1, :] # 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj() b0resb0 = np.sum(self.resbb.dot(radiusb0) * radiusb0.conj(), axis = 0) RHSnorms = np.power(np.abs(b0resb0), .5) vanderBase = vanderBase[: -1, :] delta = self.S - self.N - 1 nbsEff = max(0, nbs - delta) if self.errorEstimatorKind == "SIMPLIFIED": radiusA = np.tensordot(PM, vanderBase[: nAs, :], 0) if delta == 0: radiusb = np.abs(self.Q[-1]) * radiusb0[: -1, :] else: #if self.errorEstimatorKind == "EXACT": momentQ = np.zeros(nbsEff, dtype = np.complex) momentQu = np.zeros((self.S, nAs), dtype = np.complex) radiusbTen = np.zeros((nbsEff, nbsEff, len(mus)), dtype = np.complex) radiusATen = np.zeros((nAs, nAs, len(mus)), dtype = np.complex) if nbsEff > 0: momentQ[0] = self.Q[-1] radiusbTen[0, :, :] = vanderBase[: nbsEff, :] momentQu[:, 0] = self.P[:, -1] radiusATen[0, :, :] = vanderBase[: nAs, :] Qvals = self.getQVal(self.mus) for k in range(1, max(nAs, nbs * (nbsEff > 0))): Qvals = Qvals * radiussmus if k > delta and k < nbs: momentQ[k - delta] = self._fitinv.dot(Qvals) radiusbTen[k - delta, k :, :] = ( radiusbTen[0, : delta - k, :]) if k < nAs: momentQu[:, k] = Qvals * self._fitinv radiusATen[k, k :, :] = radiusATen[0, : - k, :] if self.POD and nAs > 1: momentQu[:, 1 :] = self.samplingEngine.RPOD.dot( momentQu[:, 1 :]) radiusA = np.tensordot(momentQu, radiusATen, 1) if nbsEff > 0: radiusb = np.tensordot(momentQ, radiusbTen, 1) if ((self.errorEstimatorKind == "SIMPLIFIED" and delta == 0) or (self.errorEstimatorKind == "EXACT" and nbsEff > 0)): # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj() ff = np.sum(self.resbb[delta + 1 :, delta + 1 :].dot(radiusb) * radiusb.conj(), axis = 0) # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj() Lf = np.sum(np.tensordot(self.resAb[delta :, :, :], radiusA, 2) * radiusb.conj(), axis = 0) else: ff, Lf = 0., 0. # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj() LL = np.sum(np.tensordot(self.resAA, radiusA, 2) * radiusA.conj(), axis = (0, 1)) jOpt = np.power(np.abs(ff - 2. * np.real(Lf) + LL), .5) return self._domcoeff(self.S - 1) * jOpt * np.abs(num / den) / RHSnorms def setupApprox(self): """ Compute Pade' interpolant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeScaleFactor() self.S = len(self.mus) self._M = self.S - 1 self._N = self.S - 1 if self.Delta < 0: self._M += self.Delta else: self._N -= self.Delta if min(self.M, self.N) < 0: if self.verbosity >= 5: verbosityDepth("MAIN", "Minimal sample size not achieved.") self.Q = np.empty(max(self.N, 0) + 1, dtype = np.complex) self.P = np.empty((len(self.mus), max(self.M, 0) + 1), dtype = np.complex) self.Q[:] = np.nan self.P[:] = np.nan self.lastApproxParameters = copy(self.approxParameters) if self.verbosity >= 5: verbosityDepth("DEL", ("Aborting computation of " "approximant.\n")) return self.greedy() if self.N > 0: if self.verbosity >= 7: verbosityDepth("INIT", ("Starting computation of " "denominator.")) TS = self._vander(self.radiusPade(self.mus), self.S - 1) while self.N > 0: RHS = np.zeros(self.S) RHS[-1] = 1. fitOut = customFit(TS.T, RHS, full = True, rcond = self.interpRcond) if self.verbosity >= 2: verbosityDepth("MAIN", ("Fitting {} samples with " "degree {} through {}... " "Conditioning of system: " "{:.4e}.").format(self.S, self.S - 1, self._fitname, fitOut[1][2][0] / fitOut[1][2][-1])) if fitOut[1][1] < self.S: warn(("Polyfit is poorly conditioned. Starting " "preemptive termination of computation of " "approximant.")) self.Q = np.empty(max(self.N, 0) + 1, dtype = np.complex) self.P = np.empty((len(self.mus), max(self.M, 0) + 1), dtype = np.complex) self.Q[:] = np.nan self.P[:] = np.nan self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 7: verbosityDepth("DEL", ("Aborting computation of " "denominator.")) if self.verbosity >= 5: verbosityDepth("DEL", ("Aborting computation of " "approximant.\n")) return self._fitinv = fitOut[0] G = (TS[:, : self.N + 1].T * fitOut[0]).T if self.POD: if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square " "root of gramian matrix.")) G = self.samplingEngine.RPOD.dot(G) _, s, eV = np.linalg.svd(G, full_matrices = False) ev = s[::-1] eV = eV[::-1, :].conj().T if self.verbosity >= 2: try: condev = s[0] / s[-1] except: condev = np.inf verbosityDepth("MAIN", ("Solved svd problem of " "size {} x {} with " "condition number " "{:.4e}.").format( self.S, self.N + 1, condev)) else: if self.verbosity >= 10: verbosityDepth("INIT", "Building gramian matrix.", end = "") G = self.samplingEngine.samples.dot(G) G2 = self.HFEngine.innerProduct(G, G) if self.verbosity >= 10: verbosityDepth("DEL", "Done building gramian.", inline = True) if self.verbosity >= 7: verbosityDepth("INIT", ("Solving eigenvalue " "problem for gramian " "matrix.")) ev, eV = np.linalg.eigh(G2) if self.verbosity >= 2: try: condev = ev[-1] / ev[0] except: condev = np.inf verbosityDepth("MAIN", ("Solved eigenvalue " "problem of size {} with " "condition number " "{:.4e}.").format( self.N + 1, condev)) if self.verbosity >= 7: verbosityDepth("DEL", ("Done solving eigenvalue " "problem.")) newParameters = checkRobustTolerance(ev, self.M, self.robustTol) if not newParameters: break self._N = newParameters["N"] self._M = newParameters["E"] if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) self.Q = eV[:, 0] if self.verbosity >= 7: verbosityDepth("DEL", "Done computing denominator.") else: self.Q = np.ones(1, dtype = np.complex) if self.verbosity >= 7: verbosityDepth("INIT", "Starting computation of numerator.") self.lastApproxParameters = copy(self.approxParameters) Qevaldiag = np.diag(self.getQVal(self.mus)) while self.M >= 0: fitVander = self._vander(self.radiusPade(self.mus), self.M) fitOut = customFit(fitVander, Qevaldiag, full = True, rcond = self.interpRcond) if fitOut[1][1] == self.M + 1: P = fitOut[0].T break warn(("Polyfit is poorly conditioned. Reducing M from {} to " "{}. Exact snapshot interpolation not guaranteed.")\ .format(self.M, fitOut[1][1] - 1)) self._M = fitOut[1][1] - 1 if self.M <= 0: raise Exception(("Instability in computation of numerator. " "Aborting.")) self.P = np.atleast_2d(P) if self.POD: self.P = self.samplingEngine.RPOD.dot(self.P) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing numerator.") self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def assembleReducedResidualBlocks(self): """Build affine blocks of reduced linear system through projections.""" self.initEstNormer() if self.As is None or self.bs is None: if self.verbosity >= 7: verbosityDepth("INIT", "Computing Taylor blocks of system.") nAs = self.HFEngine.nAs - 1 nbs = max(self.HFEngine.nbs, (nAs + 1) * self.homogeneized) self.As = [self.HFEngine.A(self.mu0, j + 1) for j in range(nAs)] self.bs = [self.HFEngine.b(self.mu0, j, self.homogeneized) for j in range(nbs)] if self.verbosity >= 7: verbosityDepth("DEL", "Done computing Taylor blocks.") computeResbb = self.resbb is None computeResAb = self.resAb is None or self.resAb.shape[1] != self.S computeResAA = self.resAA is None or self.resAA.shape[0] != self.S samples = self.samplingEngine.samples if computeResbb or computeResAb or computeResAA: if self.verbosity >= 7: verbosityDepth("INIT", "Projecting Taylor terms of residual.") nAs = len(self.As) nbs = len(self.bs) - 1 if computeResbb: self.resbb = np.empty((nbs + 1, nbs + 1), dtype = np.complex) for i in range(nbs + 1): Mbi = self.scaleFactor ** i * self.bs[i] for j in range(i): Mbj = self.scaleFactor ** j * self.bs[j] self.resbb[i, j] = self.estNormer.innerProduct(Mbj, Mbi) self.resbb[i, i] = self.estNormer.innerProduct(Mbi, Mbi) for i in range(nbs + 1): for j in range(i + 1, nbs + 1): self.resbb[i, j] = self.resbb[j][i].conj() if computeResAb: if self.resAb is None: self.resAb = np.empty((nbs, self.S, nAs), dtype = np.complex) for i in range(nbs): Mbi = self.scaleFactor ** (i + 1) * self.bs[i + 1] for j in range(nAs): MAj = (self.scaleFactor ** (j + 1) * self.As[j].dot(samples)) self.resAb[i, :, j] = self.estNormer.innerProduct( MAj, Mbi) else: Sold = self.resAb.shape[1] if Sold > self.S: self.resAb = self.resAb[:, : self.S, :] else: resAbNew = np.empty((nbs, self.S, nAs), dtype = np.complex) resAbNew[:, : Sold, :] = self.resAb self.resAb = resAbNew for i in range(nbs): Mbi = self.scaleFactor ** (i + 1) * self.bs[i + 1] for j in range(nAs): MAj = (self.scaleFactor ** (j + 1) * self.As[j].dot(samples[:, Sold :])) self.resAb[i, Sold :, j] = ( self.estNormer.innerProduct(MAj, Mbi)) if computeResAA: if self.resAA is None: self.resAA = np.empty((self.S, nAs, self.S, nAs), dtype = np.complex) for i in range(nAs): MAi = (self.scaleFactor ** (i + 1) * self.As[i].dot(samples)) for j in range(i): MAj = (self.scaleFactor ** (j + 1) * self.As[j].dot(samples)) self.resAA[:, i, :, j] = ( self.estNormer.innerProduct(MAj, MAi)) self.resAA[:, i, :, i] = self.estNormer.innerProduct( MAi, MAi) for i in range(nAs): for j in range(i + 1, nAs): self.resAA[:, i, :, j] = ( - self.resAA[j, :, :, i].conj()) + self.resAA[:, j, :, i].conj()) else: Sold = self.resAA.shape[0] if Sold > self.S: self.resAA = self.resAA[: self.S, :, : self.S, :] else: resAANew = np.empty((self.S, nAs, self.S, nAs), dtype = np.complex) resAANew[: Sold, :, : Sold, :] = self.resAA self.resAA = resAANew for i in range(nAs): MAi = (self.scaleFactor ** (i + 1) * self.As[i].dot(samples)) for j in range(i): MAj = (self.scaleFactor ** (j + 1) * self.As[j].dot(samples)) self.resAA[: Sold, i, Sold :, j] = ( self.estNormer.innerProduct(MAj[:, Sold :], MAi[:, : Sold])) self.resAA[Sold :, i, : Sold, j] = ( self.estNormer.innerProduct(MAj[:, : Sold], MAi[:, Sold :])) self.resAA[Sold :, i, Sold :, j] = ( self.estNormer.innerProduct(MAj[:, Sold :], MAi[:, Sold :])) self.resAA[: Sold, i, Sold :, i] = ( self.estNormer.innerProduct(MAi[:, Sold :], MAi[:, : Sold])) self.resAA[Sold :, i, : Sold, i] = ( self.resAA[: Sold, i, Sold :, i].conj().T) self.resAA[Sold :, i, Sold :, i] = ( self.estNormer.innerProduct(MAi[:, Sold :], MAi[:, Sold :])) for i in range(nAs): for j in range(i + 1, nAs): self.resAA[:, i, :, j] = ( self.resAA[:, j, :, i].conj()) if self.verbosity >= 7: verbosityDepth("DEL", ("Done setting up Taylor " "decomposition of residual.")) + diff --git a/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py index fa7936d..81a58b3 100644 --- a/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py +++ b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py @@ -1,489 +1,476 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import ( GenericApproximantLagrange) from rrompy.utilities.base.types import DictAny, HFEng, Tuple, List from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['GenericApproximantLagrangeGreedy'] class GenericApproximantLagrangeGreedy(GenericApproximantLagrange): """ ROM greedy Lagrange interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'muBounds': list of bounds for parameter values; defaults to [[0], [1]]; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of training points to be exhausted before training set refinement; defaults to 0.2; - 'nTrainingPoints': number of training points; defaults to maxIter / refinementRatio; - 'nTestPoints': number of starting test points; defaults to 1; - 'trainingSetGenerator': training sample points generator; defaults to uniform sampler within muBounds; - 'robustTol': tolerance for robust Pade' denominator management; defaults to 0. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'muBounds': list of bounds for parameter values; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of training points to be exhausted before training set refinement; - 'nTrainingPoints': number of training points; - 'nTestPoints': number of starting test points; - 'trainingSetGenerator': training sample points generator. - 'robustTol': tolerance for robust Pade' denominator management. extraApproxParameters: List of approxParameters keys in addition to mother class's. POD: whether to compute POD of snapshots. muBounds: list of bounds for parameter values. greedyTol: uniform error tolerance for greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTrainingPoints: number of training points. nTestPoints: number of starting test points. trainingSetGenerator: training sample points generator. robustTol: tolerance for robust Pade' denominator management. samplingEngine: Sampling engine. 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. """ TOL_INSTABILITY = 1e-6 def __init__(self, HFEngine:HFEng, mu0 : complex = 0., approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["muBounds","greedyTol", "maxIter", "refinementRatio", "nTrainingPoints", "nTestPoints", "trainingSetGenerator"]) super(GenericApproximantLagrange, self).__init__( HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """Value of approximant parameters. Its assignment may change S.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["muBounds","greedyTol", "maxIter", "refinementRatio", "nTrainingPoints", "nTestPoints", "trainingSetGenerator"], True, True, baselevel = 1) GenericApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "muBounds" in keyList: self.muBounds = approxParameters["muBounds"] elif hasattr(self, "muBounds"): self.muBounds = self.muBounds else: self.muBounds = [[0.], [1.]] if "greedyTol" in keyList: self.greedyTol = approxParameters["greedyTol"] elif hasattr(self, "greedyTol"): self.greedyTol = self.greedyTol else: self.greedyTol = 1e-2 if "maxIter" in keyList: self.maxIter = approxParameters["maxIter"] elif hasattr(self, "maxIter"): self.maxIter = self.maxIter else: self.maxIter = 1e2 if "refinementRatio" in keyList: self.refinementRatio = approxParameters["refinementRatio"] elif hasattr(self, "refinementRatio"): self.refinementRatio = self.refinementRatio else: self.refinementRatio = 0.2 if "nTrainingPoints" in keyList: self.nTrainingPoints = approxParameters["nTrainingPoints"] elif hasattr(self, "nTrainingPoints"): self.nTrainingPoints = self.nTrainingPoints else: self.nTrainingPoints = np.int(np.ceil(self.maxIter / self.refinementRatio)) if "nTestPoints" in keyList: self.nTestPoints = approxParameters["nTestPoints"] elif hasattr(self, "nTestPoints"): self.nTestPoints = self.nTestPoints else: self.nTestPoints = 1 if "trainingSetGenerator" in keyList: self.trainingSetGenerator = ( approxParameters["trainingSetGenerator"]) elif hasattr(self, "trainingSetGenerator"): self.trainingSetGenerator = self.trainingSetGenerator else: from rrompy.utilities.parameter_sampling import QuadratureSampler self.trainingSetGenerator = QuadratureSampler(self.muBounds, "UNIFORM") del QuadratureSampler @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): self._S = S @property def mus(self): """Value of mus.""" return self._mus @mus.setter def mus(self, mus): self._mus = np.array(mus, dtype = np.complex) @property def muBounds(self): """Value of muBounds.""" return self._muBounds @muBounds.setter def muBounds(self, muBounds): if len(muBounds) != 2: raise Exception("2 limits must be specified.") try: muBounds = muBounds.tolist() except: muBounds = list(muBounds) for j in range(2): try: len(muBounds[j]) except: muBounds[j] = np.array([muBounds[j]]) if len(muBounds[0]) != len(muBounds[1]): raise Exception("The bounds must have the same length.") self._muBounds = muBounds @property def greedyTol(self): """Value of greedyTol.""" return self._greedyTol @greedyTol.setter def greedyTol(self, greedyTol): if greedyTol < 0: raise ArithmeticError("greedyTol must be non-negative.") if hasattr(self, "greedyTol"): greedyTolold = self.greedyTol else: greedyTolold = -1 self._greedyTol = greedyTol self._approxParameters["greedyTol"] = self.greedyTol if greedyTolold != self.greedyTol: self.resetSamples() @property def maxIter(self): """Value of maxIter.""" return self._maxIter @maxIter.setter def maxIter(self, maxIter): if maxIter <= 0: raise ArithmeticError("maxIter must be positive.") if hasattr(self, "maxIter"): maxIterold = self.maxIter else: maxIterold = -1 self._maxIter = maxIter self._approxParameters["maxIter"] = self.maxIter if maxIterold != self.maxIter: self.resetSamples() @property def refinementRatio(self): """Value of refinementRatio.""" return self._refinementRatio @refinementRatio.setter def refinementRatio(self, refinementRatio): if refinementRatio <= 0. or refinementRatio > 1.: raise ArithmeticError(("refinementRatio must be between 0 " "(excluded) and 1.")) if hasattr(self, "refinementRatio"): refinementRatioold = self.refinementRatio else: refinementRatioold = -1 self._refinementRatio = refinementRatio self._approxParameters["refinementRatio"] = self.refinementRatio if refinementRatioold != self.refinementRatio: self.resetSamples() @property def nTrainingPoints(self): """Value of nTrainingPoints.""" return self._nTrainingPoints @nTrainingPoints.setter def nTrainingPoints(self, nTrainingPoints): if nTrainingPoints <= 1: raise ArithmeticError("nTrainingPoints must be greater than 1.") if not np.isclose(nTrainingPoints, np.int(nTrainingPoints)): raise ArithmeticError("nTrainingPoints must be an integer.") nTrainingPoints = np.int(nTrainingPoints) if hasattr(self, "nTrainingPoints"): nTrainingPointsold = self.nTrainingPoints else: nTrainingPointsold = -1 self._nTrainingPoints = nTrainingPoints self._approxParameters["nTrainingPoints"] = self.nTrainingPoints if nTrainingPointsold != self.nTrainingPoints: self.resetSamples() @property def nTestPoints(self): """Value of nTestPoints.""" return self._nTestPoints @nTestPoints.setter def nTestPoints(self, nTestPoints): if nTestPoints <= 0: raise ArithmeticError("nTestPoints must be positive.") if not np.isclose(nTestPoints, np.int(nTestPoints)): raise ArithmeticError("nTestPoints must be an integer.") nTestPoints = np.int(nTestPoints) if hasattr(self, "nTestPoints"): nTestPointsold = self.nTestPoints else: nTestPointsold = -1 self._nTestPoints = nTestPoints self._approxParameters["nTestPoints"] = self.nTestPoints if nTestPointsold != self.nTestPoints: self.resetSamples() @property def trainingSetGenerator(self): """Value of trainingSetGenerator.""" return self._trainingSetGenerator @trainingSetGenerator.setter def trainingSetGenerator(self, trainingSetGenerator): if 'generatePoints' not in dir(trainingSetGenerator): raise Exception("trainingSetGenerator type not recognized.") if hasattr(self, '_trainingSetGenerator'): trainingSetGeneratorOld = self.trainingSetGenerator self._trainingSetGenerator = trainingSetGenerator self._approxParameters["trainingSetGenerator"] = ( self.trainingSetGenerator) if (not 'trainingSetGeneratorOld' in locals() or trainingSetGeneratorOld != self.trainingSetGenerator): self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() self._mus = [] def initEstNormer(self): """Initialize estimator norm engine.""" if not hasattr(self, "estNormer"): from rrompy.hfengines.base import ProblemEngineBase as PEB -# class L2InverseNormer(PEB): -# def innerProduct(self, u:Np1D, v:Np1D) -> float: -# if not hasattr(self, "energyNormMatrix"): -# self.buildEnergyNormForm() -# return v.conj().T.dot(self.energyNormMatrix.solve(u)) - -# def buildEnergyNormForm(self): -# super().buildEnergyNormForm() -# from scipy.sparse import csc_matrix, linalg as sla -# Mass = csc_matrix(self.energyNormMatrix, -# dtype = np.complex) -# self.energyNormMatrix = sla.spilu(Mass) -# self.estNormer = L2InverseNormer() # inverse L2 norm self.estNormer = PEB() # L2 norm self.estNormer.V = self.HFEngine.V self.estNormer.buildEnergyNormForm() def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]: """ Standard residual-based error estimator with explicit residual computation. """ self.setupApprox() nmus = len(mus) err = np.empty(nmus) if self.HFEngine.nbs == 1: RHS = self.getRHS(mus[0], homogeneized = self.homogeneized) RHSNorm = self.estNormer.norm(RHS) for j in range(nmus): res = self.getRes(mus[j], homogeneized = self.homogeneized) err[j] = self.estNormer.norm(res) / RHSNorm else: for j in range(nmus): res = self.getRes(mus[j], homogeneized = self.homogeneized) RHS = self.getRHS(mus[j], homogeneized = self.homogeneized) err[j] = self.estNormer.norm(res) / self.estNormer.norm(RHS) return np.abs(err) def getMaxErrorEstimator(self) -> Tuple[float, int]: """ Compute maximum of (and index of maximum of) error estimator over training set. """ errorEstTrain = self.errorEstimator(self.muTrain) idxMaxEst = np.argmax(errorEstTrain) maxEst = errorEstTrain[idxMaxEst] return maxEst, idxMaxEst def greedyNextSample(self, muidx:int, plotEst : bool = False): """Compute next greedy snapshot of solution map.""" mu = self.muTrain[muidx] if self.verbosity >= 2: verbosityDepth("MAIN", ("Adding {}-th sample point at {} to test " "set.").format( self.samplingEngine.nsamples + 1, mu)) self.mus = np.append(self.mus, mu) idxs = np.arange(len(self.muTrain)) mask = np.ones_like(idxs, dtype = bool) mask[muidx] = False idxs = idxs[mask] self.muTrain = self.muTrain[idxs] self.samplingEngine.nextSample(mu, homogeneized = self.homogeneized) errorEstTrain = self.errorEstimator(self.muTrain) muidx = np.argmax(errorEstTrain) maxErrorEst = errorEstTrain[muidx] mu = self.muTrain[muidx] if plotEst and not np.all(np.isinf(errorEstTrain)): from matplotlib import pyplot as plt plt.figure() plt.semilogy(np.real(self.muTrain), errorEstTrain, 'k') plt.semilogy(np.real(self.muTrain), self.greedyTol * np.ones(len(self.muTrain)), 'r--') plt.semilogy(np.real(self.mus), 2. * self.greedyTol * np.ones(len(self.mus)), '*m') plt.semilogy(np.real(mu), maxErrorEst, 'xr') plt.grid() plt.show() plt.close() return errorEstTrain, muidx, maxErrorEst, mu def greedy(self, plotEst : bool = False): """Compute greedy snapshots of solution map.""" if self.samplingEngine.samples is None: if self.verbosity >= 2: verbosityDepth("INIT", "Starting computation of snapshots.") self.resetSamples() self.mus, _ = self.trainingSetGenerator.generatePoints( self.nTestPoints) self.mus = np.array([x[0] for x in self.mus], dtype = np.complex) nTrain = self.nTrainingPoints muTrainBase, _ = self.trainingSetGenerator.generatePoints(nTrain) self.muTrain = np.empty(len(muTrainBase) + 1, dtype = np.complex) j = 0 for mu in muTrainBase: if not np.any(np.isclose(self.mus, mu)): self.muTrain[j] = mu[0] j += 1 self.muTrain[j] = self.mus[-1] self.muTrain = self.muTrain[: j + 1] self.mus = self.mus[:-1] for j in range(len(self.mus)): if self.verbosity >= 2: verbosityDepth("MAIN", ("Adding {}-th sample point at {} " "to test set.").format( self.samplingEngine.nsamples + 1, self.mus[j])) self.samplingEngine.nextSample(self.mus[j], homogeneized = self.homogeneized) errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample(-1, plotEst) if self.verbosity >= 2: verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\ .format(maxErrorEst)) while (self.samplingEngine.nsamples < self.maxIter and (maxErrorEst > self.greedyTol or self.samplingEngine.nsamples < self.nTestPoints)): if (1. - self.refinementRatio) * nTrain > len(self.muTrain): muTrainExtra, _ = self.trainingSetGenerator.refine(nTrain) self.muTrain = np.sort(np.append(self.muTrain, muTrainExtra)) nTrain += len(muTrainExtra) if self.verbosity >= 5: verbosityDepth("MAIN", ("Enriching training set by {} " "elements.").format( len(muTrainExtra))) muTrainOld, maxErrorEstOld = self.muTrain, maxErrorEst errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample( muidx, plotEst) if self.verbosity >= 2: verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\ .format(maxErrorEst)) if (np.isnan(maxErrorEst) or np.isinf(maxErrorEst) or maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY): warn(("Instability in a posteriori estimator. Starting " "preemptive greedy loop termination.")) maxErrorEst = maxErrorEstOld self.muTrain = muTrainOld self.mus = self.mus[:-1] self.samplingEngine.popSample() self.setupApprox() break if self.verbosity >= 2: verbosityDepth("DEL", ("Done computing snapshots (final " "snapshot count: {}).").format( self.samplingEngine.nsamples)) 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, "_S") and self.S == len(self.mus) and super().checkComputedApprox()) def computeScaleFactor(self): """Compute parameter rescaling factor.""" self.scaleFactor = .5 * np.abs(self.HFEngine.rescaling( self.trainingSetGenerator.lims[0][0]) - self.HFEngine.rescaling( self.trainingSetGenerator.lims[1][0])) diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py index 647a1d5..54f9e09 100644 --- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py +++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py @@ -1,546 +1,558 @@ # 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 numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_approximant_taylor import GenericApproximantTaylor from rrompy.sampling.base.pod_engine import PODEngine from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, DictAny from rrompy.utilities.base.types import HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantTaylorPade'] class ApproximantTaylorPade(GenericApproximantTaylor): """ ROM single-point fast Pade' approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'rho': weight for computation of original Pade' approximant; defaults to np.inf, i.e. fast approximant; - 'M': degree of Pade' approximant numerator; defaults to 0; - 'N': degree of Pade' approximant denominator; defaults to 0; - 'E': total number of derivatives current approximant relies upon; defaults to Emax; - 'Emax': total number of derivatives of solution map to be computed; defaults to E; - 'robustTol': tolerance for robust Pade' denominator management; defaults to 0; - 'sampleType': label of sampling type; available values are: - 'ARNOLDI': orthogonalization of solution derivatives through Arnoldi algorithm; - 'KRYLOV': standard computation of solution derivatives. Defaults to 'KRYLOV'. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'rho': weight for computation of original Pade' approximant; - 'M': degree of Pade' approximant numerator; - 'N': degree of Pade' approximant denominator; - 'E': total number of derivatives current approximant relies upon; - 'Emax': total number of derivatives of solution map to be computed; - 'robustTol': tolerance for robust Pade' denominator management; - 'sampleType': label of sampling type. 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. initialHFData: HF problem initial data. uHF: High fidelity solution with wavenumber lastSolvedHF as numpy complex vector. lastSolvedHF: Wavenumber corresponding to last computed high fidelity solution. G: Square Numpy 2D vector of size (N+1) corresponding to Pade' denominator matrix (see paper). Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. uApp: Last evaluated approximant as numpy complex vector. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["M", "N", "robustTol", "rho"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["M", "N", "robustTol", "rho"], True, True, baselevel = 1) keyList = list(approxParameters.keys()) if "rho" in keyList: self._rho = approxParameters["rho"] elif hasattr(self, "rho"): self._rho = self.rho else: self._rho = np.inf GenericApproximantTaylor.approxParameters.fset(self, approxParametersCopy) self.rho = self._rho if "robustTol" in keyList: self.robustTol = approxParameters["robustTol"] elif hasattr(self, "robustTol"): self.robustTol = self.robustTol else: self.robustTol = 0 self._ignoreParWarnings = True if "M" in keyList: self.M = approxParameters["M"] elif hasattr(self, "M"): self.M = self.M else: self.M = 0 del self._ignoreParWarnings if "N" in keyList: self.N = approxParameters["N"] elif hasattr(self, "N"): self.N = self.N else: self.N = 0 @property def rho(self): """Value of rho.""" return self._rho @rho.setter def rho(self, rho): self._rho = np.abs(rho) self._approxParameters["rho"] = self.rho @property def M(self): """Value of M. Its assignment may change Emax and E.""" return self._M @M.setter def M(self, M): if M < 0: raise ArithmeticError("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if not hasattr(self, "_ignoreParWarnings"): self.checkMNEEmax() @property def N(self): """Value of N. Its assignment may change Emax.""" return self._N @N.setter def N(self, N): if N < 0: raise ArithmeticError("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N self.checkMNEEmax() def checkMNEEmax(self): """Check consistency of M, N, E, and Emax.""" M = self.M if hasattr(self, "_M") else 0 N = self.N if hasattr(self, "_N") else 0 E = self.E if hasattr(self, "_E") else M + N Emax = self.Emax if hasattr(self, "_Emax") else M + N if self.rho == np.inf: if Emax < max(N, M): warn(("Prescribed Emax is too small. Updating Emax to " "max(M, N).")) self.Emax = max(N, M) if E < max(N, M): warn("Prescribed E is too small. Updating E to max(M, N).") self.E = max(N, M) else: if Emax < N + M: warn("Prescribed Emax is too small. Updating Emax to M + N.") self.Emax = self.N + M if E < N + M: warn("Prescribed E is too small. Updating E to M + N.") 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.: warn("Overriding prescribed negative robustness tolerance to 0.") 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: warn("Prescribed Emax is too small. Updating Emax to E.") 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: warn(("Prescribed Emax is too small. Updating Emax to " "max(N, M, E).")) Emax = max(N, M, E) else: if max(N + M, E) > Emax: warn(("Prescribed Emax is too small. Updating Emax to " "max(N + M, E).")) Emax = max(N + M, E) self._Emax = Emax self._approxParameters["Emax"] = Emax if EmaxOld >= self.Emax and self.samplingEngine.samples is not None: self.samplingEngine.samples = self.samplingEngine.samples[:, : self.Emax + 1] if (self.sampleType == "ARNOLDI" and self.samplingEngine.HArnoldi is not None): self.samplingEngine.HArnoldi = self.samplingEngine.HArnoldi[ : self.Emax + 1, : self.Emax + 1] self.samplingEngine.RArnoldi = self.samplingEngine.RArnoldi[ : self.Emax + 1, : self.Emax + 1] def setupApprox(self): """ Compute Pade' approximant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeDerivatives() if self.N > 0: if self.verbosity >= 7: verbosityDepth("INIT", ("Starting computation of " "denominator.")) while self.N > 0: if self.POD: ev, eV = self.findeveVGQR() else: ev, eV = self.findeveVGExplicit() newParameters = checkRobustTolerance(ev, self.E, self.robustTol) if not newParameters: break self.approxParameters = newParameters if self.N <= 0: eV = np.ones((1, 1)) self.Q = np.poly1d(eV[:, 0]) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing denominator.") else: self.Q = np.poly1d([1]) if self.verbosity >= 7: verbosityDepth("INIT", "Starting computation of numerator.") self.P = np.zeros((self.Emax + 1, self.M + 1), dtype = np.complex) for i in range(self.Q.order): rng = np.arange(self.M + 1 - i) self.P[rng, - 1 - rng - i] = self.Q[i] if self.sampleType == "ARNOLDI": self.P = self.samplingEngine.RArnoldi.dot(self.P) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing numerator.") self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def rescaleParameter(self, R:Np2D, A : Np2D = None, exponent : float = 1.) -> Np2D: """ Prepare parameter rescaling. Args: R: Matrix whose columns need rescaling. A(optional): Matrix whose diagonal defines scaling factor. If None, previous value of scaleFactor is used. Defaults to None. exponent(optional): Exponent of scaling factor in matrix diagonal. Defaults to 1. Returns: Rescaled matrix. """ if A is not None: aDiag = np.diag(A) scaleCoeffs = np.polyfit(np.arange(A.shape[1]), np.log(aDiag), 1) self.scaleFactor = np.exp(- scaleCoeffs[0] / exponent) return np.multiply(R, np.power(self.scaleFactor,np.arange(R.shape[1]))) def buildG(self): """Assemble Pade' denominator matrix.""" self.computeDerivatives() if self.rho == np.inf: Nmin = self.E - self.N else: Nmin = self.M - self.N + 1 if self.sampleType == "KRYLOV": DerE = self.samplingEngine.samples[:, Nmin : self.E + 1] G = self.HFEngine.innerProduct(DerE, DerE) DerE = self.rescaleParameter(DerE, G, 2.) G = self.HFEngine.innerProduct(DerE, DerE) else: RArnE = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1] RArnE = self.rescaleParameter(RArnE, RArnE[Nmin :, :]) G = RArnE.conj().T.dot(RArnE) if self.rho == np.inf: self.G = G else: Gbig = G self.G = np.zeros((self.N + 1, self.N + 1), dtype = np.complex) for k in range(self.E - self.M): 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. """ if self.verbosity >= 10: verbosityDepth("INIT", "Building gramian matrix.") self.buildG() if self.verbosity >= 10: verbosityDepth("DEL", "Done building gramian.") if self.verbosity >= 7: verbosityDepth("INIT", "Solving eigenvalue problem for gramian matrix.") ev, eV = np.linalg.eigh(self.G) eV = self.rescaleParameter(eV.T).T if self.verbosity >= 5: try: condev = ev[-1] / ev[0] except: condev = np.inf verbosityDepth("MAIN", ("Solved eigenvalue problem of size {} " "with condition number {:.4e}.").format( self.N + 1, condev)) if self.verbosity >= 7: verbosityDepth("DEL", "Done solving eigenvalue problem.") 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.rho == np.inf: Nmin = self.E - self.N else: Nmin = self.M - self.N + 1 if self.sampleType == "KRYLOV": A = copy(self.samplingEngine.samples[:, Nmin : self.E + 1]) self.PODEngine = PODEngine(self.HFEngine) if self.verbosity >= 10: verbosityDepth("INIT", "Orthogonalizing samples.", end = "") R = self.PODEngine.QRHouseholder(A, only_R = True) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) else: R = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1] R = self.rescaleParameter(R, R[Nmin :, :]) if self.rho == np.inf: if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square root of " "gramian matrix.")) sizeI = R.shape[0] _, s, V = np.linalg.svd(R, full_matrices = False) else: if self.verbosity >= 10: verbosityDepth("INIT", ("Building matrix stack for square " "root of gramian."), end = "") 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]) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square root of " "gramian matrix.")) sizeI = Rtower.shape[0] _, s, V = np.linalg.svd(Rtower, full_matrices = False) eV = V.conj().T[:, ::-1] eV = self.rescaleParameter(eV.T).T if self.verbosity >= 5: try: condev = s[0] / s[-1] except: condev = np.inf verbosityDepth("MAIN", ("Solved svd problem of size {} x {} with " "condition number {:.4e}.").format( sizeI, self.N + 1, condev)) if self.verbosity >= 7: verbosityDepth("DEL", "Done solving eigenvalue problem.") return s[::-1], eV 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 self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0) def getPVal(self, mu:List[complex]): """ Evaluate Pade' numerator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating numerator at mu = {}.".format(mu), end = "") try: len(mu) except: mu = [mu] powerlist = np.vander(self.radiusPade(mu), self.M + 1).T p = self.P.dot(powerlist) if len(mu) == 1: p = p.flatten() if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return p def getQVal(self, mu:List[complex]): """ Evaluate Pade' denominator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating denominator at mu = {}.".format(mu), end = "") q = self.Q(self.radiusPade(mu)) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return q def evalApproxReduced(self, mu:complex): """ Evaluate Pade' approximant at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Evaluating approximant at mu = {}.".format(mu)) self.uAppReduced = self.getPVal(mu) / self.getQVal(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done evaluating approximant.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.samplingEngine.samples.dot(self.uAppReduced) + def getResidues(self) -> Np1D: + """ + Obtain norm of approximant residues. + + Returns: + Numpy vector of norms of residues. + """ + poles = self.getPoles() + Pvals = self.samplingEngine.samples.dot(self.getPVal(poles)) + Qder = - self.Q.deriv(1)(self.radiusPade(poles)) + return Pvals / Qder + def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() return self.HFEngine.rescalingInv(self.Q.r + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py index fc7cac1..ec50f00 100644 --- a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py +++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py @@ -1,244 +1,245 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.base.types import DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['GenericApproximantTaylor'] class GenericApproximantTaylor(GenericApproximant): """ ROM single-point approximant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'E': total number of derivatives current approximant relies upon; defaults to Emax; - 'Emax': total number of derivatives of solution map to be computed; defaults to E; - 'sampleType': label of sampling type; available values are: - 'ARNOLDI': orthogonalization of solution derivatives through Arnoldi algorithm; - 'KRYLOV': standard computation of solution derivatives. Defaults to 'KRYLOV'. Defaults to empty dict. homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'E': total number of derivatives current approximant relies upon; - 'Emax': total number of derivatives of solution map to be computed; - 'sampleType': label of sampling type. 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. initialHFData: HF problem initial data. samplingEngine: Sampling engine. 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. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self._addParametersToList(["E", "Emax", "sampleType"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity) self._postInit() def setupSampling(self): """Setup sampling engine.""" if not hasattr(self, "sampleType"): return if self.sampleType == "ARNOLDI": - from rrompy.sampling.scipy.sampling_engine_arnoldi import ( - SamplingEngineArnoldi) + from rrompy.sampling.linear_problem.sampling_engine_arnoldi \ + import SamplingEngineArnoldi super().setupSampling(SamplingEngineArnoldi) elif self.sampleType == "KRYLOV": - from rrompy.sampling.scipy.sampling_engine_krylov import ( - SamplingEngineKrylov) + from rrompy.sampling.linear_problem.sampling_engine_krylov \ + import SamplingEngineKrylov super().setupSampling(SamplingEngineKrylov) else: raise Exception("Sample type not recognized.") @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change E and Emax. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["E", "Emax", "sampleType"], True, True, baselevel = 1) GenericApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "E" in keyList: self._E = approxParameters["E"] self._approxParameters["E"] = self.E if "Emax" in keyList: self.Emax = approxParameters["Emax"] else: if not hasattr(self, "Emax"): self.Emax = self.E else: self.Emax = self.Emax else: if "Emax" in keyList: self._E = approxParameters["Emax"] self._approxParameters["E"] = self.E self.Emax = self.E else: if not (hasattr(self, "Emax") and hasattr(self, "E")): raise Exception("At least one of E and Emax must be set.") if "sampleType" in keyList: self.sampleType = approxParameters["sampleType"] elif hasattr(self, "sampleType"): self.sampleType = self.sampleType else: self.sampleType = "KRYLOV" @property def E(self): """Value of E. Its assignment may change Emax.""" return self._E @E.setter def E(self, E): if E < 0: raise ArithmeticError("E must be non-negative.") self._E = E self._approxParameters["E"] = self.E if hasattr(self, "Emax") and self.Emax < self.E: warn("Prescribed E is too large. Updating Emax to E.") 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: warn("Prescribed Emax is too small. Updating Emax to E.") self.Emax = self.E else: self._approxParameters["Emax"] = self.Emax if (EmaxOld >= self.Emax and self.samplingEngine.samples is not None): self.samplingEngine.samples = self.samplingEngine.samples[:, : self.Emax + 1] if (self.sampleType == "ARNOLDI" and self.samplingEngine.HArnoldi is not None): self.samplingEngine.HArnoldi= self.samplingEngine.HArnoldi[ : self.Emax + 1, : self.Emax + 1] self.samplingEngine.RArnoldi= self.samplingEngine.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 Exception("Sample type not recognized.") self._sampleType = sampleType except: warn(("Prescribed sampleType not recognized. Overriding to " "'KRYLOV'.")) 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.samplingEngine.samples is None: if self.verbosity >= 5: verbosityDepth("INIT", "Starting computation of derivatives.") self.samplingEngine.iterSample(self.mu0, self.Emax + 1, homogeneized = self.homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Done computing derivatives.") 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.samplingEngine.samples is not None and super().checkComputedApprox()) def normApp(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ if self.sampleType != "ARNOLDI" or self.homogeneized != homogeneized: return super().normApp(mu, homogeneized) return np.linalg.norm(self.getAppReduced(mu, homogeneized)) + diff --git a/rrompy/sampling/scipy/__init__.py b/rrompy/sampling/linear_problem/__init__.py similarity index 100% rename from rrompy/sampling/scipy/__init__.py rename to rrompy/sampling/linear_problem/__init__.py diff --git a/rrompy/sampling/scipy/sampling_engine_arnoldi.py b/rrompy/sampling/linear_problem/sampling_engine_arnoldi.py similarity index 100% rename from rrompy/sampling/scipy/sampling_engine_arnoldi.py rename to rrompy/sampling/linear_problem/sampling_engine_arnoldi.py diff --git a/rrompy/sampling/scipy/sampling_engine_krylov.py b/rrompy/sampling/linear_problem/sampling_engine_krylov.py similarity index 100% rename from rrompy/sampling/scipy/sampling_engine_krylov.py rename to rrompy/sampling/linear_problem/sampling_engine_krylov.py diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange.py b/rrompy/sampling/linear_problem/sampling_engine_lagrange.py similarity index 100% rename from rrompy/sampling/scipy/sampling_engine_lagrange.py rename to rrompy/sampling/linear_problem/sampling_engine_lagrange.py diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py b/rrompy/sampling/linear_problem/sampling_engine_lagrange_pod.py similarity index 100% rename from rrompy/sampling/scipy/sampling_engine_lagrange_pod.py rename to rrompy/sampling/linear_problem/sampling_engine_lagrange_pod.py diff --git a/rrompy/utilities/base/__init__.py b/rrompy/utilities/base/__init__.py index 90a85f8..00051cd 100644 --- a/rrompy/utilities/base/__init__.py +++ b/rrompy/utilities/base/__init__.py @@ -1,46 +1,44 @@ # 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 .find_dict_str_key import findDictStrKey from .get_new_filename import getNewFilename from .purge_dict import purgeDict from .purge_list import purgeList from .number_theory import (squareResonances, primeFactorize, getLowestPrimeFactor) from .sobol import sobolGenerate from . import types as Types -from . import fenics as Fenics from .verbosity_depth import verbosityDepth from .custom_fit import customFit __all__ = [ 'findDictStrKey', 'getNewFilename', 'purgeDict', 'purgeList', 'squareResonances', 'primeFactorize', 'getLowestPrimeFactor', 'sobolGenerate', 'Types', - 'Fenics', 'verbosityDepth', 'customFit' ] diff --git a/rrompy/utilities/base/fenics.py b/rrompy/utilities/base/fenics.py index 955b774..271e8c5 100644 --- a/rrompy/utilities/base/fenics.py +++ b/rrompy/utilities/base/fenics.py @@ -1,30 +1,32 @@ # 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'] +__all__ = ['fenZERO','fenZEROS','fenONE','fenONES','bdrTrue','bdrFalse'] # CONSTANTS fenZERO = fen.Constant(0.) +fenZEROS = lambda n: fen.Constant(tuple([0.] * n)) fenONE = fen.Constant(1.) +fenONES = lambda n: fen.Constant(tuple([1.] * n)) # BOUNDARY BOOLS bdrTrue = lambda x, on_boundary : on_boundary bdrFalse = lambda x, on_boundary : False