Page MenuHomec4science
Files F78332730

HelmholtzTaylorApproximantsSweep.py
No OneTemporary
Actions

File Metadata

Created
Mon, Aug 19, 21:24

HelmholtzTaylorApproximantsSweep.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Example homogeneous Dirichlet forcing wave SWEEP
from __future__ import print_function
import numpy as np
from context import FenicsHelmholtzEngine as HFEngine
from context import FenicsHSEngine as HSEngine
from context import ROMApproximantTaylorPade as Pade
from context import ROMApproximantTaylorRB as RB
from context import ROMApproximantSweeper as Sweeper
z0 = 12 + .5j
npoints = 31
ktars = np.linspace(7, 16, npoints)
from Fenics_snippets import SquareHomogeneousBubble
bd, msh, f = SquareHomogeneousBubble(kappa = 12**.5, theta = np.pi / 3, n = 10)
solver = HFEngine(mesh = msh, wavenumber = z0**.5, forcingTerm = f,
FEDegree = 3, DirichletBoundary = bd)
plotter = HSEngine(solver.V)
shift = 5
nsets = 3
stride = 2
Emax = stride * (nsets - 1) + shift + 2
params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True}
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+2}
appPade = Pade(solver, plotter, k0 = z0, w = np.real(z0**.5),
approxParameters = params)
appRB = RB(solver, plotter, k0 = z0, w = np.real(z0**.5),
approxParameters = params)
sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = appPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('../Data/HelmholtzBubbleTaylorPadeFE.dat',
outputs = 'ALL')
sweeper.ROMEngine = appRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep('../Data/HelmholtzBubbleTaylorRBFE.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},
['kRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'E':nDerivatives},
['kRe', 'AppNorm', 'ErrNorm'])
ktarsF = PadeOutput['kRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['kRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['ErrNorm']
RBktarsF = RBOutput['kRe']
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