Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/examples/6_boundary_value_problem_1D/boundary_value_problem_1D.py
ca25fdbab062develop

Last Change
Raw File

boundary_value_problem_1D.py
View Options

import warnings
import numpy as np
import matplotlib.pyplot as plt
from boundary_value_problem_1D_engine import BVP1DEngine as engine, BVP1DPoles
from rrompy.reduction_methods import (RationalInterpolant as RI,
RationalInterpolantGreedy as RIG)
from rrompy.parameter.parameter_sampling import (QuadratureBoxSampler as QBS,
QuadratureSampler as QS)
from rrompy.utilities.numerical import potential
ks, shift = [0., 5e3], 200.j
ksRatio = 2. * np.abs(shift) / np.abs(ks[1] - ks[0])
k0, gamma, n = np.mean(ks), .1, 10000
solver = engine(gamma, n, shift)
kEffMult = .1
kEff = np.real([ks[0] - kEffMult * (ks[1] - ks[0]),
ks[1] + kEffMult * (ks[1] - ks[0])])
methods = ["25", "30", "35", "10+", "15+"]
for kind in ["SEGMENT", "CLOUD"]:
print("Testing sampling kind {}...".format(kind))
errAbs, errRel = [], []
for method in methods:
print("S = {}".format(method))
if len(method) == 2:
if kind == "SEGMENT":
sampler = QS(ks, "CHEBYSHEV")
else: #if kind == "CLOUD":
sampler = QBS(ks, "CHEBYSHEV", ksRatio)
params = {'S':int(method), 'POD':True,
'polybasis':"CHEBYSHEV", 'sampler':sampler}
algo = RI
if method[-1] == "+":
if kind == "SEGMENT":
sampler = QS(ks, "UNIFORM")
samplerTrainSet = QS(ks, "CHEBYSHEV")
else: #if kind == "CLOUD":
sampler = QBS(ks, "UNIFORM", ksRatio)
samplerTrainSet = QBS(ks, "CHEBYSHEV", ksRatio)
params = {'S':int(method[: -1]), 'POD':True,
'polybasis':"LEGENDRE", 'greedyTol':1e-3,
'sampler':sampler, 'errorEstimatorKind':"LOOK_AHEAD",
'nTestPoints':1000, 'samplerTrainSet':samplerTrainSet}
algo = RIG
approx = algo(solver, mu0 = k0, approxParameters = params,
verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
polesEff = poles[(np.real(poles) >= kEff[0])
* (np.real(poles) <= kEff[1])]
polesEx = BVP1DPoles(gamma, kEff[0], kEff[1], shift)
polesExTop = np.sort(polesEx[(np.real(polesEx) >= np.real(ks[0]))
* (np.real(polesEx) <= np.real(ks[1]))
* (np.imag(polesEx) > - np.imag(shift))])
polesExBot = np.sort(polesEx[(np.real(polesEx) >= np.real(ks[0]))
* (np.real(polesEx) <= np.real(ks[1]))
* (np.imag(polesEx) < - np.imag(shift))])
polesEx = np.append(np.append(polesExTop, np.inf), polesExBot)
polesExErr = np.abs(np.tile(polesEx.reshape(-1, 1), [1, len(poles)])
- poles.reshape(1,-1))
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
if method[-1] == "+":
print("S_eff = {}".format(approx.S))
ax.plot(approx.muTest.re.data.flatten(),
approx.muTest.im.data.flatten(), 'k,', alpha = 0.25)
ax.plot(approx.mus.re.data.flatten(),
approx.mus.im.data.flatten(), 'k.')
ax.plot(np.real(polesEff), np.imag(polesEff), 'r+')
ax.plot(np.real(polesEx), np.imag(polesEx), 'bx')
ax.set_xlim(kEff)
ax.grid()
plt.tight_layout()
plt.show()
errAbs += [np.min(polesExErr, axis = 1)]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
potEx = potential(approx.trainedModel.centerNormalize(polesEx)(0),
sampler.normalFoci()) / sampler.groundPotential()
potEx[potEx < 1.] = 1.
potEx[len(potEx) // 2] = np.nan
errRel += [2. / np.abs(ks[1] - ks[0]) * errAbs[-1] / potEx]
symbols = '.x+dso^v<>'
fig = plt.figure(figsize = plt.figaspect(.5))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
for j in range(len(methods)):
ax1.semilogy(np.real(polesEx), errAbs[j], '{}-'.format(symbols[j]))
ax2.semilogy(np.real(polesEx), errRel[j], '{}-'.format(symbols[j]))
ax1.set_title("Absolute")
ax1.grid()
ax2.set_title("Relative")
ax2.legend(["S = {}".format(m) for m in methods])
ax2.grid()
plt.tight_layout()
plt.show()
print("\n")
c4science · Help