Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/examples/fenics/HelmholtzTaylorApproximantsSweep.py
6513a8e1212adevelop

Last Change
Raw File

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()
c4science · Help