Page MenuHomec4science
Files F92911747

HelmholtzTaylorApproximantsSweep.py
No OneTemporary
Actions

File Metadata

Created
Sun, Nov 24, 18:20

HelmholtzTaylorApproximantsSweep.py
View Options

import numpy as np
from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.hsengines.fenics import HSEngine as HS
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
from rrompy.reduction_methods.base import ParameterSweeper as Sweeper
testNo = 1
z0 = 12 + .25j
npoints = 31
shift = 5
nsets = 3
stride = 2
Emax = stride * (nsets - 1) + shift + 2
if testNo == 1:
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 10)
params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True}
ktars = np.linspace(7, 16, npoints)
filenamebase = '../data/output/HelmholtzBubbleTaylor'
elif testNo == 2:
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10)
params = {'Emax':Emax, 'sampleType':'KRYLOV', 'POD':True}
ktars = np.linspace(11, 13, npoints)
filenamebase = '../data/output/HelmholtzBoxTaylor'
plotter = HS(solver.V)
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift,
'E':stride*i+shift+1}
paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+1}
appPade = Pade(solver, plotter, mu0 = z0, approxParameters = params)
appRB = RB(solver, plotter, mu0 = z0, approxParameters = params)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = appPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase + 'PadeFE.dat', outputs = 'ALL')
sweeper.ROMEngine = appRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase + 'RBFE.dat', outputs = 'ALL')
####################
from matplotlib import pyplot as plt
for i in range(nsets):
nDerivatives = stride*i+shift+1
PadeOutput = sweeper.read(filenamePade, {'E':nDerivatives},
['muRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'E':nDerivatives},
['muRe', 'AppNorm', 'ErrNorm'])
ktarsF = PadeOutput['muRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['muRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['ErrNorm']
RBktarsF = RBOutput['muRe']
RBNormF = RBOutput['AppNorm']
RBErrorF = RBOutput['ErrNorm']
plt.figure()
plt.semilogy(ktarsF, solNormF, 'k-', label='Sol norm')
plt.semilogy(PadektarsF, PadeNormF, 'b.--',
label='Pade'' norm, E = {}'.format(nDerivatives))
plt.semilogy(RBktarsF, RBNormF, 'g.--',
label='RB norm, E = {}'.format(nDerivatives))
plt.legend()
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(PadektarsF, PadeErrorF, 'b',
label='Pade'' error, E = {}'.format(nDerivatives))
plt.semilogy(RBktarsF, RBErrorF, 'g',
label='RB error, E = {}'.format(nDerivatives))
plt.legend()
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(ktarsF, PadeErrorF / solNormF, 'b',
label='Pade'' relative error, E = {}'.format(nDerivatives))
plt.semilogy(RBktarsF, RBErrorF / solNormF, 'g',
label='RB relative error, E = {}'.format(nDerivatives))
plt.legend()
plt.grid()
plt.show()
plt.close()

Event Timeline

c4science · Help