Hello,
I am trying to run a multiphase flow problem in Dolfinx using the momentum, gas phase transport, and continuity equations. When solving the transient problem, the residual seems to change after the first iteration, but then does not change and ultimately does not converge.
Are there any steps I can take to fix the static residual problem? I changed the solver settings a few times, which changed the static residual from nan to a finite value, but still does not update.
Here is my code:
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx
from dolfinx.fem import Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological, set_bc
from dolfinx import default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx import mesh
from basix.ufl import element, mixed_element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from dolfinx.mesh import GhostMode, locate_entities_boundary
# Define domain size
width, height = 0.5, 1.0
inlet_width = 0.2
nx, ny = 100, 100
# Create mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0], [width, height]], [nx, ny], mesh.CellType.triangle)
# Inlet: middle portion of bottom boundary
def inlet_boundary(x):
return np.logical_and(
np.isclose(x[1], 0.0), # bottom edge
np.logical_and(
x[0] > (width / 2 - inlet_width),
x[0] < (width / 2 + inlet_width)
)
)
# Walls: left and right edges
def wall_boundary(x):
return np.logical_or(
np.isclose(x[0], 0.0), # left
np.isclose(x[0], width) # right
)
# No-slip top boundary
def top_noslip_boundary(x):
return np.isclose(x[1], height) # top
with XDMFFile(domain.comm, "domain.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
t = 0 # Initial time
final_t = 2 # Final time
dt = Constant(domain, 1/1000)
D_gc = 0.01 # Drift diffusion coeff
u_slip = Constant(domain, np.array((0.0, 1e-2), dtype=np.float64))
p_ref = 1e5
M_hydrogen = 2.02e-3 # kg/mol
T = 500
R = 8.31
rho_l = 2729.3 - 0.73 * T
mg_l = Constant(domain, 1e-5)
mu_l = Constant(domain, 4E-4 * np.exp(4170 / T))
g = Constant(domain,[0.0, -9.81]) # Example for 2D gravity vector
# Define the elements
V_ele = element("CG", domain.topology.cell_name(), 2, shape=(domain.geometry.dim, ))
Q_ele = element("CG", domain.basix_cell(), 1) # Continuous Lagrange, degree 1
# Create a mixed function space
W = functionspace(domain, mixed_element([V_ele, Q_ele, Q_ele]))
upphi = Function(W)
upphi_n = Function(W)
vu, vp, vphi = ufl.TestFunctions(W)
u_l, p, phi_g = ufl.split(upphi)
u_l_n, p_n, phi_g_n = ufl.split(upphi_n)
phi_l = 1 - phi_g
phi_l_n = 1 - phi_g_n
# Definign sub spaces
V_sub, _ = W.sub(0).collapse()
Q1_sub, _ = W.sub(1).collapse()
Q2_sub, _ = W.sub(2).collapse()
# Inlet velocity conditions
inlet_facets = locate_entities_boundary(domain, 1, inlet_boundary)
# y-component = 1.0
u_inlety = Constant(domain, 1.0)
bcu_inlet1 = dirichletbc(u_inlety,
locate_dofs_topological(V_sub, 1, inlet_facets),
W.sub(0).sub(1))
# x-component = 0.0
u_inletx = Constant(domain, 0.0)
bcu_inlet2 = dirichletbc(u_inletx,
locate_dofs_topological(V_sub, 1, inlet_facets),
W.sub(0).sub(0))
# No-slip boundary condition for velocity
noslip = np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType)
wall_facets = locate_entities_boundary(domain, 1, wall_boundary)
bc_wall = dirichletbc(noslip,
locate_dofs_topological(V_sub, 1, wall_facets),
V_sub)
top_facets = locate_entities_boundary(domain, 1, top_noslip_boundary)
bc_top = dirichletbc(noslip,
locate_dofs_topological(V_sub, 1, top_facets),
V_sub)
phi_val = Constant(domain, 0.1)
# Locate inlet dofs for phi_g
inlet_facets = locate_entities_boundary(domain, 1, inlet_boundary)
inlet_dofs_phi = locate_dofs_topological(Q2_sub, 1, inlet_facets)
# Apply DirichletBC
bc_phi_inlet = dirichletbc(phi_val, inlet_dofs_phi, Q2_sub)
# Combine all velocity BCs
bcs = [bcu_inlet1, bcu_inlet2, bc_wall, bc_top, bc_phi_inlet]
u_drift = -D_gc * ufl.grad(phi_g) / (phi_g + 1e-8)
u_g = u_l
rho_g = (p_ref + p) * M_hydrogen / (R * T)
continuity = (
rho_l * (phi_l - phi_l_n) + rho_g * (phi_g - phi_g_n)
) / dt * vp * ufl.dx + ufl.inner(
ufl.div(rho_l * phi_l * u_l + rho_g * phi_g * u_g), vp
) * ufl.dx
gas_phase_transport = (
(rho_g * (phi_g - phi_g_n)) / dt * vphi * ufl.dx
+ ufl.div(rho_g * phi_g * u_l) * vphi * ufl.dx
+ ufl.inner(D_gc * rho_g * ufl.grad(phi_g), ufl.grad(vphi)) * ufl.dx
+ mg_l * vphi * ufl.dx
)
momentum = (
phi_l * rho_l * ufl.inner((u_l - u_l_n) / dt, vu) * ufl.dx
+ phi_l * mu_l * ufl.inner(ufl.grad(u_l), ufl.grad(vu)) * ufl.dx
+ phi_l * rho_l * ufl.inner(ufl.grad(u_l), ufl.grad(vu)) * ufl.dx # Fixed
+ p * ufl.div(vu) * ufl.dx
- phi_l * rho_l * ufl.inner(g, vu) * ufl.dx
)
F = momentum + continuity + gas_phase_transport
problem = NonlinearProblem(F, upphi, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-8
solver.report = True
# PETSc solver options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-10
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
import os
from dolfinx.io import VTXWriter
# Save directory
save_dir = "results/"
os.makedirs(save_dir, exist_ok=True)
# Initialize VTX writers
vtx_u = VTXWriter(MPI.COMM_WORLD, f"{save_dir}u.bp", [Function(V_sub)])
vtx_p = VTXWriter(MPI.COMM_WORLD, f"{save_dir}p.bp", [Function(Q1_sub)])
vtx_phi = VTXWriter(MPI.COMM_WORLD, f"{save_dir}phi_g.bp", [Function(Q2_sub)])
vtx_u.write(t)
vtx_p.write(t)
vtx_phi.write(t)
# Time loop
while t < float(final_t):
n, converged = solver.solve(upphi)
assert converged
u, p, phi_g = upphi.split()
print("phi_g min/max:", phi_g.x.array.min(), phi_g.x.array.max())
vtx_u.write(t)
vtx_p.write(t)
vtx_phi.write(t)
print(f"Time = {t:.2f}, iterations = {n}, u[0] = {u.x.array[0]:.3e}")
upphi_n.x.array[:] = upphi.x.array
t += float(dt)
vtx_u.close()
vtx_p.close()
vtx_phi.close()
and the unchanging output until it reaches the max number of iterations:
[2025-04-24 14:19:30.102] [info] Newton iteration 0: r (abs) = 7142.13717233418 (tol = 1e-10), r (rel) = inf (tol = 1e-08)
[2025-04-24 14:19:30.331] [info] PETSc Krylov solver starting to solve system.
[2025-04-24 14:19:31.013] [info] Newton iteration 1: r (abs) = 41517.80332674378 (tol = 1e-10), r (rel) = 614.353097327223 (tol = 1e-08)
[2025-04-24 14:19:31.234] [info] PETSc Krylov solver starting to solve system.
[2025-04-24 14:19:31.445] [info] Newton iteration 2: r (abs) = 41517.80332674378 (tol = 1e-10), r (rel) = 614.353097327223 (tol = 1e-08)
[2025-04-24 14:19:31.663] [info] PETSc Krylov solver starting to solve system.
[2025-04-24 14:19:31.865] [info] Newton iteration 3: r (abs) = 41517.80332674378 (tol = 1e-10), r (rel) = 614.353097327223 (tol = 1e-08)
[2025-04-24 14:19:32.081] [info] PETSc Krylov solver starting to solve system.
[2025-04-24 14:19:32.283] [info] Newton iteration 4: r (abs) = 41517.80332674378 (tol = 1e-10), r (rel) = 614.353097327223 (tol = 1e-08)