Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/from_papers/pod_scatteringAirfoil.py
76c35e029615feature_greedy_coarseni…

Last Change
Raw File

pod_scatteringAirfoil.py
View Options

from copy import deepcopy as copy
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HSP
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.utilities.parameter_sweeper import ParameterSweeper as Sweeper
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
from rrompy.utilities.base.fenics import fenONE
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
####################
homog = True
#homog = False
####################
test = "solve"
test = "Taylor"
test = "Lagrange"
test = "TaylorSweep"
test = "LagrangeSweep"
plotSamples = True
k0 = 10
kLeft, kRight = 8 + 0.j, 12 + 0.j
ktar = 11
ktars = np.linspace(8, 12, 21) + 0.j
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
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
###
solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
###
if test == "solve":
uinc = solver.liftDirichletData(k0)
if homog:
uhtot = solver.solve(k0, homogeneized = homog)
uh = uhtot + uinc
else:
uh = solver.solve(k0, homogeneized = homog)
uhtot = uh - uinc
print(solver.norm(uh))
print(solver.norm(uhtot))
solver.plot(fen.project(a, solver.V).vector(), what = 'Real',
name = 'a')
solver.plot(uinc, what = 'Real', name = 'u_inc')
solver.plot(uh, what = 'ABS')
solver.plot(uhtot, what = 'ABS', name = 'u + u_inc')
elif test in ["Taylor", "Lagrange"]:
if test == "Taylor":
params = {'N':8, 'M':8, 'R':8, 'E':8,
'sampleType':'Arnoldi', 'POD':True}
parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
approxPade = TP(solver, mu0 = k0, approxParameters = parPade,
verbosity = verb, homogeneized = homog)
approxRB = TRB(solver, mu0 = k0, approxParameters = parRB,
verbosity = verb, homogeneized = homog)
else:
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler', 'basis'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parPade, verbosity = verb,
homogeneized = homog)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB, verbosity = verb,
homogeneized = homog)
approxPade.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxPade.plotSamples()
approxPade.plotHF(ktar, name = 'u_HF')
approxPade.plotApprox(ktar, name = 'u_Pade''')
approxPade.plotErr(ktar, name = 'err_Pade''')
approxPade.plotRes(ktar, name = 'res_Pade''')
approxRB.plotApprox(ktar, name = 'u_RB')
approxRB.plotErr(ktar, name = 'err_RB')
approxRB.plotRes(ktar, name = 'res_RB')
HFNorm, RHSNorm = approxPade.normHF(ktar), approxPade.normRHS(ktar)
PadeRes, PadeErr = approxPade.normRes(ktar), approxPade.normErr(ktar)
RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar)
print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm))
print('PadeRes:\t{}\nPadeErr:\t{}'.format(PadeRes, PadeErr))
print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr))
print('\nPoles Pade'':')
print(approxPade.getPoles())
elif test in ["TaylorSweep", "LagrangeSweep"]:
if test == "TaylorSweep":
shift = 10
nsets = 2
stride = 3
Emax = stride * (nsets - 1) + shift + 1
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift + 1, 'M':stride*i+shift+1,
'E':stride*i+shift + 1}
paramsSetsRB[i] = {'E':stride*i+shift + 1,'R':stride*i+shift + 2}
approxPade = TP(solver, mu0 = kappa,approxParameters = params,
verbosity = verb, homogeneized = homog)
approxRB = TRB(solver, mu0 = kappa, approxParameters = params,
verbosity = verb, homogeneized = homog)
else:
shift = 10
nsets = 2
stride = 3
Smax = stride * (nsets - 1) + shift + 2
paramsPade = {'S':Smax, 'POD':True, 'basis':"CHEBYSHEV",
'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+1,
'S': stride*i+shift+2}
paramsSetsRB[i] = {'R': stride*i+shift+2, 'S': stride*i+shift+2}
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsPade, verbosity = verb,
homogeneized = homog)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsRB, verbosity = verb,
homogeneized = homog)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
homogMSG = ""
if not homog: homogMSG = "Non"
filenamebase = '../data/output/airfoil' + test[:-5] + homogMSG + "Homog"
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase +'Pade.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase +'RB.dat')
if test == "TaylorSweep":
constr = ['E']
else:
constr = ['S']
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'],
['normHF', 'normApp'], constr, onePlot = True,
save = filenamebase + 'Norm',
saveFormat = "png", labels = ["Pade'", "RB"])
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normResRel'],
constr, save = filenamebase + 'Res',
saveFormat = "png", labels = ["Pade'", "RB"])
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normErrRel'],
constr, save = filenamebase + 'Err',
saveFormat = "png", labels = ["Pade'", "RB"])
c4science · Help