Adapted from PyMFEM/ex9.py. Compare with the original Example 9 in MFEM.
This example code solves the time-dependent advection equation
∂u∂t+v⋅∇u=0,where v is a given fluid velocity, and
u0(x)=u(0,x)is a given initial condition.
The example demonstrates the use of Discontinuous Galerkin (DG) bilinear forms in MFEM (face integrators), the use of explicit ODE time integrators, the definition of periodic boundary conditions through periodic meshes, as well as the use of GLVis for persistent visualization of a time-evolving solution.
from __future__ import print_function
from os.path import expanduser, join
import time
import numpy as np
from numpy import sqrt, pi, cos, sin, hypot, arctan2
from scipy.special import erfc
import mfem.ser as mfem
from mfem.ser import intArray
from glvis import glvis
# 1. Setup problem parameters
problem = 0
order = 3
vis_steps = 5
ode_solver_type = 4
sleep_interval = .05
if problem == 0:
# Translating hump
mesh = "periodic-hexagon.mesh"
ref_levels = 2
dt = 0.01
t_final = 10
elif problem == 1:
# Rotating power sines
mesh = "periodic-square.mesh"
ref_levels = 2
dt = 0.005
t_final = 9
elif problem == 2:
# Rotating trigs
mesh = "disc-nurbs.mesh"
ref_levels = 3
dt = 0.005
t_final = 9
elif problem == 3:
# Twisting sines
mesh = "periodic-square.mesh"
ref_levels = 4
dt = 0.0025
t_final = 9
vis_steps = 20
# 2. Download the mesh from GitHub and read it.
# We can handle geometrically periodic meshes in this code.
mesh_url = f"https://raw.githubusercontent.com/mfem/mfem/master/data/{mesh}"
mesh_str = !curl -s {mesh_url}
with open(mesh, "w") as f:
f.write("\n".join(mesh_str))
mesh = mfem.Mesh(mesh, 1,1)
dim = mesh.Dimension()
# 3. Define the ODE solver used for time integration.
# Several explicit Runge-Kutta methods are available.
ode_solver = None
if ode_solver_type == 1: ode_solver = mfem.ForwardEulerSolver()
elif ode_solver_type == 2: ode_solver = mfem.RK2Solver(1.0)
elif ode_solver_type == 3: ode_solver = mfem.RK3SSolver()
elif ode_solver_type == 4: ode_solver = mfem.RK4Solver()
elif ode_solver_type == 6: ode_solver = mfem.RK6Solver()
else: print(f"Unknown ODE solver type: {ode_solver_type}")
# 4. Refine the mesh to increase the resolution. In this example we do
# 'ref_levels' of uniform refinement. If the mesh is of NURBS type,
# we convert it to a (piecewise-polynomial) high-order mesh.
for lev in range(ref_levels):
mesh.UniformRefinement();
if mesh.NURBSext:
mesh.SetCurvature(max(order, 1))
bb_min, bb_max = mesh.GetBoundingBox(max(order, 1));
# 5. Define the discontinuous DG finite element space of the given
# polynomial order on the refined mesh.
fec = mfem.DG_FECollection(order, dim)
fes = mfem.FiniteElementSpace(mesh, fec)
print("Number of unknowns: " + str(fes.GetVSize()))
# Define coefficient using VecotrPyCoefficient and PyCoefficient
# A user needs to define EvalValue method
class velocity_coeff(mfem.VectorPyCoefficient):
def EvalValue(self, x):
dim = len(x)
center = (bb_min + bb_max)/2.0
# map to the reference [-1,1] domain
X = 2 * (x - center) / (bb_max - bb_min)
if problem == 0:
if dim == 1: v = [1.0,]
elif dim == 2: v = [sqrt(2./3.), sqrt(1./3)]
elif dim == 3: v = [sqrt(3./6.), sqrt(2./6), sqrt(1./6.)]
elif (problem == 1 or problem == 2):
# Clockwise rotation in 2D around the origin
w = pi/2
if dim == 1: v = [1.0,]
elif dim == 2: v = [w*X[1], - w*X[0]]
elif dim == 3: v = [w*X[1], - w*X[0], 0]
elif (problem == 3):
# Clockwise twisting rotation in 2D around the origin
w = pi/2
d = max((X[0]+1.)*(1.-X[0]),0.) * max((X[1]+1.)*(1.-X[1]),0.)
d = d ** 2
if dim == 1: v = [1.0,]
elif dim == 2: v = [d*w*X[1], - d*w*X[0]]
elif dim == 3: v = [d*w*X[1], - d*w*X[0], 0]
return v
class u0_coeff(mfem.PyCoefficient):
def EvalValue(self, x):
dim = len(x)
center = (bb_min + bb_max)/2.0
# map to the reference [-1,1] domain
X = 2 * (x - center) / (bb_max - bb_min)
if (problem == 0 or problem == 1):
if dim == 1:
return exp(-40. * (X[0]-0.5)**2)
elif (dim == 2 or dim == 3):
rx = 0.45; ry = 0.25; cx = 0.; cy = -0.2; w = 10.
if dim == 3:
s = (1. + 0.25*cos(2 * pi * x[2]))
rx = rx * s
ry = ry * s
return ( erfc( w * (X[0]-cx-rx)) * erfc(-w*(X[0]-cx+rx)) *
erfc( w * (X[1]-cy-ry)) * erfc(-w*(X[1]-cy+ry)) )/16
elif problem == 2:
rho = hypot(x[0], x[1])
phi = arctan2(x[1], x[0])
return (sin(pi * rho) **2) * sin(3*phi)
elif problem == 3:
return sin(pi * X[0]) * sin(pi * X[1])
return 0.0
# Inflow boundary condition (zero for the problems considered in
# this example)
class inflow_coeff(mfem.PyCoefficient):
def EvalValue(self, x):
return 0
# 6. Set up and assemble the bilinear and linear forms corresponding to
# the DG discretization. The DGTraceIntegrator involves integrals over
# mesh interior faces.
velocity = velocity_coeff(dim)
inflow = inflow_coeff()
u0 = u0_coeff()
m = mfem.BilinearForm(fes)
m.AddDomainIntegrator(mfem.MassIntegrator())
k = mfem.BilinearForm(fes)
k.AddDomainIntegrator(mfem.ConvectionIntegrator(velocity, -1.0))
k.AddInteriorFaceIntegrator(mfem.TransposeIntegrator(mfem.DGTraceIntegrator(velocity, 1.0, -0.5)))
k.AddBdrFaceIntegrator(mfem.TransposeIntegrator(mfem.DGTraceIntegrator(velocity, 1.0, -0.5)))
b = mfem.LinearForm(fes)
b.AddBdrFaceIntegrator(mfem.BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5))
m.Assemble()
m.Finalize()
skip_zeros = 0
k.Assemble(skip_zeros)
k.Finalize(skip_zeros)
b.Assemble()
# Initial conditions
u = mfem.GridFunction(fes)
u.ProjectCoefficient(u0)
g = glvis((mesh, u))
g.render()
class FE_Evolution(mfem.PyTimeDependentOperator):
def __init__(self, M, K, b):
mfem.PyTimeDependentOperator.__init__(self, M.Size())
self.K = K
self.M = M
self.b = b
self.z = mfem.Vector(M.Size())
self.zp = np.zeros(M.Size())
self.M_prec = mfem.DSmoother()
self.M_solver = mfem.CGSolver()
self.M_solver.SetPreconditioner(self.M_prec)
self.M_solver.SetOperator(M)
self.M_solver.iterative_mode = False
self.M_solver.SetRelTol(1e-9)
self.M_solver.SetAbsTol(0.0)
self.M_solver.SetMaxIter(100)
self.M_solver.SetPrintLevel(0)
def Mult(self, x, y):
self.K.Mult(x, self.z);
self.z += b;
self.M_solver.Mult(self.z, y)
adv = FE_Evolution(m.SpMat(), k.SpMat(), b)
ode_solver.Init(adv)
t = 0.0; ti = 0;
time.sleep(1)
while True:
if t > t_final - dt/2: break
t, dt = ode_solver.Step(u, t, dt);
ti = ti + 1
if ti % vis_steps == 0:
g.update((mesh, u))
time.sleep(sleep_interval)
print(f"time step: {ti}, time: {t:.2e}", end="\r")