Page MenuHomec4science
Files F62283066

double_slit.py
No OneTemporary
Actions

File Metadata

Created
Sun, May 12, 04:08

double_slit.py
View Options

import numpy as np
from double_slit_engine import DoubleSlitEngine as engine
from rrompy.reduction_methods import (NearestNeighbor as NN,
RationalInterpolant as RI,
RationalInterpolantGreedy as RIG)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
EmptySampler as ES)
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 = {'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_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_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
normSol, normErr, normErr / normSol))
print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format(
normRHS, normRes, normRes / normRHS))
print("--- Closest snapshot ---")
approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0,
approxParameters = {'S':approx.S, 'POD':True,
'sampler':ES()})
approxNN.setSamples(approx.storeSamples())
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):
solver.plot(res)
print("pole = {:.5e}".format(pol))
print("\n")

Event Timeline

c4science · Help