Hello,
I am trying to set up a 2D Bubbly flow system in Dolfinx, but having trouble setting up the solver settings. I’m attempting to solve for each variable (liquid velocity, pressure, and gas fraction) by following the 2D channel example, but am unsure if my procedure is correct.
- If my weak form equations are equal to 0, how can I solve the system with all three variables/steps (no longer bilinear)?
- Is it simpler to solve one mixed function space with all three variables?
I’ve attached code snippets below and would appreciate any guidance - thank you!
Mesh setup:
gmsh.initialize()
# Mesh definition
width, height = 0.5, 1
nx, ny = 100, 100
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
domain = gmsh.model.occ.addRectangle(0, 0, 0, width, height, tag=1)
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 = 2, 3, 4
# Domain size
inlet_width = 0.025 # Half-width of inlet (e.g., 5% of width)
inlet, outlet, walls = [], [], []
if mesh_comm.rank == model_rank:
# Get all 1D boundary entities from 2D geometry
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
dim, tag = boundary
com = gmsh.model.occ.getCenterOfMass(dim, tag)
# Check for INLET (centered at bottom middle)
if np.isclose(com[1], 0.0):
if (width / 2 - inlet_width) < com[0] < (width / 2 + inlet_width):
inlet.append(tag)
else:
walls.append(tag)
# Check for OUTLET (top boundary)
elif np.isclose(com[1], height):
outlet.append(tag)
# All other sides are WALLS
else:
walls.append(tag)
# Tag physical groups
gmsh.model.addPhysicalGroup(1, inlet, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outlet, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
gmsh.model.mesh.generate(gdim)
domain, cell_markers, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ft.name = "Facet markers"
Here is how I define my function spaces for the three variables:
v_cg2 = element("Lagrange", domain.topology.cell_name(), 2, shape=(domain.geometry.dim, ))
s_cg1 = element("Lagrange", domain.topology.cell_name(), 1)
V = functionspace(domain, v_cg2)
Q = functionspace(domain, s_cg1)
u_l = Function(V)
p = Function(Q)
phi_g = Function(Q)
u_l_n = Function(V)
p_n = Function(Q)
phi_g_n = Function(Q)
vu = ufl.TestFunction(V)
vp, vphi = ufl.TestFunction(Q), ufl.TestFunction(Q)
Here are my variational equations, where I moved all terms to one side (LHS=0):
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
)
My attempt to assemble and solve the transient problem:
a1 = form(lhs(momentum))
L1 = form(rhs(momentum))
A1 = assemble_matrix(a1, bcs=bcu)
b1 = create_vector(L1)
a2 = form(lhs(continuity))
L2 = form(rhs(continuity))
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)
a3 = form(lhs(gas_phase_transport))
L3 = form(rhs(gas_phase_transport))
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)
# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")
# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
from pathlib import Path
from dolfinx.io import VTXWriter
from petsc4py import PETSc
from tqdm.autonotebook import tqdm
# Output setup
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, str(folder / "u_l.bp"), [u_l], engine="BP4")
vtx_p = VTXWriter(mesh.comm, str(folder / "p.bp"), [p], engine="BP4")
vtx_phi = VTXWriter(mesh.comm, str(folder / "phi_g.bp"), [phi_g], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
vtx_phi.write(t)
progress = tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):
progress.update(1)
t += dt
inlet_velocity.t = t
u_inlet.interpolate(inlet_velocity)
# --- Step 1: Momentum equation for u_l ---
A1.zeroEntries()
assemble_matrix(A1, a1, bcs=bcu)
A1.assemble()
with b1.localForm() as loc:
loc.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_l.x.petsc_vec)
u_l.x.scatter_forward()
# --- Step 2: Pressure equation for p ---
with b2.localForm() as loc:
loc.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, p.x.petsc_vec)
p.x.scatter_forward()
# --- Step 3: Gas transport equation for phi_g ---
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, phi_g.x.petsc_vec)
phi_g.x.scatter_forward()
But immediately run into the problem:
Traceback (most recent call last):
File "/mnt/c/Users/ckhur/sparge_test_v2.py", line 187, in <module>
A1 = assemble_matrix(a1, bcs=bcu)
File "/home/ckhurana/anaconda3/envs/fenicsx-env/lib/python3.13/functools.py", line 934, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
File "/home/ckhurana/anaconda3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/petsc.py", line 444, in assemble_matrix
A = _cpp.fem.petsc.create_matrix(a._cpp_object)
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.