Dear Community
I am using the following code to simulate Stokes flow in a 3D spherical annulus, and calculate the drag force. For small values of the external radius (Re), the code runs fine when executed using mpirun:
mpirun -np 36 python3 calc_vel_3d_stokes.py ulimit -s unlimited ulimit -c unlimited
However, as I go to larger values of Re, the execution terminates with the following error messages:
--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
Read -1, expected 181549016, errno = 3
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 0 on node exited on signal 9 (Killed).
--------------------------------------------------------------------------
The relevant FENICS script is given below
import numpy as np
from dolfin import *
import numpy as np
from mshr import *
Re = 10.
Ri = 1.
bound = Sphere(Point(0., 0., 0.), Re)
ball = Sphere(Point(0., 0., 0.), Ri)
domain = bound - ball
# Making a mark dictionary
# Note: the values should be UNIQUE identifiers.
mark = {"generic": 0,
"inner": 1,
"outer": 2}
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomains.set_all(mark["generic"])
r_mid = Ri + (Re - Ri) / 2
class inner_boun(SubDomain):
def inside(self, x, on_boundary):
return sqrt(x[0]**2 + x[1]**2 + x[2]**2) < r_mid and on_boundary
class outer_boun(SubDomain):
def inside(self, x, on_boundary):
return sqrt(x[0]**2 + x[1]**2 + x[2]**2) > r_mid and on_boundary
## Marking boundaries
in_boun = inner_boun()
in_boun.mark(subdomains, mark["inner"])
out_boun = outer_boun()
out_boun.mark(subdomains, mark["outer"])
geom_name="geometry.pvd"
File(geom_name) << subdomains
## Define function spaces
V = VectorElement("CG", mesh.ufl_cell(), 2)
P = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*P)
## Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains) # Volume integration
ds = Measure("ds", domain=mesh, subdomain_data=subdomains) # Surface integration
## Applying Dirichlet BCs at the inner and outer boundaries
# No-slip at inner boundary
# Constant velocity at outer boundary
u_outer_boun = Constant((0.0, 0.0, U_0))
bc_outer_boun = DirichletBC(W.sub(0), u_outer_boun, subdomains, mark["outer"])
noslip = Constant((0.0, 0.0, 0.0))
bc_inner_boun = DirichletBC(W.sub(0), noslip, subdomains, mark["inner"])
bc_p = DirichletBC(W.sub(1), 0, "x[0] < -1.0 + DOLFIN_EPS", "pointwise")
bcs = [bc_outer_boun,bc_inner_boun,bc_p]
# Body force:
force = Constant((0.0, 0.0, 0.0))
##Setting up variational formulation
a = inner(grad(u), grad(v))*dx - p*div(v)*dx + q*div(u)*dx L = inner(force, v) * dx
### Solve variational problem
w = Function(W)
solve(a == L, w, bcs)
### Split using deep copy
(u, p) = w.split(True)
ufile = XDMFFile("velocity.xdmf")
ufile.write(u)
pfile = XDMFFile("pressure.xdmf")
pfile.write(p)
# Traction calculation on inner boundary
# Viscous stress
# sym returns the symmetric part of a matrix
stress_visc = 2*sym(grad(u))
# Total stress
stress = -p*Identity(3) + stress_visc
# Surface normal
n = FacetNormal(mesh)
# Imposed flow direction:
n_flow = Constant((0.0,0.0,1.0))
n_flow = project(n_flow, W.sub(0).collapse())
# Compute traction vector
traction = dot(stress, n)
# Integrate decomposed traction to get total F_drag and F_lift
F_drag = assemble(dot(traction, n_flow) * ds(mark["inner"]))
print ("F_drag = ")
print (F_drag)
# Writing the result of drag calculation to file
fname="Ro{}_msh{}_f_drag.dat".format(Re,mesh_resolution)
e_file=open(fname,"w")
print("%.12f\t%d\t%.12f" %(Re,mesh_resolution,F_drag), file=e_file)
e_file.close()
Could you please give some suggestions for troubleshooting the issue?
Thank You
Warm Regards