diff --git a/examples/1_symmetric_disk/symmetric_disk.py b/examples/1_symmetric_disk/symmetric_disk.py index ca9e8a7..e5d9026 100644 --- a/examples/1_symmetric_disk/symmetric_disk.py +++ b/examples/1_symmetric_disk/symmetric_disk.py @@ -1,65 +1,86 @@ import numpy as np from symmetric_disk_engine import SymmetricDiskEngine as engine -from rrompy.reduction_methods.standard import (RationalInterpolant as RI, +from rrompy.sampling.standard import SamplingEngineStandard as SES +from rrompy.reduction_methods.standard import (NearestNeighbor as NN, + RationalInterpolant as RI, ReducedBasis as RB) from rrompy.reduction_methods.greedy import (RationalInterpolantGreedy as RIG, ReducedBasisGreedy as RBG) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS ks = [10., 20.] k0, n = np.mean(np.power(ks, 2.)) ** .5, 150 solver = engine(k0, n) k = 12. for method in ["RI", "RB", "RI_GREEDY", "RB_GREEDY"]: print("Testing {} method".format(method)) if method == "RI": params = {'N':39, 'M':39, 'S':40, 'POD':True, 'polybasis':"CHEBYSHEV", 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RI if method == "RB": params = {'R':40, 'S':40, 'POD':True, 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RB if method == "RI_GREEDY": params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, 'sampler':QS(ks, "UNIFORM", scalingExp = 2.), 'errorEstimatorKind':"DISCREPANCY"} algo = RIG if method == "RB_GREEDY": params = {'S':10, 'POD':True, 'greedyTol':1e-2, 'sampler':QS(ks, "UNIFORM", scalingExp = 2.)} algo = RBG approx = algo(solver, mu0 = k0, approx_state = True, approxParameters = params, verbosity = 20) if len(method) == 2: approx.setupApprox() else: approx.setupApprox("LAST") + print("--- Approximant ---") approx.plotApprox(k, name = 'u_app') approx.plotHF(k, name = 'u_HF') - approx.plotErr(k, name = 'err') - approx.plotRes(k, name = 'res') + approx.plotErr(k, name = 'err_app') + approx.plotRes(k, name = 'res_app') normErr = approx.normErr(k)[0] normSol = approx.normHF(k)[0] normRes = approx.normRes(k)[0] normRHS = approx.normRHS(k)[0] - print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format( normSol, normErr, normErr / normSol)) - print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + print("--- Closest snapshot ---") + eng = SES(solver, verbosity = 0) + approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0, + approxParameters = {'S':approx.S, 'POD':True}) + approxNN.setSamples(approx.samplingEngine) + approxNN.plotApprox(k, name = 'u_close') + approxNN.plotHF(k, name = 'u_HF') + approxNN.plotErr(k, name = 'err_close') + approxNN.plotRes(k, name = 'res_close') + normErr = approxNN.normErr(k)[0] + normSol = approxNN.normHF(k)[0] + normRes = approxNN.normRes(k)[0] + normRHS = approxNN.normRHS(k)[0] + print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format( + normSol, normErr, normErr / normSol)) + print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format( + normRHS, normRes, normRes / normRHS)) + if method[:2] == "RI": poles, residues = approx.getResidues() if method[:2] == "RB": poles = approx.getPoles() print("Poles:\n{}".format(poles)) if method[:2] == "RI": for pol, res in zip(poles, residues.T): solver.plot(res) print("pole = {:.5e}".format(pol)) print("\n") diff --git a/examples/2_double_slit/double_slit.py b/examples/2_double_slit/double_slit.py index e379bc3..85d2adb 100644 --- a/examples/2_double_slit/double_slit.py +++ b/examples/2_double_slit/double_slit.py @@ -1,62 +1,83 @@ import numpy as np from double_slit_engine import DoubleSlitEngine as engine -from rrompy.reduction_methods.standard import RationalInterpolant as RI +from rrompy.sampling.standard import SamplingEngineStandard as SES +from rrompy.reduction_methods.standard import (NearestNeighbor as NN, + RationalInterpolant as RI) from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG from rrompy.parameter.parameter_sampling import QuadratureSampler as QS from rrompy.solver.fenics import interp_project ks = [10., 15.] k0, n = np.mean(ks), 150 solver = engine(k0, n) k = 11. for method in ["RI", "RI_GREEDY"]: print("Testing {} method".format(method)) if method == "RI": params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV", 'sampler':QS(ks, "CHEBYSHEV")} algo = RI if method == "RI_GREEDY": params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, 'sampler':QS(ks, "UNIFORM"), 'errorEstimatorKind':"LOOK_AHEAD", 'trainSetGenerator':QS(ks, "CHEBYSHEV")} algo = RIG approx = algo(solver, mu0 = k0, approx_state = True, approxParameters = params, verbosity = 20) if len(method) == 2: approx.setupApprox() else: approx.setupApprox("LAST") + print("--- Approximant ---") approx.plotApprox(k, name = 'u_app') approx.plotHF(k, name = 'u_HF') - approx.plotErr(k, name = 'err') - approx.plotRes(k, name = 'res') + approx.plotErr(k, name = 'err_app') + approx.plotRes(k, name = 'res_app') normErr = approx.normErr(k)[0] normSol = approx.normHF(k)[0] normRes = approx.normRes(k)[0] normRHS = approx.normRHS(k)[0] - print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format( normSol, normErr, normErr / normSol)) - print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + print("--- Closest snapshot ---") + eng = SES(solver, verbosity = 0) + approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0, + approxParameters = {'S':approx.S, 'POD':True}) + approxNN.setSamples(approx.samplingEngine) + approxNN.plotApprox(k, name = 'u_close') + approxNN.plotHF(k, name = 'u_HF') + approxNN.plotErr(k, name = 'err_close') + approxNN.plotRes(k, name = 'res_close') + normErr = approxNN.normErr(k)[0] + normSol = approxNN.normHF(k)[0] + normRes = approxNN.normRes(k)[0] + normRHS = approxNN.normRHS(k)[0] + print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format( + normSol, normErr, normErr / normSol)) + print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format( + normRHS, normRes, normRes / normRHS)) + uIncR, uIncI = solver.getDirichletValues(k) uIncR = interp_project(uIncR, solver.V) uIncI = interp_project(uIncI, solver.V) uInc = np.array(uIncR.vector()) + 1.j * np.array(uIncI.vector()) uEx = approx.getHF(k)[0] - uInc uApp = approx.getApprox(k)[0] - uInc solver.plot(uEx, name = 'uex_tot') solver.plot(uApp, name = 'u_app_tot') poles, residues = approx.getResidues() print("Poles:\n{}".format(poles)) for pol, res in zip(poles, residues.T): solver.plot(res) print("pole = {:.5e}".format(pol)) print("\n") diff --git a/examples/3_sector_angle/sector_angle.py b/examples/3_sector_angle/sector_angle.py index 6ffd28d..eb73ecf 100644 --- a/examples/3_sector_angle/sector_angle.py +++ b/examples/3_sector_angle/sector_angle.py @@ -1,91 +1,114 @@ import numpy as np import matplotlib.pyplot as plt from sector_angle_engine import SectorAngleEngine as engine +from rrompy.sampling.standard import SamplingEngineStandard as SES +from rrompy.reduction_methods.standard import NearestNeighbor as NN from rrompy.reduction_methods.pivoted import ( RationalInterpolantPivoted as RIP, RationalInterpolantGreedyPivoted as RIGP) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS ks, ts = [10., 15.], [.4, .6] k0, t0, n = np.mean(np.power(ks, 2.)) ** .5, np.mean(ts), 50 solver = engine(k0, t0, n) murange = [[ks[0], ts[0]], [ks[-1], ts[-1]]] mu = [12., .535] +fighandles = [] for method in ["RI", "RI_GREEDY"]: print("Testing {} method".format(method)) if method == "RI": - params = {'N':19, 'M':19, 'S':20, 'MMarginal':1, 'SMarginal':7, + params = {'N':19, 'M':19, 'S':20, 'MMarginal':3, 'SMarginal':11, 'POD':True, 'polybasis':"CHEBYSHEV", 'polybasisMarginal':"MONOMIAL_GAUSSIAN", - 'matchingWeight':1., 'cutOffTolerance': 1., + 'radialDirectionalWeightsMarginal': 10., + 'matchingWeight':1., 'cutOffTolerance': 2., 'samplerPivot':QS(ks, "CHEBYSHEV", 2.), 'samplerMarginal':QS(ts, "UNIFORM")} algo = RIP if method == "RI_GREEDY": - params = {'S':10, 'MMarginal':1, 'SMarginal':7, 'POD':True, + params = {'S':10, 'MMarginal':3, 'SMarginal':11, 'POD':True, 'polybasis':"LEGENDRE", 'polybasisMarginal':"MONOMIAL_GAUSSIAN", - 'matchingWeight':1., 'cutOffTolerance': 1., + 'radialDirectionalWeightsMarginal': 10., + 'matchingWeight':1., 'cutOffTolerance': 2., 'samplerPivot':QS(ks, "UNIFORM", 2.), 'greedyTol':1e-3, 'errorEstimatorKind':"INTERPOLATORY", 'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.), 'samplerMarginal':QS(ts, "UNIFORM")} algo = RIGP approx = algo([0], solver, mu0 = [k0, t0], approx_state = True, - approxParameters = params, verbosity = 20) + approxParameters = params, verbosity = 10) if len(method) == 2: approx.setupApprox() else: approx.setupApprox("LAST") + print("--- Approximant ---") approx.plotApprox(mu, name = 'u_app') approx.plotHF(mu, name = 'u_HF') - approx.plotErr(mu, name = 'err') - approx.plotRes(mu, name = 'res') + approx.plotErr(mu, name = 'err_app') + approx.plotRes(mu, name = 'res_app') normErr = approx.normErr(mu)[0] normSol = approx.normHF(mu)[0] normRes = approx.normRes(mu)[0] normRHS = approx.normRHS(mu)[0] - print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format( normSol, normErr, normErr / normSol)) - print("RHSNorm:\t{:.5e}\nRes:\t{:.5e}\nResRel:\t{:.5e}".format( + print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) + print("--- Closest snapshot ---") + eng = SES(solver, verbosity = 0) + eng.nsamples = approx.samplingEngine.nsamplesCoalesced + eng.mus = approx.samplingEngine.musCoalesced + eng.samples = approx.samplingEngine.samples_fullCoalesced + paramsNN = {'S':eng.nsamples, 'POD':False, + 'sampler':QS([[ks[0], ts[0]], [ks[-1], ts[-1]]], "UNIFORM")} + approxNN = NN(solver, mu0 = [k0, t0], approx_state = True, + approxParameters = paramsNN, verbosity = 0) + approxNN.setSamples(eng) + approxNN.plotApprox(mu, name = 'u_close') + approxNN.plotHF(mu, name = 'u_HF') + approxNN.plotErr(mu, name = 'err_close') + approxNN.plotRes(mu, name = 'res_close') + normErr = approxNN.normErr(mu)[0] + normSol = approxNN.normHF(mu)[0] + normRes = approxNN.normRes(mu)[0] + normRHS = approxNN.normRHS(mu)[0] + print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format( + normSol, normErr, normErr / normSol)) + print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format( + normRHS, normRes, normRes / normRHS)) + verb = approx.verbosity approx.verbosity = 0 tspace = np.linspace(ts[0], ts[-1], 100) for j, t in enumerate(tspace): pls = approx.getPoles([None, t]) pls[np.abs(np.imag(pls ** 2.)) > 1e-5] = np.nan if j == 0: poles = np.empty((len(tspace), len(pls))) poles[j] = np.real(pls) approx.verbosity = verb - plt.figure(figsize = (12, 5)) - plt.subplot(121) - plt.plot(poles, tspace) - plt.ylim(ts) - plt.xlabel('mu_1') - plt.ylabel('mu_2') - plt.grid() - plt.subplot(122) - plt.plot(poles, tspace) - plt.xlim(ks) - plt.ylim(ts) - plt.xlabel('mu_1') - plt.ylabel('mu_2') - plt.grid() + fighandles += [plt.figure(figsize = (12, 5))] + ax1 = fighandles[-1].add_subplot(1, 2, 1) + ax2 = fighandles[-1].add_subplot(1, 2, 2) + ax1.plot(poles, tspace) + ax1.set_ylim(ts) + ax1.set_xlabel('mu_1') + ax1.set_ylabel('mu_2') + ax1.grid() + ax2.plot(poles, tspace) + for mm in approx.musMarginal: + ax2.plot(ks, [mm[0, 0]] * 2, 'k--', linewidth = 1) + ax2.set_xlim(ks) + ax2.set_ylim(ts) + ax2.set_xlabel('mu_1') + ax2.set_ylabel('mu_2') + ax2.grid() plt.show() - plt.close() - poles, residues = approx.getResidues([None, mu[1]]) - for pol, res in zip(poles, residues.T): - if (np.abs(np.imag(pol ** 2.)) < 1e-5 and np.real(pol) >= ks[0] - and np.real(pol) <= ks[-1]): - solver.plot(res) - print("pole = {:.5e}".format(pol)) - print("\n") diff --git a/examples/4_funnel_output/funnel_output.py b/examples/4_funnel_output/funnel_output.py index 23b01d6..b41419f 100644 --- a/examples/4_funnel_output/funnel_output.py +++ b/examples/4_funnel_output/funnel_output.py @@ -1,41 +1,58 @@ import numpy as np from funnel_output_engine import FunnelOutputEngine as engine -from rrompy.reduction_methods.standard import RationalInterpolant as RI +from rrompy.sampling.standard import SamplingEngineStandard as SES +from rrompy.reduction_methods.standard import (NearestNeighbor as NN, + RationalInterpolant as RI) from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG from rrompy.parameter.parameter_sampling import QuadratureSampler as QS ks = [5., 10.] k0, n = np.mean(ks), 50 solver = engine(k0, n) k = 6.5 for method in ["RI", "RI_OUTPUT", "RI_GREEDY", "RI_GREEDY_OUTPUT"]: print("Testing {} method".format(method)) if "GREEDY" not in method: params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV", 'sampler':QS(ks, "CHEBYSHEV")} algo = RI if "GREEDY" in method: params = {'S':2, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-1, 'maxIter':25, 'sampler':QS(ks, "UNIFORM"), 'errorEstimatorKind':"LOOK_AHEAD_OUTPUT"} algo = RIG approx = algo(solver, mu0 = k0, approx_state = method[-7 :] != "_OUTPUT", approxParameters = params, verbosity = 5) if "GREEDY" not in method: approx.setupApprox() else: approx.setupApprox("LAST") + print("--- Approximant ---") approx.plotApprox(k, name = 'u_app') approx.plotHF(k, name = 'u_HF') - approx.plotErr(k, name = 'err') + approx.plotErr(k, name = 'err_app') normErr = approx.normErr(k)[0] normSol = approx.normHF(k)[0] - print("SolNorm:\t{:.5e}\nErr:\t{:.5e}\nErrRel:\t{:.5e}".format( + print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format( normSol, normErr, normErr / normSol)) + + print("--- Closest snapshot ---") + eng = SES(solver, verbosity = 0) + approxNN = NN(solver, mu0 = k0, approx_state = method[-7 :] != "_OUTPUT", + approxParameters = {'S':approx.S, 'POD':True}, verbosity = 0) + approxNN.setSamples(approx.samplingEngine) + approxNN.plotApprox(k, name = 'u_close') + approxNN.plotHF(k, name = 'u_HF') + approxNN.plotErr(k, name = 'err_close') + normErr = approxNN.normErr(k)[0] + normSol = approxNN.normHF(k)[0] + print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format( + normSol, normErr, normErr / normSol)) + print("Poles:\n{}".format(approx.getPoles())) print("\n")