My dear, thank you very much for your previous assistance. The issue regarding the brain mesh has been resolved. However, the simulation failed due to the excessively dense mesh. To facilitate a better understanding of my program framework, I am providing the code supplied by Solving 3D Stokes problem is slower on dolfinx than dolfin, which follows the same program logic. Moreover, I have fixed a little bit of pressure. When executing with “mpirun -n 24”, the following error message occurs. I look forward to your reply.
import os
os.environ['OMP_NUM_THREADS'] = '1'
from mpi4py import MPI
comm = MPI.COMM_WORLD
import numpy as np
import time
import ufl, basix
from dolfinx import fem, io, mesh
from ufl import grad, div, nabla_grad, dx, inner, dot
from dolfinx.fem.petsc import LinearProblem
# build mesh
msh = mesh.create_box(MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([2, 1, 0.3])],
[60, 60, 60], cell_type=mesh.CellType.tetrahedron)
# define element
Pv = basix.ufl.element("CG", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = basix.ufl.element("CG", msh.basix_cell(), 1)
TH = basix.ufl.mixed_element([Pv, P1])
# create functionspace
W = fem.functionspace(msh, TH)
V = fem.functionspace(msh, Pv)
# define functions
dU = ufl.TrialFunction(W)
U = fem.Function(W)
velocity = fem.Function(V)
f = fem.Function(V)
# set boundary conditions
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
def LEFT(xx):
return np.isclose(xx[0], 0)
def RIGHT(xx):
return np.isclose(xx[0], 2)
def corner(x):
return np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))
boundary_surface = mesh.exterior_facet_indices(msh.topology)
left_surface = mesh.locate_entities_boundary(msh, fdim, LEFT)
right_surface = mesh.locate_entities_boundary(msh, fdim, RIGHT)
side = np.concatenate((left_surface, right_surface))
other_surface = np.setdiff1d(boundary_surface, side)
W_v_sub, _ = W.sub(0).collapse()
Q0 = W.sub(1)
V_in = fem.Function(W_v_sub)
V_noslip = fem.Function(W_v_sub)
V_noslip.x.array[:] = 0
V_in.sub(0).interpolate(lambda x: 1 * 4 * x[1] * (1 - x[1]) / (1) ** 2 * 4 * x[2] * (0.3 - x[2]) / (0.3) ** 2)
bcu = fem.dirichletbc(V_in, fem.locate_dofs_topological([W.sub(0), W.sub(0).collapse()[0]], fdim, left_surface),
W.sub(0))
bcu1 = fem.dirichletbc(V_in, fem.locate_dofs_topological([W.sub(0), W.sub(0).collapse()[0]], fdim, right_surface),
W.sub(0))
bcu2 = fem.dirichletbc(V_noslip, fem.locate_dofs_topological([W.sub(0), W.sub(0).collapse()[0]], fdim, other_surface),
W.sub(0))
# fix pressure in one corner
corner_vertex = mesh.locate_entities_boundary(msh, 0, corner)
dofs = fem.locate_dofs_topological(Q0, 0, corner_vertex)
bc3 = fem.dirichletbc(fem.Constant(msh, 0.0), dofs, Q0)
bcs = [bcu, bcu1, bcu2]
# define the linear problem
h = 2 * ufl.Circumradius(msh)
v, q = ufl.TestFunctions(W)
u, p = ufl.split(U)
du, dp = ufl.split(dU)
F_a = (
-inner(dp, div(v)) * dx
+ inner(2 * ufl.sym(nabla_grad(du)), ufl.sym(nabla_grad(v))) * dx
+ inner(div(du), q) * dx
+ 0.01 * h * h * inner(grad(dp), grad(q)) * dx
)
F_L = (
inner(f, v) * dx
)
problem_U = LinearProblem(F_a, F_L, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
# solve
ta = time.time()
Uh = problem_U.solve()
tb = time.time()
if MPI.COMM_WORLD.rank == 0:
print("Time: ", tb - ta)
Error
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= PID 8415 RUNNING AT LAPTOP-BBJKVL4R
= EXIT CODE: 9
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Killed (signal 9)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions
Version information: DOLFINx version: 0.9.0