Hello everyone, I’m working on the DG solver of the CHNS (Cahn-Hilliard Navier-Stokes) for the multiphase flow. The SIPG is used to treat the diffusion term and an upwind flux scheme is used for the convection term. Now my solver can bear the large ratio of the density and viscosity while it can also bear the small absolute viscosity. As you can see, when I use the parameters from the hysingtest2, it can run, but when I use the real dam-break test parameters, the Gibbs oscillations emerge and finally the solution diverges or blows up. For the rising bubble test, I didn’t observe the oscillation of the velocity field but the velocity singularity emerges, which causes the time step to become so small and leads to failure of the computation.
Any suggestion is appreciated. Thanks a lot！
#!/usr/bin/env python3
"""
Cahn-Hilliard Navier-Stokes rising bubble problem - DG implementation
- BDF1 time discretization (both NS and CH)
- DG SIPG + lagged upwind advection
- Physical coefficients and Korteweg force evaluated at current fields
- Upwind direction lagged from previous step to keep Newton Jacobian differentiable
"""
import numpy as np
import ufl
import basix.ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh, fem, log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter
import dolfinx
from pathlib import Path
print("dolfinx version :", dolfinx.__version__)
print("basix version :", basix.__version__)
# ====================== Geometry and mesh ======================
Lx, Ly = 1.5, 0.75
Nx, Ny = 150, 75
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([Lx, Ly])],
[Nx, Ny],
cell_type=mesh.CellType.triangle
)
boundary_facets = mesh.locate_entities_boundary(
domain, domain.topology.dim - 1, lambda x: np.full(x.shape[1], True, dtype=bool)
)
facet_tags = mesh.meshtags(
domain, domain.topology.dim - 1, boundary_facets, np.full(len(boundary_facets), 1, dtype=np.int32)
)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dS = ufl.dS(domain=domain)
# ====================== DG function spaces ======================
cell_name = domain.topology.cell_name()
k = 2
ve = basix.ufl.element("Discontinuous Lagrange", cell_name, k, shape=(domain.geometry.dim,))
pe = basix.ufl.element("Discontinuous Lagrange", cell_name, k - 1)
ce = basix.ufl.element("Discontinuous Lagrange", cell_name, k - 1)
me = basix.ufl.mixed_element([ve, pe, ce, ce])
W = fem.functionspace(domain, me)
w = fem.Function(W)
w_n = fem.Function(W)
u, p, c, mu_c = ufl.split(w)
v, q, phi, psi = ufl.TestFunctions(W)
# ====================== Lagged fields ======================
# For BDF1 extrapolation and nonlinear coefficient evaluation
V_u, _ = W.sub(0).collapse()
V_c, _ = W.sub(2).collapse()
u_n = fem.Function(V_u, name="u_n")
c_n = fem.Function(V_c, name="c_n")
# ====================== Physical constants ======================
mesh_size = min(Lx/Nx, Ly/Ny)
dt_val = mesh_size * 0.05
print(dt_val)
dt = fem.Constant(domain, default_scalar_type(dt_val))
hx = Lx / Nx
hy = Ly / Ny
h_max = np.sqrt(hx**2 + hy**2)
alpha = 0.75
beta = 0.05
eps_val = alpha * h_max
M_val = beta * (eps_val ** 2)
eps = fem.Constant(domain, default_scalar_type(eps_val))
M = fem.Constant(domain, default_scalar_type(M_val))
sigma = fem.Constant(domain, default_scalar_type(0.0752))
g_vec = fem.Constant(domain, default_scalar_type((0.0, -9.8)))
lambda_val = 3.0 * np.sqrt(2.0) * sigma.value * eps.value
lmbda = fem.Constant(domain, default_scalar_type(lambda_val))
rho_liquid, rho_gas = 1000.0, 1.0
mu_liquid, mu_gas = 0.5, 0.005
def rho(c_var):
c_b = ufl.max_value(0.0, ufl.min_value(1.0, c_var))
return c_b * rho_liquid + (1.0 - c_b) * rho_gas
def mu_kin(c_var):
c_b = ufl.max_value(0.0, ufl.min_value(1.0, c_var))
return c_b * mu_liquid + (1.0 - c_b) * mu_gas
# ====================== Initial condition ======================
def initial_water_column(x):
phi_x = 0.5 * (1.0 + np.tanh((0.5 - x[0]) / (np.sqrt(2.0) * eps_val)))
phi_y = 0.5 * (1.0 + np.tanh((0.5 - x[1]) / (np.sqrt(2.0) * eps_val)))
return phi_x * phi_y
w_n.sub(2).interpolate(initial_water_column)
w.x.array[:] = w_n.x.array[:]
u_n.interpolate(w_n.sub(0))
c_n.interpolate(w_n.sub(2))
# ====================== DG helper quantities ======================
h = ufl.CellDiameter(domain)
n = ufl.FacetNormal(domain)
alpha_sipg_ns = fem.Constant(domain, default_scalar_type(100.0 * k**2))
alpha_sipg_ch = fem.Constant(domain, default_scalar_type(100.0 * k**2))
alpha_bd = alpha_sipg_ns * 50
alpha_p = fem.Constant(domain, default_scalar_type(2.0))
beta_ch = 20
def jump(phi, n):
return ufl.outer(phi("+"), n("+")) + ufl.outer(phi("-"), n("-"))
# ====================== Residual forms (BDF1 + lagged upwind) ======================
dfdc = 2.0 * c * (1.0 - c) * (1.0 - 2.0 * c)
def make_F_NS():
term_time = ufl.inner( (rho(c)*u - rho(c_n)*u_n) / dt , v) * ufl.dx
# Lag only the upwind *direction* from the previous step (u_n).
# The advected value uses the current unknown u; this keeps the Newton Jacobian differentiable.
lmbda_upw = ufl.conditional(ufl.gt(ufl.dot(u_n, n), 0), 1, 0)
u_uw = lmbda_upw("+") * u("+") + lmbda_upw("-") * u("-")
term_conv_vol = -rho(c) * ufl.inner(ufl.outer(u,u), ufl.grad(v)) * ufl.dx
term_conv_surf = (
rho(c("+")) * ufl.inner(ufl.dot(u, n)("+") * u_uw, v("+")) * dS +
rho(c("-")) * ufl.inner(ufl.dot(u, n)("-") * u_uw, v("-")) * dS
)
term_conv_bd = 0
F_conv = term_conv_vol + term_conv_surf + term_conv_bd
mu = mu_kin(c)
def epsilon(u_var):
return ufl.sym(ufl.grad(u_var))
F_visc = (
ufl.inner(2.0 * mu * epsilon(u), epsilon(v)) * ufl.dx
- ufl.inner(ufl.avg(2.0 * mu * epsilon(u)), jump(v, n)) * dS
- ufl.inner(jump(u, n), ufl.avg(2.0 * mu * epsilon(v))) * dS
+ (alpha_sipg_ns / ufl.avg(h)) * ufl.inner(jump(u, n), jump(v, n)) * ufl.avg(mu) * dS
)
un = ufl.dot(u, n) * n
vn = ufl.dot(v, n) * n
F_visc += (
- ufl.inner(2.0 * mu * epsilon(u) * n, vn) * ds
- ufl.inner(2.0 * mu * epsilon(v) * n, un) * ds
+ (alpha_bd / h) * mu_liquid * ufl.inner(un, vn) * ds
)
F_pressure = (
-ufl.inner(p, ufl.div(v)) * ufl.dx
+ ufl.inner(ufl.avg(p), ufl.jump(v, n)) * dS
+ (alpha_p / ufl.avg(h)) * ufl.dot(ufl.jump(p, n), ufl.jump(q, n)) * dS
)
F_continuity = (
-ufl.inner(ufl.div(u), q) * ufl.dx
+ ufl.inner(ufl.jump(u, n), ufl.avg(q)) * dS
+ ufl.inner(ufl.dot(u, n), q) * ds
)
F_gravity = -ufl.inner(rho(c) * g_vec, v) * ufl.dx
# Korteweg force (phase-field capillary term)
F_korteweg = (
- ufl.inner(mu_c * ufl.grad(c), v) * ufl.dx
- ufl.inner(ufl.avg(mu_c) * ufl.jump(c, n), ufl.avg(v)) * dS
)
return term_time + F_conv + F_visc + F_pressure + F_continuity + F_gravity + F_korteweg
def make_F_CH():
term_ch_time = ufl.inner((c - c_n) / dt, phi) * ufl.dx
# CH advection: lag only the upwind direction (consistent with NS)
lmbda_upw = ufl.conditional(ufl.gt(ufl.dot(u_n, n), 0), 1, 0)
c_uw = lmbda_upw("+") * c("+") + lmbda_upw("-") * c("-")
term_ch_adv_vol = -ufl.inner(ufl.outer(c, u), ufl.grad(phi)) * ufl.dx
term_ch_adv_surf = (
ufl.inner(ufl.dot(u, n)("+") * c_uw, phi("+")) * dS +
ufl.inner(ufl.dot(u, n)("-") * c_uw, phi("-")) * dS
)
term_ch_adv_bd = ufl.inner(ufl.dot(u_n, n) * lmbda_upw * c, phi) * ds
F_ch_adv = term_ch_adv_vol + term_ch_adv_surf + term_ch_adv_bd
F_ch_mu_diff = M * (
ufl.inner(ufl.grad(mu_c), ufl.grad(phi)) * ufl.dx
- ufl.inner(ufl.avg(ufl.grad(mu_c)), ufl.jump(phi, n)) * dS
- ufl.inner(ufl.jump(mu_c, n), ufl.avg(ufl.grad(phi))) * dS
+ (alpha_sipg_ch / ufl.avg(h)) * ufl.inner(ufl.jump(mu_c, n), ufl.jump(phi, n)) * dS
+ beta_ch * ufl.inner(ufl.jump(mu_c, n), ufl.jump(phi, n)) * dS
)
F_ch_c_lap = -lmbda * (
ufl.inner(ufl.grad(c), ufl.grad(psi)) * ufl.dx
- ufl.inner(ufl.avg(ufl.grad(c)), ufl.jump(psi, n)) * dS
- ufl.inner(ufl.jump(c, n), ufl.avg(ufl.grad(psi))) * dS
+ (alpha_sipg_ch / ufl.avg(h)) * ufl.inner(ufl.jump(c, n), ufl.jump(psi, n)) * dS
)
F_ch_mu = ufl.inner(mu_c, psi) * ufl.dx - (lmbda / eps**2) * ufl.inner(dfdc, psi) * ufl.dx
return term_ch_time + F_ch_adv + F_ch_mu_diff + F_ch_c_lap + F_ch_mu
F_total = make_F_NS() + make_F_CH()
J_total = ufl.derivative(F_total, w)
# ====================== Boundary conditions ======================
W_p, _ = W.sub(1).collapse()
def origin(x):
return np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 0.0))
dofs_p = fem.locate_dofs_geometrical((W.sub(1), W_p), origin)
bcs = []
if len(dofs_p[0]) > 0:
p_bc = fem.dirichletbc(default_scalar_type(0.0), dofs_p[0], W.sub(1))
bcs.append(p_bc)
# ====================== Solver settings ======================
petsc_options = {
"snes_type": "newtonls",
"snes_rtol": 1e-8,
"snes_atol": 1e-10,
"snes_stol": 0.0,
"snes_max_it": 50,
"snes_monitor": None,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
problem = NonlinearProblem(F_total, w, bcs=bcs, J=J_total, petsc_options_prefix="rising_bubble_dg_", petsc_options=petsc_options)
# ====================== Output settings ======================
W_u_out = fem.functionspace(domain, ("DG", k, (domain.geometry.dim,)))
W_p_out = fem.functionspace(domain, ("DG", k-1))
W_c_out = fem.functionspace(domain, ("DG", k-1))
u_out = fem.Function(W_u_out, name="Velocity")
p_out = fem.Function(W_p_out, name="Pressure")
c_out = fem.Function(W_c_out, name="Phase")
folder = Path("rising_bubble_DG_results_BDF1_lagged_dir")
folder.mkdir(exist_ok=True, parents=True)
u_writer = VTXWriter(domain.comm, folder / "velocity.bp", u_out)
p_writer = VTXWriter(domain.comm, folder / "pressure.bp", p_out)
c_writer = VTXWriter(domain.comm, folder / "phase.bp", c_out)
def update_output_functions():
u_out.interpolate(w.sub(0))
p_out.interpolate(w.sub(1))
c_out.interpolate(w.sub(2))
update_output_functions()
u_writer.write(0.0)
p_writer.write(0.0)
c_writer.write(0.0)
def compute_cfl(u_func, dt_val, h_char):
u_vec = u_func
u_mag = ufl.sqrt(ufl.dot(u_vec, u_vec))
u_max = fem.assemble_scalar(fem.form(u_mag * ufl.dx))
cfl = u_max * dt_val / h_char
return u_max, cfl
# ====================== Time marching loop ======================
T_end = 3.0
step = 0
output_frequency = 2
print("=== DG simulation started (BDF1 + lagged upwind) ===")
t = 0.0
while t < T_end:
t += dt_val
step += 1
_ = problem.solve()
converged_reason = problem.solver.getConvergedReason()
assert converged_reason > 0
num_iterations = problem.solver.getIterationNumber()
print(f"time {t} Step {step}: {converged_reason=} {num_iterations=}")
u_max, cfl = compute_cfl(w.sub(0), dt_val, mesh_size)
print(f"t={t:.4f} max|u|={u_max:.4e} CFL={cfl:.3f}")
w_n.x.array[:] = w.x.array[:]
u_n.interpolate(w_n.sub(0))
c_n.interpolate(w_n.sub(2))
if step % output_frequency == 0:
update_output_functions()
u_writer.write(t)
p_writer.write(t)
c_writer.write(t)
u_writer.close()
p_writer.close()
c_writer.close()
modify for the hysing2parameter
modify for the bubble rising
@nate, hello!, sorry to bother you. Since you’re an expert on DG, would you mind taking a look? Thanks a lot
Unfortunately I don’t have much experience here, and I don’t have time to read your code in detail. But based on what you’ve discussed in your post perhaps a simple upwind flux approximation is insufficient? This may be one of the problems that require careful stabilisation using a higher order, e.g., WENO scheme, or a sophisticated flux limiter.
This is really tough to do for unstructured grids. You could take a look at https://www.ocellaris.org/ to see if it can help you at all.