Page MenuHomec4science
Files F60783427

sector_angle.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 14:02

sector_angle.py
View Options

import numpy as np
import matplotlib.pyplot as plt
from sector_angle_engine import SectorAngleEngine as engine
from rrompy.reduction_methods import (NearestNeighbor as NN,
RationalInterpolantPivoted as RIP,
RationalInterpolantGreedyPivoted as RIGP)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
EmptySampler as ES)
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 = {'S':20, "paramsMarginal":{"MMarginal": 3}, 'SMarginal':11,
'POD':True, 'polybasis':"CHEBYSHEV",
'polybasisMarginal':"MONOMIAL_GAUSSIAN",
'radialDirectionalWeightsMarginal': 100.,
'matchingWeight':1., 'samplerPivot':QS(ks, "CHEBYSHEV", 2.),
'samplerMarginal':QS(ts, "UNIFORM")}
algo = RIP
if method == "RI_GREEDY":
params = {'S':10, "paramsMarginal":{"MMarginal": 3}, 'SMarginal':11,
'POD':True, 'polybasis':"LEGENDRE",
'polybasisMarginal':"MONOMIAL_GAUSSIAN",
'radialDirectionalWeightsMarginal': 100.,
'matchingWeight':1., 'samplerPivot':QS(ks, "UNIFORM", 2.),
'greedyTol':1e-3, 'errorEstimatorKind':"LOOK_AHEAD_RES",
'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.),
'samplerMarginal':QS(ts, "UNIFORM")}
algo = RIGP
approx = algo([0], solver, mu0 = [k0, t0], approx_state = True,
approxParameters = params, verbosity = 10,
storeAllSamples = True)
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_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_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 ---")
paramsNN = {'S':len(approx.mus), 'POD':True, 'sampler':ES()}
approxNN = NN(solver, mu0 = [k0, t0], approx_state = True,
approxParameters = paramsNN, verbosity = 0)
approxNN.setSamples(approx.storedSamplesFilenames)
approx.purgeStoredSamples()
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
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()
print("\n")

Event Timeline

c4science · Help