Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/helmholtz_resonator.py
5bd9689f4473feature_greedy_coarseni…

Last Change
Raw File

helmholtz_resonator.py
View Options

from matplotlib import pyplot as plt
import fenics as fen
import mshr
import ufl
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
p = plt.jet()
n = 50
boundary = mshr.Polygon([fen.Point(0, 0), fen.Point(6, 0), fen.Point(6, 1),
fen.Point(1.3, 1), fen.Point(1.3, 1.2),
fen.Point(1.65, 1.2), fen.Point(1.65, 2.2),
fen.Point(.65, 2.2), fen.Point(.65, 1.2),
fen.Point(1, 1.2), fen.Point(1, 1), fen.Point(0, 1)])
mesh = mshr.generate_mesh(boundary, n)
class Resonator(fen.SubDomain):
def inside(self, x, on_boundary):
return fen.between(x[1], (1.25, 2.25))
class Noslip(fen.SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class Inlet(fen.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fen.near(x[0], 0.)
class Outlet(fen.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fen.near(x[0], 6.)
class Liner(fen.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fen.near(x[1], 2.25)
resonator = Resonator()
noslip = Noslip()
inlet = Inlet()
outlet = Outlet()
liner = Liner()
sub_domains = fen.MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
resonator.mark(sub_domains, 1)
boundaries = fen.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
noslip.mark(boundaries, 0)
inlet.mark(boundaries, 1)
outlet.mark(boundaries, 2)
liner.mark(boundaries, 3)
for k in range(10):
kappa = .25 + .05 * k
ZR = 10.
x, y = fen.SpatialCoordinate(mesh)[:]
solver = HPE()
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.omega = kappa
solver.RobinBoundary = lambda x, on_b: on_b and (fen.near(x[0], 6.)
or fen.near(x[1], 2.25))
solver.NeumannBoundary = "REST"
solver.signR = + 1.
solver.NeumannDatum = [fen.Constant(0.),
ufl.conditional(ufl.And(ufl.gt(y, 0.), ufl.lt(y, 1.)),
fen.Constant(kappa), fen.Constant(0.))]
solver.RobinDatumH = [fen.Constant(0.),
ufl.conditional(ufl.gt(y, 1.25),
fen.Constant(kappa / ZR),
fen.Constant(kappa))]
uh = solver.solve(kappa)
solver.plot(uh, name = "k={}".format(kappa))
print(solver.norm(uh))
c4science · Help