Infinite parallel apply

1

I have problem, i parallel my script and it freeze in apply:

if comm.rank == 0:
    print('A e 3')


# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu_2]
[bc.apply(A2) for bc in bcp_2]

if comm.rank == 0:
    print('A e 4')

It prints ‘A e 3’ and then freezes consuming CPU resources.
Full code:

u_n and u_n_2 - were calculated earlier
Using a single core, everything works well

mesh = Mesh()
with XDMFFile("T.xdmf") as infile:
    infile.read(mesh)

V = VectorFunctionSpace(mesh, "P", 1)
P = FunctionSpace(mesh, "P", 1)

dim = mesh.topology().dim()
if comm.rank == 0:
    print('Dim:',dim)

mvc = MeshValueCollection("size_t", mesh, dim-1) 
with XDMFFile("facet_mesh_T.xdmf") as infile:
    infile.read(mvc, "name_to_read")
facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

if comm.rank == 0:
    print('Allow')
    
u_n.set_allow_extrapolation(True)
u_n_2.set_allow_extrapolation(True)

bcp_inflow_top = DirichletBC(V, u_n_2, facet_domains, 63)
bcp_inflow_bottom = DirichletBC(V, u_n, facet_domains, 64)
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), facet_domains, 66)
bcp_outflow = DirichletBC(P, Constant(0), facet_domains, 65)

bcu_2 = [bcu_walls, bcp_inflow_top, bcp_inflow_bottom]
bcp_2 = [bcp_outflow]

if comm.rank == 0:
    print('A e')

T = 1            # final time
num_steps = 10000   # number of time steps
dt = T / num_steps # time step size
k  = Constant(dt)
rho = 1000
mu = 0.001

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(P)
q = TestFunction(P)

u_n_asd = Function(V)
u_  = Function(V)
p_n = Function(P)
p_  = Function(P)

if comm.rank == 0:
    print('A e 2')

# Define expressions used in variational forms
U  = 0.5*(u_n_asd + u)
n  = FacetNormal(mesh)
f_constant  = Constant((0, 0, 0))

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n_asd) / k, v)*dx \
   + rho*dot(dot(u_n_asd, nabla_grad(u_n_asd)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f_constant, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

if comm.rank == 0:
    print('A e 3')


# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu_2]
[bc.apply(A2) for bc in bcp_2]

if comm.rank == 0:
    print('A e 4')

t = 0
iter_n = 0
out = File("time/u.pvd")

for n in range(num_steps): 

    # Update current time
    t += dt
    iter_n += 1
    
    start_time = time.time()
    
    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu_2]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp_2]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
    endtime =  time.time() - start_time

    error = np.abs(np.array(u_n_asd.vector()) - np.array(u_.vector())).max()
    max_u = np.abs(np.array(u_n_asd.vector())).max()
    u_n_asd.assign(u_)
    p_n.assign(p_)
    
    if iter_n % 100 == 0:
        out << f_n, t
    
    if comm.rank == 0:
        print('t = %.4f: error = %.6g, max_u= %.6g endtime = %.6g' % (t, error, max_u, endtime))
    
File("pvd/u_real.pvd") << u_n_asd
File("sol/u_n_real.xml.gz") << u_n_asd

Please tell me how to fix this problem.

2

Could you try to make the errro reproducable with a built-in mesh and try to reduce the problem to a minimal code. It will then be easier to pinpoint what is wrong.

3

Yes, of course, I created a problem that first calculates the speed in a unit square on one grid, then applies the calculated speed as a boundary condition to the other grid.

from fenics import *
import numpy as np
from mpi4py import MPI
import time

comm = MPI.COMM_WORLD

T = 1.0           # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size
mu = 1             # kinematic viscosity
rho = 1            # density

# Create mesh and define function spaces
mesh = UnitSquareMesh(30, 30)
V = VectorFunctionSpace(mesh, 'P', 1)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
     rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    plot(u_)

    # Compute error
    u_e = Expression(('4*x[1]*(1.0 - x[1])', '0'), degree=2)
    u_e = interpolate(u_e, V)
    error = np.abs(np.array(u_e.vector()) - np.array(u_.vector())).max()
    if comm.rank == 0:
        print('max u:', np.array(u_.vector()).max())

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    
if comm.rank == 0:
    print('Two work')
    
# Create mesh and define function spaces
mesh = UnitSquareMesh(100, 100)
V = VectorFunctionSpace(mesh, 'P', 1)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

if comm.rank == 0:
    print('Define bc for two work')

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcu_inflow  = DirichletBC(V, u_n, inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outflow]


# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U   = 0.5*(u_n_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n_n) / k, v)*dx + \
     rho*dot(dot(u_n_n, nabla_grad(u_n_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

if comm.rank == 0:
    print('Assemble two work')

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

if comm.rank == 0:
    print('Apply two work')

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    plot(u_)

    # Compute error
    u_e = Expression(('4*x[1]*(1.0 - x[1])', '0'), degree=2)
    u_e = interpolate(u_e, V)
    error = np.abs(np.array(u_e.vector()) - np.array(u_.vector())).max()
    if comm.rank == 0:
        print('max u:', np.array(u_.vector()).max())

    # Update previous solution
    u_n_n.assign(u_)
    p_n.assign(p_)

Pay attention to the setting of the boundary conditions in the second task, it seems to me that this is an error, but I do not know how to apply the cut of the velocity field from one task to another. An example like this works correctly without parallels.

4

You should interpolate u_n from the first mesh to the second mesh, as shown below:


u_n.set_allow_extrapolation(True)
# Create mesh and define function spaces
mesh = UnitSquareMesh(100, 100)
V = VectorFunctionSpace(mesh, 'P', 1)
Q = FunctionSpace(mesh, 'P', 1)
u_n_2 = interpolate(u_n, V)
# Define boundaries
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

if comm.rank == 0:
    print('Define bc for two work')

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcu_inflow  = DirichletBC(V, u_n_2, inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outflow]

then your code should run in both serial and parallel.

Please note that you should put more effort into reducing your example before you post. For instance, you do not need 500 time steps (one should suffice), and you should stop the code at the point where you obtain an error message.
In general, I would suggest you have a look at the examples of minimal code examples in:

1 Like
5

If the second geometry is very different from the first, for example, we looked at a piece of pipe to establish a flow, and then applied our flow to a curved pipe at the entrance. In 3D geometry. How then to set the boundary condition?

6

As I have already suggested, you should rather solve a simple flow (say Stokes) on your complex geometry to get an inflow condition, instead of trying to interpolate solutions from vastly different grids onto each other.

1 Like