Page MenuHomec4science
Files F62516974

HelmholtzMembraneTaylor.py
No OneTemporary
Actions

File Metadata

Created
Mon, May 13, 18:17

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)

Event Timeline

c4science · Help