Page MenuHomec4science
Files F61059219

mhd.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 06:51

mhd.py
View Options

import numpy as np
import matplotlib.pyplot as plt
from mhd_engine import MHDEngine as engine
from rrompy.reduction_methods import (RationalInterpolant as RI,
RationalInterpolantGreedy as RIG)
from rrompy.parameter.parameter_sampling import (FFTSampler as FFTS,
QuadratureCircleSampler as QCS,
QuadratureBoxSampler as QBS)
ks = [-.35 + .5j, 0. + .5j]
k0 = np.mean(ks)
solver = engine(5)
kEffDelta = .1 * (ks[1] - ks[0])
kEff = np.real([ks[0] - kEffDelta, ks[1] + kEffDelta])
iEff = kEff - .5 * np.sum(np.real(ks)) + np.imag(ks[0])
nPoles = 50
polesEx = solver.getPolesExact(nPoles, k0)
for corrector in [False, True]:
for method in ["FFT", "BOX", "GREEDY"]:
print("Testing {} method with{} corrector".format(method,
"out" * (not corrector)))
if method == "FFT":
params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL",
'sampler':FFTS(ks)}
algo = RI
if method == "BOX":
params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL",
'sampler':QBS(ks)}
algo = RI
if method == "GREEDY":
params = {'S':30, 'POD':True, 'greedyTol':1e-2,
'polybasis':"MONOMIAL", 'sampler':QCS(ks),
'errorEstimatorKind':"LOOK_AHEAD", 'nTestPoints':10000,
'trainSetGenerator':FFTS(ks)}
algo = RIG
params['correctorForce'] = corrector
approx = algo(solver, mu0 = k0, approx_state = True,
approxParameters = params, verbosity = 10)
approx.setupApprox()
poles, residues = approx.getResidues()
inRange = np.logical_and(
np.logical_and(np.real(poles) >= kEff[0], np.real(poles) <= kEff[1]),
np.logical_and(np.imag(poles) >= iEff[0], np.imag(poles) <= iEff[1]))
polesEff = poles[inRange]
resNormEff = np.linalg.norm(residues, axis = 1)[inRange]
rLm = np.min(np.log(resNormEff))
rLmM = np.max(np.log(resNormEff)) - rLm
fig = plt.figure(figsize = (10, 10))
ax = fig.add_subplot(1, 1, 1)
if method == "GREEDY":
ax.plot(approx.muTest.re.data.flatten(),
approx.muTest.im.data.flatten(), 'k,', alpha = 0.25)
for pl, rN in zip(polesEff, resNormEff):
if corrector:
alpha = .35 + .4 * (np.log(rN) - rLm) / rLmM
else:
alpha = .1 + .65 * (np.log(rN) - rLm) / rLmM
ax.annotate("{0:.0e}".format(rN), (np.real(pl), np.imag(pl)),
alpha = alpha)
ax.plot(np.real(pl), np.imag(pl), 'r+', alpha = alpha + .25)
ax.plot(approx.mus.re.data.flatten(),
approx.mus.im.data.flatten(), 'k.')
ax.plot(np.real(polesEx), np.imag(polesEx), 'bx')
ax.set_xlim(kEff)
ax.set_ylim(iEff)
ax.grid()
plt.tight_layout()
plt.show()
print("\n")

Event Timeline

c4science · Help