Page MenuHomec4science
Files F63027529

HelmholtzAirfoilTaylor.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 17, 06:41

HelmholtzAirfoilTaylor.py
View Options

#!/usr/bin/env python3
import numpy as np
from context import FreeFemHelmholtzScatteringEngine as HFS
from context import FreeFemHelmholtzScatteringAugmentedEngine as HFSA
from context import FreeFemHSEngine as HS
from context import FreeFemHSAugmentedEngine 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
kappa = 10
theta = PI * -30 / 180.
mu = 1.1
epsilon = .1
V = ("real epsilon = {}, mu = {};\n"
"macro sign(f) (f >= 0 ? 1. : -1.)//\n"
"macro mBase1() (x^2 - x + y^2)//\n"
"macro mBase2() (((x^2 + y^2) / ((x - 1)^2 + y^2))^.25)//\n"
"macro mBase3() (- y / mBase1)//\n"
"macro mPhi1() (atan(mBase3) / 2)//\n"
"macro mPhi2() (mPhi1 - pi/2 * sign(mBase3))//\n"
"macro mK1() ((((1 + cos(mPhi1) / mBase2 + (.5 / mBase2)^2.) / (1 - 2. * cos(mPhi1) / mBase2 + mBase2^-2.))^.5 - mu)/ epsilon)//\n"
"macro mK2() ((((1 + cos(mPhi2) / mBase2 + (.5 / mBase2)^2.) / (1 - 2. * cos(mPhi2) / mBase2 + mBase2^-2.))^.5 - mu)/ epsilon)//\n"
"macro mHep1() (.9 * .5 * (1. + mK1 + sin(mK1) / pi) + .1)//\n"
"macro mHep2() (.9 * .5 * (1. + mK2 + sin(mK2) / pi) + .1)//\n"
"int n = 100, R = 2;\n"
"border wing(t = 0, 2 * pi){{x = cos(t)/3 + 5./12 + (3*cos(t) - 3/4) / (17 - 8*cos(t)); y = sin(t)/3 - 3*sin(t) / (17 - 8*cos(t)); label = 1;}}\n"
"border Inf(t = 0, 2 * pi){{x = R * cos(t); y = R * sin(t); label = 2;}}\n"
"mesh Th = buildmesh(wing(-n) + Inf(n));\n"
"load \"Element_P3\";\n"
"fespace V(Th, P3);").format(epsilon, mu)
u0 = "- exp(1.i * {} * ({} * x + {} * y))".format(kappa, np.cos(theta),
np.sin(theta))
a = ("(mBase1 >= 0 ? (mK1>= -1. ? (mK1 <= 1. ? mHep1 : 1) : .1) : "
"(mK2 >= -1. ? (mK2 <= 1. ? mHep2 : 1) : .1))")
if ttype == "simple":
solver = HFS(V, k0, diffusivity = a, DirichletBoundary = "1",
RobinBoundary = "2", DirichletDatum = u0)
plotter = HS(solver.V)
else:
if ttype[-1] == "I": constraintType = "IDENTITY"
else: constraintType = "MASS"
solver = HFSA(V, k0, diffusivity = a, DirichletBoundary = "1",
RobinBoundary = "2", DirichletDatum = u0,
constraintType = constraintType)
plotter = HSA(solver.V, d = 2)
uinc = solver.liftDirichletData()
if ttype[: -1] == "augmented":
uinc = [uinc[: int(len(uinc)/2)], kappa * uinc[: int(len(uinc)/2)]]
if test == "solve":
av = plotter.convExpression(a)
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]+'PadeFF.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep('../Data/HelmholtzAirfoil'+test[:-5]+'RBFF.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.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.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)

Event Timeline

c4science · Help