Parallelizing the filter code in OpenFOAM: Topology Optimization

[Ananth-Narayan-IITM]

Hello all,

I am writing a filter function for Topology optimization. I am trying to verify the code by generating a impulse (one cell in entire domain will have higher value (here its 1.0) rest of the cells have 0.0). Ideally this should produce diffused result as this is weighted average filtering.

I do have a problem as when I run in parallel, each local domain creates a impulse thus leading to 8 different impulse (as I used 8 processors). But when I run in serial, this code works good.

How do I modify the code such that I expect the same result in Parallel as in serial. I am trying to lookup things such as Pstream, syncTools etc., but required some support or clue.

Attaching results for serial and parallel
Serial Output:

Parallel Output:

dictionary topoDict = daOptionPtr_->getAllOptions().subDict("topologyOptimization");
scalar rMin = readScalar(topoDict.lookup("filterRadius"));
scalar beta = readScalar(topoDict.lookup("regularizationCurvature"));

// Debug start
label centerCell = mesh.nCells()/2; // Local for each processors
rhoN[centerCell] = 1.0;
// Debug end

forAll(muE, cellI)
{
    scalar numerator = 0.0;
    scalar denominator = 0.0;

    const vector& xc = mesh.C()[cellI];

    forAll(rhoN, j)
    {
        scalar d = mag(mesh.C()[j] - xc);

        if (d <= rMin)
        {
            scalar w = (rMin - d)/rMin;

            numerator += rhoN[j]*w;
            denominator += w;
        }
    }
    muE[cellI] = numerator/denominator;
    // Debug only for the impulse cell
    if (printToScreen_)
    {
        if (cellI == mesh.nCells()/2)
        {
            Info << "----------------------------------" << endl;
            Info << "Cell = " << cellI << endl;
            Info << "rMin = " << rMin << endl;

            label nSelected = 0;

            forAll(rhoN, j)
            {
                scalar d = mag(mesh.C()[j] - xc);

                if (d <= rMin)
                {
                    scalar w = (rMin - d)/rMin;

                    Info << "Cell = " << j
                        << " distance = " << d
                        << " weight = " << w
                        << " rhoN = " << rhoN[j]
                        << endl;

                    nSelected++;
                }
            }
            Info << "Number of selected cells = "
                << nSelected << endl;
        }
    }
    
}

if (printToScreen_)
{
    Info<< "Proc: " << Pstream::myProcNo()
    << "  Local cells: " << rhoN.size()
    << "  Mesh cells: " << mesh.nCells()
    << endl;
    Info<< "Proc patches:" << endl;

    forAll(mesh.boundaryMesh(), patchI)
    {
        Info<< patchI << " "
            << mesh.boundaryMesh()[patchI].type()
            << endl;
    }
}
alpha = 1.0 - (Foam::exp(-beta*muE)) + (muE*Foam::exp(-beta));