Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python/HelmholtzMembraneTaylor.py
4c3d184239eemaster

Last Change
Raw File

HelmholtzMembraneTaylor.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import fenics as fen
import numpy as np
import sympy as sp
from context import FenicsHelmholtzEngine as HF
from context import FenicsHSEngine as HS
from context import ROMApproximantTaylorPade as TP
from context import ROMApproximantSweeper as Sweeper
from matplotlib import pyplot as plt
plt.jet()
####################
test = "solve"
#test = "Taylor"
test = "TaylorSweep"
ttype = "simple"
plotSamples = 'ABS'
#plotSamples = []
z0 = 100 + 1.j
ktar = 97.5**.5
ktars = np.sort(np.concatenate((np.linspace(78, 122, 101), np.linspace(100.9, 101.6, 5))))
ktars = np.array([z0])
####################
def boundaryNeumann(x, on_boundary):
return on_boundary and x[1] > .25 and x[0] > 0.995 and x[0] < 1.005
PI = np.pi
nu = 10
x0, y0 = .5, .5
A, B = 1., 0
Rr, Ri = .1, .1
x, y = sp.symbols('x[0] x[1]', real=True)
fex = (A * sp.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Rr**2)
+ 1.j * B * sp.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Ri**2))
meshname = '../Data/crack_coarse.xml'
#meshname = '../Data/crack_fine.xml'
mesh = fen.Mesh(meshname)
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
plt.figure()
fen.plot(mesh, title = 'mesh')
#plt.gcf().savefig('./Membrane/mesh.eps', format='eps', dpi=1000)
solver = HF(mesh = mesh, wavenumber = ktar, forcingTerm = forcingTerm, FEDegree = 3,
NeumannBoundary = boundaryNeumann, DirichletBoundary = 'rest')
plotter = HS(solver.V)
def crack_disp(u):
uRe = fen.Function(solver.V)
uIm = fen.Function(solver.V)
uRe.vector()[:] = np.array(np.real(u), dtype = float)
uIm.vector()[:] = np.array(np.imag(u), dtype = float)
return fen.assemble((uRe * uRe + uIm * uIm) * solver.ds(0)) ** .5
solver.functional = crack_disp
if test == "solve":
uh = solver.solve()
plotter.plot(uh, what = 'REAL')
#plt.gcf().savefig('./Membrane/u0.eps', format='eps', dpi=1000)
print(solver.functional(uh))
elif test == "Taylor":
params = {'N':6, 'M':5, 'E':6, 'sampleType':'Arnoldi', 'POD':True}
approxPade = TP(solver, plotter, k0 = z0, w = nu, plotDer = plotSamples,
approxParameters = params)
PadeErr, solNorm = approxPade.approxError(ktar), approxPade.HFNorm(ktar)
print('SolNorm:\t{}\nErrPade:\t{}\nErrRelPade:\t{}'.format(
solNorm, PadeErr, np.divide(PadeErr, solNorm)))
print('\nPoles Pade'':')
print(approxPade.getPoles(True))
uHF = approxPade.getHF(ktar)
plotter.plot(uHF, name = 'u_ex')
approxPade.plotApp(ktar, name = 'u_Pade''')
approxPade.plotErr(ktar, name = 'errPade''')
elif test == "TaylorSweep":
shift = 2
nsets = 3
stride = 3
Emax = stride * (nsets - 1) + shift + 1
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift,
'E':stride*i+shift+1}
approxPade = TP(solver, plotter, k0 = z0, approxParameters = params)
sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('../Data/HelmholtzMembraneTP.dat', outputs = 'ALL')
for i in range(nsets):
val = stride*i+shift + 1
PadeOutput = sweeper.read(filenamePade, {'E':val},
['kRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
ktarsF = PadeOutput['kRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['kRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['ErrNorm']
plt.figure(figsize = (10, 5))
plt.semilogy(ktarsF, solNormF, 'k-', label='Sol norm')
plt.semilogy(PadektarsF, PadeNormF, 'b--',
label='Pade'' norm, E = {0}'.format(val))
plt.legend()
plt.grid()
p = plt.legend(loc = 'upper left')
#plt.savefig('./Membrane/normA' + str(i) + '.eps', format='eps', dpi=1000)
plt.figure(figsize = (10, 5))
plt.semilogy(PadektarsF, PadeErrorF, 'b',
label='Pade'' error, E = {0}'.format(val))
plt.legend()
plt.grid()
p = plt.legend(loc = 'upper left')
#plt.savefig('./Membrane/errorA' + str(i) + '.eps', format='eps', dpi=1000)
plt.figure(figsize = (10, 5))
plt.semilogy(ktarsF, np.divide(PadeErrorF, solNormF), 'b',
label='Pade'' relative error, E = {0}'.format(val))
plt.legend()
plt.grid()
p = plt.legend(loc = 'upper left')
#plt.savefig('./Membrane/errorAR' + str(i) + '.eps', format='eps', dpi=1000)
c4science · Help