Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python/HelmholtzRBLagrangeApproximant.py
7e2d435a8d1bmaster

Last Change
Raw File

HelmholtzRBLagrangeApproximant.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
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 ROMApproximantLagrangeRB as RB
PI = np.pi
testNo = 3
if testNo == 1:
params = {'S':5, 'nodesType':'CHEBYSHEV', 'POD':True}
nu = 12**.5
theta = PI / 3
ztar = 11
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 = 40
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.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = 'all', DirichletDatum = 0)
plotter = HSEngine.FenicsHSEngine(solver.V)
approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [10 + .5j, 14 + .5j], w = np.real(nu),
approxParameters = params)
approx.plotApp(ztar, name = 'u_RB')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(appErr, solNorm, appErr/solNorm)
print(approx.getPoles(True))
############
elif testNo == 2:
def Dboundary(x, on_boundary):
return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
A = 10
kappa = 4
theta = PI * 90 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - A * sp.exp(1.j * phiex)
nx = ny = 40
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
params = {'S':30, 'nodesType':'CHEBYSHEV', 'POD':True}
solver = HFEngine.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm)
plotter = HSEngine.FenicsHSEngine(solver.V)
approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [0, 8],
approxParameters = params,
plotSnapshots = False)
ktar = 4.5
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
print(appErr, solNorm, appErr/solNorm)
print(approx.getPoles(True))
############
elif testNo == 3:
def Dboundary(x, on_boundary):
return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
A = 10
kappa = 4
theta = PI * 90 / 180
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - A * sp.exp(1.j * phiex)
nx = ny = 40
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
params = {'S':30, 'nodesType':'CHEBYSHEV', 'POD':True}
solver = HFEngine.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
constraintType = 'IDENTITY')
plotter = HSEngine.FenicsHSAugmentedEngine(solver.V, 2)
approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [0, 8],
approxParameters = params,
plotSnapshots = False)
ktar = 4.5
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
appErr, solNorm = approx.approxError(ktar, kappa), approx.HFNorm(ktar, kappa)
print(appErr, solNorm, np.divide(appErr, solNorm))
print(approx.getPoles(True))
c4science · Help