Here is my code and I want to run it in parallel, however when I run with
mpirun -n 80 python NS.py
it doesn’t work at all with no response after the output
Process 0: Solving linear variational problem.
Process 1: Solving linear variational problem.
Process 65: Solving linear variational problem.
Process 25: Solving linear variational problem.
Process 50: Solving linear variational problem.
Process 15: Solving linear variational problem.
Process 61: Solving linear variational problem.
Process 64: Solving linear variational problem.
Process 39: Solving linear variational problem.
Process 68: Solving linear variational problem.
Process 31: Solving linear variational problem.
Process 42: Solving linear variational problem.
Process 30: Solving linear variational problem.
Process 47: Solving linear variational problem.
Process 79: Solving linear variational problem.
Process 14: Solving linear variational problem.
Process 48: Solving linear variational problem.
Process 55: Solving linear variational problem.
Process 63: Solving linear variational problem.
Process 18: Solving linear variational problem.
Process 12: Solving linear variational problem.
Process 56: Solving linear variational problem.
Process 2: Solving linear variational problem.
Process 46: Solving linear variational problem.
Process 58: Solving linear variational problem.
Process 69: Solving linear variational problem.
Process 70: Solving linear variational problem.
Process 57: Solving linear variational problem.
Process 33: Solving linear variational problem.
Process 78: Solving linear variational problem.
Process 5: Solving linear variational problem.
Process 43: Solving linear variational problem.
Process 60: Solving linear variational problem.
Process 10: Solving linear variational problem.
Process 49: Solving linear variational problem.
Process 26: Solving linear variational problem.
Process 74: Solving linear variational problem.
Process 36: Solving linear variational problem.
Process 71: Solving linear variational problem.
Process 7: Solving linear variational problem.
Process 9: Solving linear variational problem.
Process 32: Solving linear variational problem.
Process 13: Solving linear variational problem.
Process 72: Solving linear variational problem.
Process 73: Solving linear variational problem.
Process 38: Solving linear variational problem.
Process 22: Solving linear variational problem.
Process 59: Solving linear variational problem.
Process 44: Solving linear variational problem.
Process 19: Solving linear variational problem.
Process 77: Solving linear variational problem.
Process 3: Solving linear variational problem.
Process 17: Solving linear variational problem.
Process 37: Solving linear variational problem.
Process 4: Solving linear variational problem.
Process 75: Solving linear variational problem.
Process 8: Solving linear variational problem.
Process 29: Solving linear variational problem.
Process 20: Solving linear variational problem.
Process 67: Solving linear variational problem.
Process 24: Solving linear variational problem.
Process 40: Solving linear variational problem.
Process 28: Solving linear variational problem.
Process 27: Solving linear variational problem.
Process 76: Solving linear variational problem.
Process 23: Solving linear variational problem.
Process 6: Solving linear variational problem.
Process 34: Solving linear variational problem.
Process 52: Solving linear variational problem.
Process 66: Solving linear variational problem.
Process 51: Solving linear variational problem.
Process 16: Solving linear variational problem.
Process 54: Solving linear variational problem.
Process 53: Solving linear variational problem.
Process 35: Solving linear variational problem.
Process 41: Solving linear variational problem.
Process 45: Solving linear variational problem.
Process 21: Solving linear variational problem.
Process 62: Solving linear variational problem.
Process 11: Solving linear variational problem.
My code is like this:
"""
This is an example for lid-driven cavity.
Apply Oseen-iteration for nonlinear term u·grad(u), i.e., linearize it as u_n·grad(u_{n+1}) in time step n+1.
For time derivative, we apply Back-Euler scheme.
"""
import os
os.environ['OMP_NUM_THREADS'] = '1'
from mpi4py import MPI
from fenics import *
from dolfin import *
import meshio
import matplotlib.pyplot as plt
import numpy as np
A =0.0035
B =0.16
D =8.2
E =0.64
n =0.2128
a0=88.4237
a1=58.4368
b1=129.496
a2=-57.2075
b2=84.2232
a3=-54.152
b3=-1.2649
a4=-18.8393
b4=-1.6359
a5=-22.6316
b5=-13.693
a6=-1.8443
b6=-12.8328
a7=-2.0967
b7=-4.6423
a8=-0.0095054
b8=-5.0506
w=7.9524
t = np.arange(0,1.6, 0.0001)
S=0.0014665
c0 = 1.106e+04
c1 = -1576
d1 = 1471
c2 = -937.8
d2 = -1187
c3 = 359.1
d3 = 17.34
c4 = 49.73
d4 = -50.38
c5 = 26.29
d5 = 59.15
c6 = 5.488
d6 = -13.89
c7 = -4.794
d7 = 18.28
c8 = 4.711
d8 = -9.013
w2 = 7.907
mpi_comm = MPI.comm_world
mpi_rank = mpi_comm.Get_rank()
mesh = Mesh(mpi_comm)
with XDMFFile(mpi_comm,"aorta_4branches.xdmf") as infile:
infile.read(mesh)
##
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(mpi_comm,"mf_aorta_4branches.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
boundary_markers=MeshFunction("size_t",mesh,mvc)
##
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(mpi_comm,"cf_aorta_4branches.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
domain_markers=MeshFunction("size_t",mesh,mvc)
ds_f = Measure("ds",domain=mesh, subdomain_data=boundary_markers)
dx_f = Measure("dx",domain=mesh, subdomain_data=domain_markers)
# Physical Constant
rho = 1060 # Fluid Density
mu = 0.0035 # Fluid Viscosity
T = 1.6 # Final Time
dt = 0.001 # Time Step
step = int(T/dt) + 1 # Number of Iterations
######################## Boundary Mark ########################
inlet_marker=1
bt_marker=2
lcca_marker=3
lsa_marker=4
descend_marker=5
wall_marker=6
domain_marker=7
inlet_direction=(-0.08563,-0.05819,0.99463)
bt_direction=(-0.67927,-0.08572,0.72886)
lcca_direction=(-0.20145,-0.05481,0.97796)
lsa_direction=(0.23336,0.83564,0.49723)
####################### Function Space #######################
# Taylor-Hood Elements
u_elem = VectorElement("CG", mesh.ufl_cell(), 2)
p_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
elem_f = MixedElement([u_elem, p_elem])
V_f = FunctionSpace(mesh, elem_f)
u, p = TrialFunctions(V_f)
v, q = TestFunctions(V_f)
NS_old = Function(V_f)
u_old, p_old = NS_old.split(True)
NS_present = Function(V_f)
u_fem, p_fem = NS_present.split(True)
parameters2={"linear_solver":"mumps",'preconditioner':'ilu'}
###################### Variational Form ######################
# Strain Rate Tensor
def epsilon(u):
return sym(nabla_grad(u))
n_f = FacetNormal(mesh)
# Source Term
f = Constant((0.0, 0.0,0.0))
A_inlet=0.0014665
A_bt=0.000211778
A_lcca=6.99274E-05
A_lsa=4.97932E-05
A_descend=0.0005121579
Qx_bt=1.593421206
Qx_lcca=0.915617678
Qx_lsa=0.772636993
####################### Time Iteration ######################
t = 0.0
#vtkfile_u = File('Example/velocity.pvd')
#vtkfile_p = File('Example/pressure.pvd')
xdmffile_u = XDMFFile(mpi_comm, 'results/velocity.xdmf')
xdmffile_p = XDMFFile(mpi_comm, 'results/pressure.xdmf')
with open('results.txt', 'w') as file:
file.write('t,inlet_p,bt_p,lcca_p,lsa_p,descend_p\n')
for n in range(step):
t += dt
print("** Solving tentative velocity step...\n" )
v_inlet_constant = (a0 + a1*np.cos(t*w) + b1*np.sin(t*w) +a2*np.cos(2*t*w) + b2*np.sin(2*t*w) + a3*np.cos(3*t*w) + b3*np.sin(3*t*w) + a4*np.cos(4*t*w) + b4*np.sin(4*t*w) + a5*np.cos(5*t*w) + b5*np.sin(5*t*w) + a6*np.cos(6*t*w) + b6*np.sin(6*t*w) + a7*np.cos(7*t*w) + b7*np.sin(7*t*w) + a8*np.cos(8*t*w) + b8*np.sin(8*t*w))/S*1e-6
v_bt = Constant(tuple([element * Qx_bt*v_inlet_constant for element in bt_direction]))
v_lcca = Constant(tuple([element * Qx_lcca*v_inlet_constant for element in lcca_direction]))
v_lsa = Constant(tuple([element * Qx_lsa*v_inlet_constant for element in lsa_direction]))
v_inlet = Constant(tuple([element*v_inlet_constant for element in inlet_direction]))
p_descend= c0 + c1*np.cos(t*w2) + d1*np.sin(t*w2) + c2*np.cos(2*t*w2) + d2*np.sin(2*t*w2) + c3*np.cos(3*t*w2) + d3*np.sin(3*t*w2) +c4*np.cos(4*t*w2) + d4*np.sin(4*t*w2) + c5*np.cos(5*t*w2) + d5*np.sin(5*t*w2) +c6*np.cos(6*t*w2) + d6*np.sin(6*t*w2) + c7*np.cos(7*t*w2) + d7*np.sin(7*t*w2) +c8*np.cos(8*t*w2) + d8*np.sin(8*t*w2);
bcu_inlet=DirichletBC(V_f.sub(0),v_inlet,boundary_markers,inlet_marker)
bcu_bt=DirichletBC(V_f.sub(0),v_bt,boundary_markers,bt_marker)
bcu_lcca=DirichletBC(V_f.sub(0),v_lcca,boundary_markers,lcca_marker)
bcu_lsa=DirichletBC(V_f.sub(0),v_lsa,boundary_markers,lsa_marker)
bcu_wall=DirichletBC(V_f.sub(0),Constant((0.0,0.0,0.0)),boundary_markers,wall_marker)
bc = [bcu_inlet,bcu_bt,bcu_lcca,bcu_lsa,bcu_wall]
# Solve NS equation
# Variational Form
a_f = rho*dot(u, v)*dx \
+ 2*dt*mu*inner(epsilon(u), epsilon(v))*dx \
+ dt*rho*dot(dot(u_old, nabla_grad(u)),v)*dx\
- dt*p*div(v)*dx \
+ dt*q*div(u)*dx
L_f = dt*dot(f, v)*dx \
+rho*dot(u_old,v)*dx \
- dt*dot(Constant(p_descend)*n_f, v)*ds_f(5)
solve(a_f == L_f, NS_present, bc,solver_parameters=parameters2)
u_fem, p_fem = NS_present.split(True)
if n % 2 == 0:
#vtkfile_u << (u_fem, t)
#vtkfile_p << (p_fem, t)
# Create XDMF files for visualization output
u_n.assign(u_fem)
p_n.assign(p_fem)
average_p_inlet = assemble(p_fem*ds_f(1))/assemble(Constant(A_inlet)*ds_f(1))
average_p_bt = assemble(p_fem*ds_f(2))/assemble(Constant(A_bt)*ds_f(2))
average_p_lcca = assemble(p_fem*ds_f(3))/assemble(Constant(A_lcca)*ds_f(3))
average_p_lsa = assemble(p_fem*ds_f(4))/assemble(Constant(A_lsa)*ds_f(4))
average_p_descend = assemble(p_fem*ds_f(5))/assemble(Constant(A_descend)*ds_f(5))
with open('results.txt', 'w') as file:
file.write(f'{t},{average_p_inlet},{average_p_bt},{average_p_lcca},{average_p_lsa},{average_p_descend}\n')
# Update u_old, p_old
assign(u_old, u_fem)
assign(p_old, p_fem)
print(f"------Step {n} ------")
print("Computation Complete!")
Could anyone help me solve the problem?