Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FEniCS/HelmholtzRBLagrangeApproximant.py
45608c38c688master

Last Change
Raw File

HelmholtzRBLagrangeApproximant.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import print_function
import numpy as np
from context import FenicsHelmholtzEngine as HFEngine
from context import FenicsHelmholtzScatteringEngine as HFSEngine
from context import FenicsHelmholtzScatteringAugmentedEngine as HFSAEngine
from context import FenicsHSEngine as HSEngine
from context import FenicsHSAugmentedEngine as HSAEngine
from context import ROMApproximantLagrangeRB as RB
testNo = 3
if testNo == 1:
params = {'S':5, 'nodesType':'CHEBYSHEV', 'POD':True}
z0s = [10 + .5j, 14 + .5j]
z0 = np.mean(z0s)
ztar = 11
from FEniCS_snippets import SquareHomogeneousBubble
boundary, mesh, forcingTerm = SquareHomogeneousBubble(kappa = 12 ** .5,
theta = np.pi / 3,
n = 40)
solver = HFEngine(mesh = mesh, wavenumber = z0**.5,
forcingTerm = forcingTerm, FEDegree = 3,
DirichletBoundary = boundary, DirichletDatum = 0)
plotter = HSEngine(solver.V)
approx = RB(solver, plotter, ks = z0s, w = np.real(z0**.5),
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(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles(True))
############
elif testNo == 2:
params = {'S':10, 'nodesType':'CHEBYSHEV', 'POD':True}
z0s = np.power([3.85 + .15j, 4.15 + .15j], 2.)
z0 = np.mean(z0s)
ztar = 4 ** 2.
from FEniCS_snippets import SquareTransmissionDirichlet
boundary, mesh, n, u0 = SquareTransmissionDirichlet(nT = 2, nB = 1,
theta = np.pi * 45 / 180,
kappa = 4., n = 50)
solver = HFEngine(mesh = mesh, wavenumber = z0**.5,
refractionIndex = n, FEDegree = 3,
DirichletBoundary = boundary, DirichletDatum = u0)
plotter = HSEngine(solver.V)
approx = RB(solver, plotter, ks = z0s, w = np.real(z0**.5),
approxParameters = params, plotSnap = 'ALL')
approx.plotApp(ztar, name = 'u_Pade''')
approx.plotHF(ztar, name = 'u_HF')
approx.plotErr(ztar, name = 'err')
appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles(True))
############
elif testNo in [3, 4]:
params = {'S':30, 'nodesType':'CHEBYSHEV', 'POD':True}
k0s = [0, 8]
k0 = np.mean(k0s)
ktar = 4.5
from FEniCS_snippets import SquareScatteringTB
bdrD, bdrN, mesh, forcingTerm = SquareScatteringTB(kappa = 4,
theta = np.pi / 2,
n = 40)
if testNo == 3:
solver = HFSEngine(mesh = mesh, wavenumber = k0, FEDegree = 3,
forcingTerm = forcingTerm, DirichletBoundary = bdrD,
RobinBoundary = bdrN)
plotter = HSEngine(solver.V)
else:
solver = HFSAEngine(mesh = mesh, wavenumber = k0, FEDegree = 3,
forcingTerm = forcingTerm,
DirichletBoundary = bdrD, RobinBoundary = bdrN)
plotter = HSAEngine(solver.V, 2)
approx = RB(solver, plotter, ks = k0s, approxParameters = params)
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(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print('\nPoles Pade'':')
print(approx.getPoles(True))
c4science · Help