Page MenuHomec4science
Files F73218841

HelmholtzAirfoilTaylor.py
No OneTemporary
Actions

File Metadata

Created
Fri, Jul 19, 17:28

HelmholtzAirfoilTaylor.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 FenicsHelmholtzScatteringEngine as HFS
from context import FenicsHelmholtzScatteringAugmentedEngine as HFSA
from context import FenicsHSEngine as HS
from context import FenicsHSAugmentedEngine as HSA
from context import ROMApproximantTaylorPade as TP
from context import ROMApproximantTaylorRB as TRB
from context import ROMApproximantLagrangePade as LP
from context import ROMApproximantLagrangeRB as LRB
from context import ROMApproximantSweeper as Sweeper
from matplotlib import pyplot as plt
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
####################
test = "solve"
test = "Taylor"
#test = "Lagrange"
test = "TaylorSweep"
#test = "LagrangeSweep"
ttype = "simple"
#ttype = "augmentedI"
#ttype = "augmentedM"
plotSamples = 'ALL'
#plotSamples = []
k0 = 10 + 1.j
kLeft, kRight = 8 + 1.j, 12 + 1.j
ktar = 11
ktars = np.linspace(8, 12, 33) - .5j
ktars = np.linspace(12.125, 12.5, 4) - .5j
####################
PI = np.pi
R = 2
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
kappa = 10
theta = PI * -30 / 180.
mu = 1.1
epsilon = .1
x, y = sp.symbols('x[0] x[1]', real=True)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
u0ex = - sp.exp(1.j * phiex)
checkReal = x**2-x+y**2
rhoroot = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = sp.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = sp.atan(-y/(x**2-x+y**2)) / 2 - PI * sp.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhoroot*sp.cos(phiroot1)+.5)**2.+(rhoroot*sp.sin(phiroot1))**2.)/
((rhoroot*sp.cos(phiroot1)-1.)**2.+(rhoroot*sp.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhoroot*sp.cos(phiroot2)+.5)**2.+(rhoroot*sp.sin(phiroot2))**2.)/
((rhoroot*sp.cos(phiroot2)-1.)**2.+(rhoroot*sp.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + sp.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + sp.sin(PI*kappam2/epsilon) / PI) + .1
mesh = fen.Mesh('../Data/airfoil.xml')
a = fen.Expression(('{0}>=0 ? ({2}>=-{1} ? ({2}<={1} ? {4} : 1) : .1) : '
'({3}>=-{1} ? ({3}<={1} ? {5} : 1) : .1)')\
.format(sp.printing.ccode(checkReal), str(epsilon),
sp.printing.ccode(kappam1),
sp.printing.ccode(kappam2),
sp.printing.ccode(Heps1),
sp.printing.ccode(Heps2)), degree = 4)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))),
sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
if ttype == "simple":
solver = HFS(mesh = mesh, wavenumber = k0, diffusivity = a,
forcingTerm = 0, FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm)
plotter = HS(solver.V)
else:
if ttype[-1] == "I": constraintType = "IDENTITY"
else: constraintType = "MASS"
solver = HFSA(mesh = mesh, wavenumber = k0, diffusivity = a,
forcingTerm = 0, FEDegree = 3, DirichletBoundary = Dboundary,
RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
constraintType = constraintType)
plotter = HSA(solver.V, d = 2)
baseRe, baseIm = DirichletTerm
baseRe = fen.project(fen.Expression(baseRe, degree = 4), solver.V)
baseIm = fen.project(fen.Expression(baseIm, degree = 4), solver.V)
uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
if ttype[: -1] == "augmented":
uinc = [uinc, kappa * uinc]
if test == "solve":
aF = fen.interpolate(a, solver.V)
av = aF.vector()
uh = solver.solve()
print(plotter.norm(uh, kappa))
if ttype == "simple":
uhtot = uh - uinc
else:
uhtot = [x - y for x, y in zip(uh, uinc)]
print(plotter.norm(uhtot, kappa))
plotter.plot(av, split = False, what = 'Real', name = 'a')
plotter.plot(uhtot - uh, what = 'Real', name = 'u_inc')
plotter.plot(uh, what = 'ABS')
plotter.plot(uhtot, what = 'ABS', name = 'u + u_inc')
elif test in ["Taylor", "Lagrange"]:
if test == "Taylor":
params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Krylov', 'POD':False}
parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
approxPade = TP(solver, plotter, k0 = k0, plotDer = plotSamples,
approxParameters = parPade)
approxRB = TRB(solver, plotter, k0 = k0, approxParameters = parRB)
else:
params = {'N':8, 'M':7, 'R':9, 'S':9, 'polyBasis':'CHEBYSHEV',
'nodesType':'CHEBYSHEV','POD':True}
parPade = subdict(params, ['N', 'M', 'S', 'polyBasis', 'POD'])
parRB = subdict(params, ['R', 'S', 'nodesType', 'POD'])
approxPade = LP(solver, plotter, ks = [kLeft, kRight], w = kappa,
plotSnap = plotSamples, approxParameters = parPade)
approxRB = LRB(solver, plotter, ks = [kLeft, kRight], w = kappa,
approxParameters = parRB)
PadeErr, solNorm = approxPade.approxError(ktar), approxPade.HFNorm(ktar)
RBErr = approxRB.approxError(ktar)
print(('SolNorm:\t{}\nErrPade:\t{}\nErrRelPade:\t{}\nErrRB:\t\t{}'
'\nErrRelRB:\t{}').format(solNorm, PadeErr,
np.divide(PadeErr, solNorm), RBErr,
np.divide(RBErr, solNorm)))
print('\nPoles Pade'':')
print(approxPade.getPoles(True))
print('\nPoles RB:')
print(approxRB.getPoles(True))
uHF = approxPade.getHF(ktar)
plotter.plot(uHF, name = 'u_ex')
approxPade.plotApp(ktar, name = 'u_Pade''')
approxPade.plotErr(ktar, name = 'errPade''')
approxRB.plotApp(ktar, name = 'u_RB')
approxRB.plotErr(ktar, name = 'errRB')
elif test in ["TaylorSweep", "LagrangeSweep"]:
if test == "TaylorSweep":
shift = 2
nsets = 4
stride = 3
Emax = stride * (nsets - 1) + shift + 1
params = {'Emax':Emax, 'sampleType':'Krylov', 'POD':False}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift,
'E':stride*i+shift+1}
paramsSetsRB[i] = {'E':stride*i+shift,'R':stride*i+shift+1}
approxPade = TP(solver, plotter, k0 = kappa, approxParameters = params)
approxRB = TRB(solver, plotter, k0 = kappa, approxParameters = params)
else:
kLeft, kRight = 8, 12
shift = 3
nsets = 3
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}
approxPade = LP(solver, plotter, ks = [kLeft, kRight], w = kappa,
approxParameters = paramsPade)
approxRB = LRB(solver, plotter, ks = [kLeft, kRight], w = kappa,
approxParameters = paramsRB)
sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars,
mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade =sweeper.sweep('../Data/HelmholtzAirfoil'+test[:-5]+'PadeFE.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep('../Data/HelmholtzAirfoil'+test[:-5]+'RBFE.dat')
for i in range(nsets):
if test == "TaylorSweep":
val = stride*i+shift
PadeOutput = sweeper.read(filenamePade, {'M':val},
['kRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'E':val},
['kRe', 'AppNorm', 'ErrNorm'])
let = 'E'
else:
val = stride*i+shift+2
PadeOutput = sweeper.read(filenamePade, {'S':val},
['kRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'S':val},
['kRe', 'AppNorm', 'ErrNorm'])
let = 'S'
ktarsF = PadeOutput['kRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['kRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['ErrNorm']
RBktarsF = RBOutput['kRe']
RBNormF = RBOutput['AppNorm']
RBErrorF = RBOutput['ErrNorm']
plt.figure(figsize = (10, 5))
plt.plot(ktarsF, solNormF, 'k-', label='Sol norm')
plt.plot(PadektarsF, PadeNormF, 'b--',
label='Pade'' norm, {1} = {0}'.format(val, let))
plt.plot(RBktarsF, RBNormF, 'g--',
label='RB norm, {1} = {0}'.format(val, let))
plt.legend()
plt.grid()
p = plt.legend(loc = 'upper left')
plt.savefig('./normA' + str(i) + '.eps',
format='eps', dpi=1000)
plt.show()
plt.close()
plt.figure(figsize = (10, 5))
plt.semilogy(PadektarsF, PadeErrorF, 'b',
label='Pade'' error, {1} = {0}'.format(val, let))
plt.semilogy(RBktarsF, RBErrorF, 'g',
label='RB error, {1} = {0}'.format(val, let))
plt.legend()
plt.grid()
p = plt.legend(loc = 'lower right')
plt.savefig('./errorA' + str(i) + '.eps',
format='eps', dpi=1000)
plt.show()
plt.close()
plt.figure(figsize = (10, 5))
plt.semilogy(ktarsF, np.divide(PadeErrorF, solNormF), 'b',
label='Pade'' relative error, {1} = {0}'.format(val, let))
plt.semilogy(RBktarsF, np.divide(RBErrorF, solNormF), 'g',
label='RB relative error, {1} = {0}'.format(val, let))
plt.legend()
plt.grid()
p = plt.legend(loc = 'lower right')
plt.savefig('./errorAR' + str(i) + '.eps',
format='eps', dpi=1000)
plt.show()
plt.close()

Event Timeline

c4science · Help