Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/greedy/squareScatteringAirfoil.py
c191c12b6f43master

Last Change
Raw File

squareScatteringAirfoil.py
View Options

import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade
from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangeRBGreedy as RB
from rrompy.utilities.base.fenics import fenONE
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
homog = True
homog = False
k0s = np.linspace(5, 20, 50) - .5j
k0 = np.mean(k0s)
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTrainingPoints':500, 'Delta':0,
'greedyTol':1e-2, 'nTestPoints':2}
#########
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.real(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
#########
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb, homogeneized = homog)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb,
homogeneized = homog)
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApp(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estNormer.norm(approx.getRes(k0s[j]))
/ approx.estNormer.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
c4science · Help