Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/Python/HelmholtzLagrangeApproximantsSweep.py
a2de4f5207cffeature_greedy_coarseni…

Last Change
Raw File

HelmholtzLagrangeApproximantsSweep.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Example homogeneous Dirichlet forcing wave SWEEP
from __future__ import print_function
import fenics as fen
import numpy as np
import sympy as sp
from context import FenicsHelmholtzEngine as HFEngine
from context import FenicsHSEngine as HSEngine
from context import ROMApproximantLagrangePade as Pade
from context import ROMApproximantLagrangeRB as RB
from context import ROMApproximantSweeper as Sweeper
PI = np.pi
nu = 12**.5
theta = PI / 3
npoints = 31
ktars = np.linspace(0, 21, npoints)
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
nx = ny = 10
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
solver = HFEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = 'all', DirichletDatum = 0)
plotter = HSEngine(solver.V)
shift = 3
nsets = 5
stride = 3
Smax = stride * (nsets - 1) + shift + 2
paramsPade = {'S':Smax, 'polyBasis':'CHEBYSHEV', 'POD':True}
paramsRB = {'S':Smax, 'nodesType':'CHEBYSHEV', 'POD':True}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N': stride * i + shift + 1, 'M': stride * i + shift,
'S': stride * i + shift + 2}
paramsSetsRB[i] = {'R': stride * i + shift + 1, 'S': stride * i + shift + 2}
appPade = Pade(solver, plotter, ks = [4 + .5j, 14 + .5j],
w = nu, approxParameters = paramsPade)
appRB = RB(solver, plotter, ks = [4 + .5j, 14 + .5j],
w = nu, approxParameters = paramsRB)
sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = appPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('../Data/HelmholtzBubbleLagrangePade.dat', outputs = 'ALL')
sweeper.ROMEngine = appRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep('../Data/HelmholtzBubbleLagrangeRB.dat', outputs = 'ALL')
####################
from matplotlib import pyplot as plt
plt.jet()
for i in range(nsets):
nSamples = stride*i+shift+2
PadeOutput = sweeper.read(filenamePade, {'S':nSamples},
['kRe', 'HFNorm', 'AppNorm', 'AppError'])
RBOutput = sweeper.read(filenameRB, {'S':nSamples},
['kRe', 'AppNorm', 'AppError'])
ktarsF = PadeOutput['kRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['kRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['AppError']
RBktarsF = RBOutput['kRe']
RBNormF = RBOutput['AppNorm']
RBErrorF = RBOutput['AppError']
plt.figure()
plt.semilogy(ktarsF, solNormF, 'k-', label='Sol norm')
plt.semilogy(PadektarsF, PadeNormF, 'b.--', label='Pade'' norm, S = {}'.format(nSamples))
plt.semilogy(RBktarsF, RBNormF, 'g.--', label='RB norm, S = {}'.format(nSamples))
plt.legend()
plt.grid()
plt.figure()
plt.semilogy(PadektarsF, PadeErrorF, 'b', label='Pade'' error, S = {}'.format(nSamples))
plt.semilogy(RBktarsF, RBErrorF, 'g', label='RB error, S = {}'.format(nSamples))
plt.legend()
plt.grid()
plt.figure()
plt.semilogy(ktarsF, PadeErrorF / solNormF, 'b', label='Pade'' relative error, S = {}'.format(nSamples))
plt.semilogy(RBktarsF, RBErrorF / solNormF, 'g', label='RB relative error, S = {}'.format(nSamples))
plt.legend()
plt.grid()
c4science · Help