Hello! I’m new to DOLFINx and FEniCS, so I hope this isn’t an overly simple question.
I’m trying to write FEniCS code to solve a Conjugate Heat Transfer problem on a 2D domain composed of air, a PCB, and components.
Since I’m dealing with forced convection, I’m trying to decouple the energy equation, which I plan to solve after the fluid problem. However, I’m currently unable to solve the steady-state Navier-Stokes problem.
I’ve written the following code (see below), but when I open the results in ParaView, I see solutions where entire fields are zero.
The console output says it solved successfully:
Info : Reading 'results/circuit_mesh_2d.msh'...
Info : 68 entities
Info : 13869 nodes
Info : 14502 elements
Info : Done reading 'results/circuit_mesh_2d.msh'
solved successfully
p wrote successfully
u wrote successfully
This is my full code:
import os
from pathlib import Path
from typing import Final, List
import dolfinx
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, io
from dolfinx.fem import Constant, Function, dirichletbc, locate_dofs_topological
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py import PETSc
COMM: Final[MPI.Intracomm] = MPI.COMM_WORLD
FLUID_TAG: Final[int] = 1
PCB_TAG: Final[int] = 2
COMPONENT_BASE_TAG: Final[int] = 3
COMPONENT_TAGS: Final[List[int]] = [3, 4, 5, 6, 7]
INLET_TAG: Final[int] = 11
OUTLET_TAG: Final[int] = 12
WALLS_TAG: Final[int] = 13
SOLID_TAG: Final[int] = 14 # Fluid-solid interface (no-slip)
NU: Final[float] = 1.48e-5 # Kinematic viscosity [m²/s]
VEL_INLET: Final[float] = 0.85 # Inlet velocity (x-direction) [m/s]
MESH_FILE: Final[str] = "circuit_mesh_2d.msh"
OUTPUT_FOLDER: Final[str] = "results"
k = 1 # Polynomial degree
def log(msg: str) -> None:
"""
Print message only on MPI rank 0.
Args:
msg: Message to print.
"""
if COMM.rank == 0:
print(msg, flush=True)
def inlet_velocity(x: np.ndarray) -> np.ndarray:
"""
Defines a constant velocity profile at an inlet boundary (x-direction only).
Returns a velocity array of shape (num_components, num_points):
[VEL_INLET, 0.0, ...]
"""
ux = np.full_like(x[0], VEL_INLET)
return np.array([ux, np.full_like(x[0], 0.0)])
def define_bcs(msh, ft, ct, W):
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
wall_facets = ft.find(WALLS_TAG)
inlet_facets = ft.find(INLET_TAG)
outlet_facets = ft.find(OUTLET_TAG)
solid_facets = ft.find(SOLID_TAG)
combined_indices = np.concatenate(
[ct.find(tag) for tag in COMPONENT_TAGS + [PCB_TAG]]
)
# DOFS
dofs_wall = locate_dofs_topological(W.sub(0), msh.topology.dim - 1, wall_facets)
dofs_inlet = locate_dofs_topological(W.sub(0), msh.topology.dim - 1, inlet_facets)
dofs_solid = locate_dofs_topological(W.sub(0), msh.topology.dim - 1, solid_facets)
dofs_solid_cell = locate_dofs_topological(
W.sub(0), msh.topology.dim, combined_indices
)
dofs_outlet = locate_dofs_topological(W.sub(1), msh.topology.dim - 1, outlet_facets)
# BC
zero_vector_func = Function(W.sub(0).collapse()[0])
zero_vector_func.x.array[:] = 0.0
bc_noslip = dirichletbc(zero_vector_func, dofs_wall)
bc_solid_u = dirichletbc(zero_vector_func, dofs_solid)
u_inlet_func = Function(W.sub(0).collapse()[0])
u_inlet_func.interpolate(inlet_velocity)
bc_inlet_x = dirichletbc(u_inlet_func, dofs_inlet)
bc_outlet_p = dirichletbc(
Constant(msh, PETSc.ScalarType(0.0)), dofs_outlet, W.sub(1)
)
bc_solid_cell_u = dirichletbc(zero_vector_func, dofs_solid_cell)
bcs = [bc_noslip, bc_inlet_x, bc_solid_u] # , bc_solid_cell_u]
return bcs, [bc_outlet_p]
folder = Path(OUTPUT_FOLDER)
folder.mkdir(exist_ok=True, parents=True)
mesh_file = os.path.join(folder, MESH_FILE)
msh = io.gmsh.read_from_msh(mesh_file, COMM, gdim=2)
msh, ct, ft = msh.mesh, msh.cell_tags, msh.facet_tags
# Constants and function spaces for the Stokes/Navier-Stokes problem
nu_const = fem.Constant(msh, dolfinx.default_scalar_type(NU))
gdim = msh.geometry.dim
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
# Function spaces for the velocity and for the pressure
v_cg2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
s_cg1 = element("Lagrange", msh.basix_cell(), 1)
V = fem.functionspace(msh, v_cg2)
Q = fem.functionspace(msh, s_cg1)
w_mixed_element = mixed_element([v_cg2, s_cg1])
W = fem.functionspace(msh, w_mixed_element)
# BC
bcs = define_bcs(msh, ft, ct, W)
all_bcs = bcs[0] + bcs[1]
w = Function(W)
(u, p) = w.split()
(v, q) = ufl.TestFunctions(W)
F = (
nu_const * ufl.inner(ufl.grad(u), ufl.grad(v)) # Viscous term
- ufl.div(v) * p # Pressure term (divergence of test function)
+ q * ufl.div(u) # Continuity term
+ ufl.inner(ufl.grad(u) * u, v) # Advection term (Non-linear)
) * ufl.dx()
petsc_options = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_linesearch_type": "bt",
}
pde = fem.petsc.NonlinearProblem(
F, w, bcs=all_bcs, petsc_options_prefix="nls", petsc_options=petsc_options
)
w_o = pde.solve()
assert pde.solver.getConvergedReason() > 0
print("solved successfully")
p_file = io.VTKFile(msh.comm, folder / "p.pvd", "w")
# p_file.write_mesh(msh)
p_collapsed_func = w_o.sub(1).collapse()
p_file.write_function(p_collapsed_func, 0.0)
p_file.close()
print("p wrote successfully")
u_file = io.VTKFile(msh.comm, folder / "u.pvd", "w")
# u_file.write_mesh(msh)
u_collapsed_func = w_o.sub(0).collapse()
element_v1 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
V1 = fem.functionspace(msh, element_v1)
u_save = fem.Function(V1)
u_save.interpolate(u_collapsed_func)
u_file.write_function(u_save, 0.0)
u_file.close()
print("u wrote successfully")
