Page MenuHomec4science
Files F61022121

active_remeshing.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 01:07

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 (
RationalInterpolantGreedyPivotedGreedyNoMatch as RIGPGNM,
RationalInterpolantGreedyPivotedGreedy 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)
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][0]
_ = solver.plot(u, what = "REAL", name = "u(z={}, t={})".format(*mu),
is_state = True, figsize = (12, 4))
print("Y(z={}, t={}) = {} (solution norm squared)".format(*mu, Y))
fighandles = []
with open("./active_remeshing_hf_samples.pkl", "rb") as f:
zspace, tspace, Yex = load(f)
# zspace = np.linspace(zs[0], zs[-1], 198)
# tspace = np.linspace(ts[0], ts[-1], 50)
# Yex = np.log10([[solver.solve([z, t]) for t in tspace] for z in zspace])
# (from a ~2h45m simulation on one node of the EPFL Helvetios cluster)
for match in [0, 1]:
params = {"POD": True, "S": 5, "SMarginal": 3, "greedyTol": 1e-4,
"nTestPoints": 500, "polybasis": "LEGENDRE",
"maxIterMarginal": 25, "trainSetGenerator": QS(zs, "UNIFORM"),
"samplerPivot":QS(zs, "CHEBYSHEV"), "samplerMarginal":SGS(ts),
"errorEstimatorKind": "LOOK_AHEAD_OUTPUT",
"errorEstimatorKindMarginal": "LOOK_AHEAD_RECOVER"}
if match:
print("\nTesting output-based matching with weight 1.")
params["greedyTolMarginal"] = 1e-2
params["matchingWeight"] = 1.
params["sharedRatio"] = .75
params["polybasisMarginal"] = "PIECEWISE_LINEAR_UNIFORM"
algo = RIGPG
else:
print("\nTesting matching-free approach.")
params["greedyTolMarginal"] = 2e-2
algo = RIGPGNM
approx = algo([0], solver, mu0 = [z0, t0], approxParameters = params,
verbosity = 15)
approx.setupApprox("ALL")
verb, verbTM = approx.verbosity, approx.trainedModel.verbosity
approx.verbosity, approx.trainedModel.verbosity = 0, 0
for j, t in enumerate(tspace):
out = approx.getApprox(np.pad(zspace.reshape(-1, 1), [(0, 0), (0, 1)],
"constant", constant_values = t))
pls = approx.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] = np.log10(out.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.verbosity, approx.trainedModel.verbosity = verb, verbTM
Ymin, Ymax = min(np.min(Ys), np.min(Yex)), max(np.max(Ys), np.max(Yex))
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.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),
Ys, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
levels = np.linspace(Ymin, Ymax, 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),
Yex, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
levels = np.linspace(Ymin, Ymax, 50))
ax2.set_xlabel('z')
ax2.set_ylabel('t')
ax2.grid()
plt.colorbar(p, ax = ax2)
plt.show()
print("Approximate and exact output\n")

Event Timeline

c4science · Help