Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/airfoil/airfoil_engine.py
f972e1cdd3b2feature_greedy_coarseni…

Last Change
Raw File

airfoil_engine.py
View Options

import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.utilities.fenics import fenONE
PI = np.pi
class AirfoilScatteringEngine(HSP):
def __init__(self, kappa:float, theta:float,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(R = 5, kappa = kappa, theta = theta, n = 1,
degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
mesh = fen.Mesh('../data/mesh/airfoil2412_1.xml')
self.V = fen.FunctionSpace(mesh, "P", 1)
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))
self.DirichletDatum = [u0R, u0I]
mu = 1.1
epsilon = .1
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)
self.diffusivity = a
c4science · Help