I’ve been trying to cook up a model that solves the (steady) compressible Euler equations with a Discontinuous Galerkin (DG) method, but for whatever reason I’m not getting any useful results. The approach that I’ve taken thus far is
\begin{equation} -\int_\Omega \mathcal{F} \colon \nabla \bf{v}\ d\Omega + \int_{\Gamma} \mathcal{F}^* \colon (\bf{v} \otimes \bf{n})\ d\Gamma + \int_{\partial \Omega} \mathcal{F}_{\partial \Omega}^* \colon [[ \bf{v} ]]\ d\partial\Omega = 0, \end{equation}
where \Gamma are the internal cell boundaries and \partial \Omega are the external cell boundaries. For flux \mathcal{F} I’m using the formulation that is given in the dolfin_dg paper (in listing 14). For the internal cell boundary flux \mathcal{F}^* I’m using an upwind scheme. The external boundary flux \mathcal{F}_{\partial \Omega}^* is used to weakly enforce the Dirichlet boundary conditions where necessary.
The geometry and mesh that are used are taken from this DFG cylinder flow case. As mentioned in eq. 3 of this paper, for a subsonic Euler flow we should prescribe the density, pressure and velocity at the inlet, and the pressure at the outlet. To enforce the slip wall boundary conditions on the cylinder and other walls I use the approach that is mentioned on page 18 of the aforementioned dolfin_dg paper. I don’t know why, but the NewtonSolver does not give any useful output (the norm of the solution is infinite).
Here is the code that I’m trying to get to work, including mesh construction:
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
from copy import deepcopy
import numpy as np
from dolfinx import default_real_type, fem, io
from dolfinx.fem.petsc import NonlinearProblem, assemble_vector
from dolfinx.nls.petsc import NewtonSolver
from basix.ufl import element, mixed_element
from ufl import (
CellDiameter,
FacetNormal,
TestFunction,
TrialFunction,
conditional,
derivative,
dS,
ds,
dx,
dot,
grad,
gt,
inner,
outer,
Measure,
split,
as_vector,
as_matrix,
)
def F_c(u_vec):
rho, u1, u2, E = u_vec[0], u_vec[1]/u_vec[0], u_vec[2]/u_vec[0], u_vec[3]/u_vec[0]
p_comp = (gamma - 1.0)*rho*(E - 0.5*(u1**2 + u2**2))
H = E + p_comp/rho
return as_matrix([[rho*u1, rho*u2 ],
[rho*u1**2 + p_comp, rho*u1*u2 ],
[rho*u1*u2, rho*u2**2 + p_comp],
[rho*H*u1, rho*H*u2 ]])
def F_prescribed_pressure(u_vec, p_inp):
rho, u1, u2, E = u_vec[0], u_vec[1]/u_vec[0], u_vec[2]/u_vec[0], u_vec[3]/u_vec[0]
H = E + p_inp/rho
return as_matrix([[rho*u1, rho*u2 ],
[rho*u1**2 + p_inp, rho*u1*u2 ],
[rho*u1*u2, rho*u2**2 + p_inp],
[rho*H*u1, rho*H*u2 ]])
def jump(phi, n):
return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))
# ----------------------------------------------- #
# Generate mesh around cylinder
# ----------------------------------------------- #
gmsh.initialize()
L = 2.2
H_domain = 0.8
r = 0.05
c_x = 0.3
c_y = H_domain / 2.
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H_domain, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
fluid_marker = 1
if mesh_comm.rank == model_rank:
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H_domain / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H_domain / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H_domain, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
# Create distance field from obstacle.
# Add threshold of mesh sizes based on the distance field
# LcMax - /--------
# /
# LcMin -o---------/
# | | |
# Point DistMin DistMax
res_min = r / 25
if mesh_comm.rank == model_rank:
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.1 * H_domain)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", H_domain)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 8)
# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
# gmsh.option.setNumber("Mesh.RecombineAll", 1)
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(1)
gmsh.model.mesh.optimize("Netgen")
msh, _, fts = io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
fts.name = "Facet markers"
h_msh = CellDiameter(msh)
n_msh = FacetNormal(msh)
print("mesh generated")
# Create connectivties required for defining integration entities
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_entities(fdim)
msh.topology.create_connectivity(fdim, tdim)
msh.topology.create_connectivity(tdim, fdim)
ds = Measure('ds', domain=msh, subdomain_data=fts)
pol_deg = 0 # polynomial degree of basis functions
rho_0 = 1.
M_0 = 0.5
p_0 = 1.
gamma = 1.4
# Inlet conditions
c_0 = abs(gamma*p_0/rho_0)**0.5
n_in = dolfinx.fem.Constant(
msh, np.array([1., 0.], dtype=np.double))
u_ref = dolfinx.fem.Constant(msh, M_0*c_0)
# u_ref = dolfinx.fem.Constant(msh, 0.5)
rho_in = dolfinx.fem.Constant(msh, rho_0)
u_in = u_ref * n_in
rhou_in = rho_in * u_in
u_in_np = u_ref.value * n_in.value
E_in = dolfinx.fem.Constant(msh, p_0 / (gamma - 1.) + (rho_0 / 2.) * np.dot(u_in_np, u_in_np))
p_out = dolfinx.fem.Constant(msh, p_0)
# define the function spaces for the various quantities
DG_scalar = element("Lagrange", msh.basix_cell(), pol_deg, discontinuous=True, dtype=default_real_type)
DG_vel = element("Lagrange", msh.basix_cell(), pol_deg, shape=(msh.geometry.dim,), discontinuous=True, dtype=default_real_type)
DG_var_space = fem.functionspace(msh, mixed_element([DG_scalar, DG_vel, DG_scalar]))
# DG_vel_space = fem.functionspace(msh, DG_vel)
u = fem.Function(DG_var_space)
v = TestFunction(DG_var_space)
rho_in_expr = dolfinx.fem.Expression(rho_in, DG_var_space.sub(0).element.interpolation_points())
rhou_in_expr = dolfinx.fem.Expression(rhou_in, DG_var_space.sub(1).sub(0).element.interpolation_points())
E_in_expr = dolfinx.fem.Expression(E_in, DG_var_space.sub(2).element.interpolation_points())
# Interpolate initial conditions for each variable
u.sub(0).interpolate(rho_in_expr)
u.sub(1).interpolate(rhou_in_expr)
u.sub(2).interpolate(E_in_expr)
u.x.scatter_forward()
rho, rhou, E = split(u)
print("Constructing weak form")
# First we define the various internal and boundary flux functions
FG_interior = F_c(u)
u_inlet = as_vector([rho_in, rhou_in[0], rhou_in[1], E])
FG_inlet = F_prescribed_pressure(u_inlet, p_out)
FG_out = F_prescribed_pressure(u, p_out)
slip_proj = as_matrix(((1, 0, 0, 0),
(0, 1-2*n_msh[0]**2, -2*n_msh[0]*n_msh[1], 0),
(0, -2*n_msh[0]*n_msh[1], 1-2*n_msh[1]**2, 0),
(0, 0, 0, 1)))
u_slipwall = slip_proj * u
FG_wall = F_c(u_slipwall)
# we construct the upwind flux on the internal cell faces
u_vel = as_vector([u[1]/u[0], u[2]/u[0]])
lmbda = conditional(gt(dot(u_vel, n_msh), 0), 1, 0)
FG_upwind = lmbda("+") * FG_interior("+") + lmbda("-") * FG_interior("-")
# we construct the weak form, starting with the volume terms
weakform = - inner(FG_interior, grad(v)) * dx
# then the internal cell faces
weakform += inner(FG_upwind, jump(v, n_msh)) * dS
# and lastly we construct the fluxes on the external cell faces
weakform += inner(FG_out, outer(v, n_msh)) * ds(outlet_marker)
weakform += inner(FG_inlet, outer(v, n_msh)) * ds(inlet_marker)
weakform += inner(FG_wall, outer(v, n_msh)) * ( ds(wall_marker) + ds(obstacle_marker) )
du = TrialFunction(DG_var_space)
J = derivative(weakform, u, du)
weak_form, J_form = fem.form(weakform), fem.form(J)
problem = NonlinearProblem(weak_form, u, J=J_form) #, bcs=bcs_dirichlet)
# # problem = NonlinearPDE_SNESProblem(weak_residual, u, [], M_mass)
solver = NewtonSolver(msh.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-16
solver.atol = 1e-12
solver.report = True
solver.error_on_nonconvergence = False
solver.max_it = 3
def updater(solver, dx, x):
tau = 1.
dx_norm = dx.norm()
if dx_norm >= 1e-10:
theta = min((2.0*tau/dx_norm)**(1./2.), 1.0)
else:
theta = 1.
print("dx norm: {}".format(dx_norm))
print("relaxation factor: {}".format(theta))
x.axpy(-theta, dx)
solver.set_update(updater)
ksp = solver.krylov_solver
opts = PETSc.Options()
ksp.setOptionsPrefix("")
option_prefix = ksp.getOptionsPrefix()
# for a direct solver
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()
print("Starting Newton solve...")
solver.solve(u)
u.x.scatter_forward()
print("Finished Newton solve...")
res_vec_new = assemble_vector(weak_form)
res_vec_norm_new = np.linalg.norm(res_vec_new.array)
print("Max abs residual vector entry: {}".format(np.abs(res_vec_new).max()))
print("Residual norm: {}".format(res_vec_norm_new))
rho_vis = fem.Function(DG_var_space.sub(0).collapse()[0])
rhou_vis = fem.Function(DG_var_space.sub(1).collapse()[0])
E_vis = fem.Function(DG_var_space.sub(2).collapse()[0])
rho_file = io.VTXWriter(msh.comm, "rho.bp", rho_vis)
rhou_file = io.VTXWriter(msh.comm, "rhou.bp", rhou_vis)
E_file = io.VTXWriter(msh.comm, "E.bp", E_vis)
rho_vis.interpolate(u.sub(0).collapse())
rhou_vis.interpolate(u.sub(1).collapse())
E_vis.interpolate(u.sub(2).collapse())
rho_file.write(0)
rhou_file.write(0)
E_file.write(0)
rho_file.close()
rhou_file.close()
E_file.close()
I’m hoping that someone might have an inkling of what I should do to get this example to run. Let me know if there’s any information missing from the description above. I’ll tag @nate here to see what he thinks, given that he likely wrote the original DG implementation of this model.