Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/9_active_remeshing/active_remeshing.py
63e560ed3acamaster

Last Change
Raw File

active_remeshing.py
View Options

import numpy as np
from pickle import load
from matplotlib import pyplot as plt
from active_remeshing_engine import ActiveRemeshingEngine
from rrompy.reduction_methods import (
RationalInterpolantGreedyPivotedNoMatch as RIGPNM,
RationalInterpolantGreedyPivotedPoleMatch as RIGPG)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
SparseGridSampler as SGS)
zs, ts = [0., 100.], [0., .5]
z0, t0, n = np.mean(zs), np.mean(ts), 150
solver = ActiveRemeshingEngine(z0, t0, n)
solver.cutOffPolesRMinRel, solver.cutOffPolesRMaxRel = -2.5, 2.5
solver.cutOffPolesIMin, solver.cutOffPolesIMax = -.01, .01
mus = [[z0, ts[0]], [z0, ts[1]]]
for mu in mus:
u = solver.solve(mu, return_state = True)[0]
Y = solver.applyC(u, mu)[0]
_ = solver.plot(u, what = "REAL", name = "u(z={}, t={})".format(*mu),
is_state = True, figsize = (12, 4))
print("Y(z={}, t={}) = {} (solution average)".format(*mu, np.real(Y)))
fighandles = []
with open("./active_remeshing_hf_samples.pkl", "rb") as f:
zspace, tspace, Yex = load(f)
# zspace, tspace = np.linspace(*zs, 200), np.linspace(*ts, 50)
# Yex = [[solver.solve([z, t]) for t in tspace] for z in zspace]
# (from a ~2h45m simulation on one node of the EPFL Helvetios cluster)
radius2Err = np.mean(np.abs(Yex) ** 2.)
YCmin, YCmax = np.quantile(Yex, .05), np.quantile(Yex, .95)
YexC = np.clip(Yex, YCmin, YCmax)
approx = []
for match in range(2):
params = {"POD": True, "S": 5, "greedyTol": 1e-4, "nTestPoints": 500,
"polybasis": "LEGENDRE", "samplerTrainSet": QS(zs, "UNIFORM"),
"samplerPivot": QS(zs, "CHEBYSHEV"), "SMarginal": 5,
"samplerMarginal": SGS(ts),
"errorEstimatorKind": "LOOK_AHEAD_OUTPUT"}
if match:
print("\nTesting output-based matching.")
params["matchingWeight"] = 1.
params["matchingShared"] = .75
params["polybasisMarginal"] = "PIECEWISE_LINEAR_UNIFORM"
algo = RIGPG
else:
print("\nTesting matching-free approach.")
algo = RIGPNM
approx += [algo([0], solver, mu0 = [z0, t0], approxParameters = params,
verbosity = 5)]
if match:
approx[match].setTrainedModel(approx[0])
else:
approx[match].setupApprox()
verb = approx[match].verbosity
verbTM = approx[match].trainedModel.verbosity
approx[match].verbosity, approx[match].trainedModel.verbosity = 0, 0
for j, t in enumerate(tspace):
out = approx[match].getApprox(np.pad(zspace.reshape(-1, 1),
[(0, 0), (0, 1)], "constant",
constant_values = t))
pls = approx[match].getPoles([None, t])
pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
if j == 0:
Ys = np.empty((len(zspace), len(tspace)))
poles = np.empty((len(tspace), len(pls)))
Ys[:, j] = out.re.data
if len(pls) > poles.shape[1]:
poles = np.pad(poles, [(0, 0), (0, len(pls) - poles.shape[1])],
"constant", constant_values = np.nan)
poles[j, : len(pls)] = np.real(pls)
approx[match].verbosity = verb
approx[match].trainedModel.verbosity = verbTM
YsC = np.clip(Ys, YCmin, YCmax)
err = (np.abs(Yex - YsC) / (np.abs(Yex) ** 2. + radius2Err) ** .5
/ (np.abs(Ys) ** 2. + radius2Err) ** .5)
fighandles += [plt.figure(figsize = (15, 5))]
ax1 = fighandles[-1].add_subplot(1, 2, 1)
ax2 = fighandles[-1].add_subplot(1, 2, 2)
if match:
ax1.plot(poles, tspace)
else:
ax1.plot(poles, tspace, "k.")
ax1.set_ylim(ts)
ax1.set_xlabel("z")
ax1.set_ylabel("t")
ax1.grid()
if match:
ax2.plot(poles, tspace)
else:
ax2.plot(poles, tspace, "k.")
for mm in approx[match].musMarginal:
ax2.plot(zs, [mm[0, 0]] * 2, "k--", linewidth = 1)
ax2.set_xlim(zs)
ax2.set_ylim(ts)
ax2.set_xlabel("z")
ax2.set_ylabel("t")
ax2.grid()
plt.show()
print("Approximate poles")
fighandles += [plt.figure(figsize = (15, 5))]
ax1 = fighandles[-1].add_subplot(1, 2, 1)
ax2 = fighandles[-1].add_subplot(1, 2, 2)
p = ax1.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
YsC, vmin = YCmin, vmax = YCmax,
levels = np.linspace(YCmin, YCmax, 50))
plt.colorbar(p, ax = ax1)
ax1.set_xlabel("z")
ax1.set_ylabel("t")
ax1.grid()
p = ax2.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
YexC, vmin = YCmin, vmax = YCmax,
levels = np.linspace(YCmin, YCmax, 50))
ax2.set_xlabel("z")
ax2.set_ylabel("t")
ax2.grid()
plt.colorbar(p, ax = ax2)
plt.show()
print("Approximate and exact output\n")
fighandles += [plt.figure(figsize = (9, 6))]
ax1 = fighandles[-1].add_subplot(1, 1, 1)
p = ax1.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
np.log10(err), vmin = -10, vmax = 0,
levels = np.linspace(-10, 0, 50))
plt.colorbar(p, ax = ax1)
ax1.set_xlabel("z")
ax1.set_ylabel("t")
ax1.grid()
plt.show()
print("Output error (log-chordal)\n")
c4science · Help