Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/fenics/HelmholtzAirfoil.py
dce0cf46688amaster

Last Change
Raw File

HelmholtzAirfoil.py
View Options

from copy import copy
import fenics as fen
import numpy as np
import sympy as sp
from rrompy.hfengines.fenics import ScatteringProblemEngine as SPE
from rrompy.hsengines.fenics import HSEngine as HS
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB
from rrompy.reduction_methods.base import ParameterSweeper as Sweeper
from rrompy.sampling import QuadratureSampler as QS
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
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 * - 45 / 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/mesh/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 = 3)
DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))),
sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
if ttype == "simple":
solver = SPE()
solver.omega = k0
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = 'rest'
solver.DirichletDatum = [fen.Expression(x, degree = 3)\
for x in DirichletTerm]
plotter = HS(solver.V)
else:
print('NOPE')
# 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 = solver.DirichletDatum
baseRe = fen.project(baseRe, solver.V)
baseIm = fen.project(baseIm, solver.V)
uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
if ttype[: -1] == "augmented":
print('NOPE')
# uinc = [uinc, kappa * uinc]
if test == "solve":
aF = fen.interpolate(a, solver.V)
av = aF.vector()
uh = solver.solve(k0)
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, 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, mu0 = k0, plotDer = plotSamples,
approxParameters = parPade)
approxRB = TRB(solver, plotter, mu0 = k0, approxParameters = parRB)
else:
params = {'N':8, 'M':7, 'R':9, 'S':9, 'POD':True,
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = LP(solver, plotter, mu0 = np.mean([kLeft, kRight]),
plotSnap = plotSamples, approxParameters = parPade)
approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
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())
print('\nPoles RB:')
print(approxRB.getPoles())
uHF = approxPade.getHF(ktar)
plotter.plot(uHF, name = 'u_ex')
approxPade.plotApp(ktar, name = 'u_Pade''')
approxRB.plotApp(ktar, name = 'u_RB')
approxPade.plotErr(ktar, name = 'errPade''')
approxRB.plotErr(ktar, name = 'errRB')
elif test in ["TaylorSweep", "LagrangeSweep"]:
if test == "TaylorSweep":
shift = 10
nsets = 3
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 + 1,'R':stride*i+shift + 1}
approxPade = TP(solver, plotter, mu0 = kappa,approxParameters = params)
approxRB = TRB(solver, plotter, mu0 = kappa, approxParameters = params)
else:
shift = 10
nsets = 3
stride = 3
Smax = stride * (nsets - 1) + shift + 2
paramsPade = {'S':Smax, 'POD':True,
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
paramsRB = copy(paramsPade)
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, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsPade)
approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsRB)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('../data/output/HelmholtzAirfoil'+test[:-5]+'PadeFE.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep('../data/output/HelmholtzAirfoil'+test[:-5]+'RBFE.dat')
for i in range(nsets):
if test == "TaylorSweep":
val = stride*i+shift+1
PadeOutput = sweeper.read(filenamePade, {'N':val},
['muRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'E':val},
['muRe', 'AppNorm', 'ErrNorm'])
let = 'E'
else:
val = stride*i+shift+2
PadeOutput = sweeper.read(filenamePade, {'S':val},
['muRe', 'HFNorm', 'AppNorm', 'ErrNorm'])
RBOutput = sweeper.read(filenameRB, {'S':val},
['muRe', 'AppNorm', 'ErrNorm'])
let = 'S'
ktarsF = PadeOutput['muRe']
solNormF = PadeOutput['HFNorm']
PadektarsF = PadeOutput['muRe']
PadeNormF = PadeOutput['AppNorm']
PadeErrorF = PadeOutput['ErrNorm']
RBktarsF = RBOutput['muRe']
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(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()
c4science · Help