Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FreeFem++/HelmholtzAirfoilTaylor.py
6e6c171cf128master

Last Change
Raw File

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)
c4science · Help