Not Getting correct value using scatter_reverse in parallel

dolfinx
1

Hello,
I am trying to convert my code to parallel programming(was working perfectly on single processor). I am not getting corrected value of mass=Function(W). After calculating mass.x.array, I am using scatter_reverse(la.InsertMode.add) to add all the value at shared nodes but it seems they are not adding up correctly.

Terminal command: mpirun -np 2 python3 test.py
Fenicsx version: 0.9

from dolfinx import mesh, fem, la
import numpy as np
import ufl
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
import gmsh_api.gmsh as gmsh

def create_cube():
    # Gmsh initialization
    gmsh.initialize()
    gmsh.model.add("GradedCube")

    gdim = 3
    fdim = 2

    # mesh refinement
    a = 1.0     # for bottom refinement
    b = 1.0     # for top refinement

    # points
    gmsh.model.geo.addPoint(0, 0, 0, lc / a, 1)
    gmsh.model.geo.addPoint(1, 0, 0, lc / a, 2)
    gmsh.model.geo.addPoint(1, 1, 0, lc / b, 3)
    gmsh.model.geo.addPoint(0, 1, 0, lc / b, 4)
    gmsh.model.geo.addPoint(0, 0, 1, lc / a, 5)
    gmsh.model.geo.addPoint(1, 0, 1, lc / a, 6)
    gmsh.model.geo.addPoint(1, 1, 1, lc / b, 7)
    gmsh.model.geo.addPoint(0, 1, 1, lc / b, 8)

    # lines
    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(2, 3, 2)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.addLine(4, 1, 4)
    gmsh.model.geo.addLine(1, 5, 5)
    gmsh.model.geo.addLine(5, 6, 6)
    gmsh.model.geo.addLine(6, 7, 7)
    gmsh.model.geo.addLine(7, 8, 8)
    gmsh.model.geo.addLine(8, 5, 9)
    gmsh.model.geo.addLine(4, 8, 10)
    gmsh.model.geo.addLine(3, 7, 11)
    gmsh.model.geo.addLine(2, 6, 12)

    # curve loop
    gmsh.model.geo.addCurveLoop([1, 12, -6, -5], 1)
    gmsh.model.geo.addCurveLoop([12, 7, -11, -2], 2)
    gmsh.model.geo.addCurveLoop([11, 8, -10, -3], 3)
    gmsh.model.geo.addCurveLoop([10, 9, -5, -4], 4)
    gmsh.model.geo.addCurveLoop([8, 9, 6, 7], 5)
    gmsh.model.geo.addCurveLoop([3, 4, 1, 2], 6)

    # surface
    gmsh.model.geo.addPlaneSurface([1], 1)
    gmsh.model.geo.addPlaneSurface([2], 2)
    gmsh.model.geo.addPlaneSurface([3], 3)
    gmsh.model.geo.addPlaneSurface([4], 4)
    gmsh.model.geo.addPlaneSurface([5], 5)
    gmsh.model.geo.addPlaneSurface([6], 6)

    # surface loop
    gmsh.model.geo.addSurfaceLoop([5, 3, 2, 1, 6, 4], 1)

    # volume
    gmsh.model.geo.addVolume([1], 1)

    # synchronize
    gmsh.model.geo.synchronize()

    gmsh.model.addPhysicalGroup(gdim, [1], 1)
    gmsh.model.addPhysicalGroup(fdim, [1], 1)
    gmsh.model.addPhysicalGroup(fdim, [2], 2)
    gmsh.model.addPhysicalGroup(fdim, [3], 3)
    gmsh.model.addPhysicalGroup(fdim, [4], 4)
    gmsh.model.addPhysicalGroup(fdim, [5], 5)
    gmsh.model.addPhysicalGroup(fdim, [6], 6)

    # generate a 3D mesh
    gmsh.model.mesh.generate(3)

    # Convert to Dolfinx
    msh, cll_tgs, facet_tgs = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim) # , partitioner=mesh.GhostMode.shared_facet)

    gmsh.finalize()

    return msh

def main():

    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    size = comm.Get_size()

    # create mesh
    domain = create_cube()

    # Dimension of mesh/domain
    dim = domain.topology.dim

    # Ensure mesh topology is initialized
    domain.topology.create_connectivity(dim, dim - 1)
    domain.topology.create_connectivity(dim, 0)
    domain.topology.create_connectivity(0, dim)

    # Get connectivity matrix
    cells_to_vertices = domain.topology.connectivity(dim, 0)

    # Number of cells
    num_cells = domain.topology.index_map(dim)

    # definition of function space
    degree = 1
    W = fem.functionspace(domain, ("P", degree))
    DG = fem.functionspace(domain, ("DG", 0))

    mass = fem.Function(W)

    # Volume of each cell
    cell_volume_form = fem.form(ufl.TestFunction(DG) * ufl.dx)
    volumes = fem.assemble_vector(cell_volume_form)

    for cell_index in range(num_cells.size_local):
        cell_vertices = cells_to_vertices.links(cell_index)
        volume = volumes.array[cell_index]

        # Calculating mass vector
        mass_per_vertex = rho * volume / len(cell_vertices)
        mass.x.array[cell_vertices] += mass_per_vertex

    mass.x.scatter_reverse(la.InsertMode.add)
    mass.x.scatter_forward()

    if rank == 0:
        fl.write(f"Number of processors: {size}\n")

    fl.write(f"Rank {rank}: Mass vector: \n {mass.x.array[:]}\n")

    mass_form = fem.form(ufl.TestFunction(W) * rho * ufl.dx)
    mass_vector = fem.assemble_vector(mass_form)
    mass_vector.scatter_forward()
    fl.write(f"Rank {rank}: Mass vector: \n {mass_vector.array[:]} \n")


if __name__ == "__main__":
    # Material properties (All in SI units)
    rho = 2.0  # Density (g/mm3)
    lc = 0.5  # mesh size

    # Post-processing .txt file
    filepath = "results_new/"
    filename = "test"
    fl = open(filepath + filename + '.txt', 'a')

    main()

I would really appreciate some help to obtain an array that correctly represent the mass vector.

Thank you in advance !

2

I found out the problem!!
I was using connectivity matrix which is giving incorrect cell to node relation for shared node

Instead of above, I used geometry.dofmap which is providing correct mapping.

Developer, Please look into it.

Thank you!!