Page MenuHomec4science
Files F61137682

pod.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 18:54

pod.py
View Options

import numpy as np
from airfoil_engine import AirfoilScatteringEngine
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
verb = 100
sol = "single"
sol = "sweep"
algo = "RI"
algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
Nsweep = 100
k0s = [x * 2 * np.pi / 340 for x in [1.0e2, 5.0e2]]
k0 = np.mean(np.power(k0s, 2.)) ** .5
ktar = k0s[0] + (k0s[1] - k0s[0]) * .7
params = {'N':29, 'M':29, 'R':30, 'S':30, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0s, "CHEBYSHEV"), 'robustTol':1e-14}
theta = - 45. * np.pi / 180
solver = AirfoilScatteringEngine(k0, theta, verbosity = verb,
degree_threshold = 8)
if algo == "RI":
params.pop('R')
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
params.pop('N')
params.pop('M')
params.pop('polybasis')
params.pop('robustTol')
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.setupApprox()
if sol == "single":
# approx.outParaviewTimeDomainSamples(filename = "out/outSamples",
# forceNewFile = False, folders = True)
approx.outParaviewTimeDomainApprox(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTApp{}".format(ktar),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainHF(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTHF{}".format(ktar),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainErr(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTErr{}".format(ktar),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainRes(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTRes{}".format(ktar),
forceNewFile = False, folder = True)
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('Poles:', approx.getPoles())
if sol == "sweep":
k0s = np.linspace(k0s[0], k0s[1], Nsweep)
kl, kr = min(k0s), max(k0s)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
err[j] = approx.normErr(k0s[j]) / norm[j]
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float),
'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()

Event Timeline

c4science · Help