Parallel solve from solve method

TherealSova · May 27, 2021, 11:35am

Hello, i am begginer in parallel fenics technology.
I have Navie-Stokes solver for 3D geometry, and i want solve in parallel. How i can parallel procedure solve? Is it possible to solve matrices in parallel when calling solve method?
My code with solve method:

...
#num_steps
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')
    
    max_u = np.abs(np.array(u_n.vector()) - np.array(u_.vector())).max()
    
    u_n.assign(u_)
    p_n.assign(p_)

    print('t = %.2f: error = %.6g' % (t, max_u))

I can attach the complete code of the problem that works for non-parallel computations

dokken · May 27, 2021, 2:47pm

FEniCS does by default run in parallel using mpirun -n num_procs program.py. The solvers interfaced are the PETSc solver suite, see PETSc: Summary of Sparse Linear Solvers Available from PETSc
which gives you an overview of which solvers run in parallel or not.

TherealSova · May 27, 2021, 3:31pm

What is mean X in a cell in the table?

dokken · May 27, 2021, 3:38pm

X means that it is supported (i.e. it supports either parallel or complex depending on the column

TherealSova · May 27, 2021, 3:41pm

Thanks, for you answer

TherealSova · June 1, 2021, 3:53pm

I make primitive fenics script, and i want use parallel solver:
I start with this command:

mpirun -n 8 3D_tube_script.py

But i see thih problem:

mpirun was unable to find the specified executable file, and therefore
did not launch the job.  This error was first reported for process
rank 0; it may have occurred for other processes as well.

NOTE: A common cause for this error is misspelling a mpirun command
      line parameter option (remember that mpirun interprets the first
      unrecognized command line token as the executable).

Node:       fenics
Executable: 3D_tube_script.py

My script:

import numpy as np
import meshio
msh = meshio.read("gmsh/3D_channel.msh")
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

tetra_mesh = create_mesh(msh, "tetra", prune_z=True)
meshio.write("mesh.xdmf", tetra_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("facet_mesh.xdmf", triangle_mesh)

print("Normal read!")

from dolfin import * 

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool((x[2] < DOLFIN_EPS or x[2] > (1.0 - DOLFIN_EPS)) \
                    and on_boundary)


    def map(self, x, y):
        y[0] = x[0] 
        y[1] = x[1]
        y[2] = x[2] - 1

pbc = PeriodicBoundary()  
V = VectorFunctionSpace(mesh, "P", 2, constrained_domain=pbc)
P = FunctionSpace(mesh, "P", 1)

dim = mesh.topology().dim()
print('Dim:',dim)

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

bcp_inflow = DirichletBC(P, Constant(8), facet_domains, 6)
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), facet_domains, 8)
bcp_outflow = DirichletBC(P, Constant(0), facet_domains, 7)

bcu = [bcu_walls]
bcp = [bcp_outflow, bcp_inflow]

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

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

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

# Define expressions used in variational forms
U  = 0.5*(u_n + 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) / 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_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

print("Go to assemble!")
# 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]

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')
    
    error = np.abs(np.array(u_n.vector()) - np.array(u_.vector())).max()
    max_u = np.abs(np.array(u_n.vector()) ).max()
    u_n.assign(u_)
    p_n.assign(p_)

    print('t = %.4f: error = %.6g, max_u = %.6g' % (t, error, max_u))

File("u_n") << u_n

I want to make my solvers run in parallel without splitting the mesh. How i can solve it?

TherealSova · June 1, 2021, 3:57pm

If i use:

mpirun -n 8 python3 3D_tube_script.py

I get this error:

Normal read!
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
  File "3D_tube_script.py", line 13, in <module>
    meshio.write("mesh.xdmf", tetra_mesh)
  File "/usr/local/lib/python3.8/dist-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 535, in write
    XdmfWriter(*args, **kwargs)
  File "/usr/local/lib/python3.8/dist-packages/meshio/xdmf/main.py", line 346, in __init__
    self.h5_file = h5py.File(self.h5_filename, "w")
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 406, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_hl/files.py", line 179, in make_fid
    fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 108, in h5py.h5f.create
OSError: Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')

dokken · June 1, 2021, 4:10pm

You should create and read your mesh in serial as neither Gmsh (by default) or meshio supports parallel execution. As reading in xdmf files is supported in parallel, i suggest you split your script into two pieces:

script creating and converting the mesh
one reading it into dolfin and solving the problem of interest.

You could Also spawn Gmsh or meshio on only one rank, as illustrated in:
https://jorgensd.github.io/dolfinx-tutorial/chapter3/subdomains.html#subdomains-defined-from-external-mesh-data

TherealSova · June 1, 2021, 4:17pm

Thank very mach! For script above, I must split:
One part:

import meshio
msh = meshio.read("gmsh/3D_channel.msh")
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

tetra_mesh = create_mesh(msh, "tetra", prune_z=True)
meshio.write("mesh.xdmf", tetra_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("facet_mesh.xdmf", triangle_mesh)

print("Normal read!")

Two part:

from dolfin import * 

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool((x[2] < DOLFIN_EPS or x[2] > (1.0 - DOLFIN_EPS)) \
                    and on_boundary)


    def map(self, x, y):
        y[0] = x[0] 
        y[1] = x[1]
        y[2] = x[2] - 1

pbc = PeriodicBoundary()  
V = VectorFunctionSpace(mesh, "P", 2, constrained_domain=pbc)
P = FunctionSpace(mesh, "P", 1)

dim = mesh.topology().dim()
print('Dim:',dim)

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

bcp_inflow = DirichletBC(P, Constant(8), facet_domains, 6)
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), facet_domains, 8)
bcp_outflow = DirichletBC(P, Constant(0), facet_domains, 7)

bcu = [bcu_walls]
bcp = [bcp_outflow, bcp_inflow]

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

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

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

# Define expressions used in variational forms
U  = 0.5*(u_n + 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) / 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_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

print("Go to assemble!")
# 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]

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')
    
    error = np.abs(np.array(u_n.vector()) - np.array(u_.vector())).max()
    max_u = np.abs(np.array(u_n.vector()) ).max()
    u_n.assign(u_)
    p_n.assign(p_)

    print('t = %.4f: error = %.6g, max_u = %.6g' % (t, error, max_u))

File("u_n") << u_n

Did I understand the example correctly?

dokken · June 1, 2021, 4:38pm

Yes, that is correct

Topic		Replies	Views
How can I solve transient problems in parallel?	3	1962	April 3, 2019
Problem with solve() in parallel computing	0	742	February 21, 2020
Help with parallel processing of a viscoelastic stokes flow system	3	718	June 16, 2021
Parallelizing heat equation tutorial	9	571	May 25, 2023
Parallel for-loop with multiple solve-calls, MPI dolfin-adjoint	13	3641	December 10, 2020

Parallel solve from solve method

Related topics