[hxaxtma]

Hi guys,

I am trying to translate a cellSet in parallel. Serial run is working properly, in parallel the cellSet is not spread over procBoundaries. cells_corresponds to a labellist of a cellSet.

Here is my code so far, any advice?

Thanks for your time.

                // get the number of processors
                List<List<label>> gatheredcells_(Pstream::nProcs());
                gatheredcells_[Pstream::myProcNo()] = cells_;
                Pstream::gatherList(gatheredcells_);
                Pstream::scatterList(gatheredcells_);

                List<label> globalCells_  = ListListOps::combine<List<label> >(gatheredcells_, accessOp<List<label> >()); // Combine all labels to one

                       forAll(cells_, i)
                        {


                            label cellI = cells_[i];                 
                            scalar dxt = C[cellI].component(0)  -   amplitude_.component(0)*(   std::sin((t)*frequency_.component(0)*2*constant::mathematical::pi)   -   std::sin((t+dt)*frequency_.component(0)*2*constant::mathematical::pi)  );
                            scalar dyt = C[cellI].component(1)  -   amplitude_.component(1)*(   std::sin((t)*frequency_.component(1)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(1)*2*constant::mathematical::pi)  );
                            scalar dzt = C[cellI].component(2)  -   amplitude_.component(2)*(   std::sin((t)*frequency_.component(2)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(2)*2*constant::mathematical::pi)  );
                            vector displacement(dxt,dyt,dzt);
                            label myCellId=mesh_.findNearestCell(displacement); //(meshSearch faster but does not wirk in parallel?)
                            if (myCellId != -1)  //if cell is found 
                            {
                                cells_[i]=myCellId;
                            }

                                      
                        }

[olesen]

Look at globalMeshData. You have gathered local cell ids, you'll need to use the global value.

[hxaxtma]

Originally Posted by olesen

Look at globalMeshData. You have gathered local cell ids, you'll need to use the global value.

Thanks for your advice.

I tried to figure it out how to use the globalMeshData and this is what I get so far.

1.) I was able to collect the global CellLabels
2.) displacement calculation is running
3.) New cellLAbels are obtained by findNearestCell
4.)How to scatter globalLabels back -> Open TODO?
5.) How to fill cells_ since I am using fvOPtions Open TODO?

Maybe you can point me again in the right direction.

 word cellSetName = ibSource_;
    cellSet myCellSet(mesh_, cellSetName);
    globalIndex globalNumbering(mesh_.nCells());
    labelList globalIds(myCellSet.sortedToc());
    for (label& id : globalIds)
    {
         id = globalNumbering.toGlobal(id);
    }
    Pout<< "Global cellIds " << globalIds.size() << nl;



    meshSearch ms(mesh_); 
    vectorField displacement;
    forAll(globalIds,i)
    {
        label cellI = globalIds[i];                 
        scalar dxt = C[cellI].component(0)  -   amplitude_.component(0)*(   std::sin((t)*frequency_.component(0)*2*constant::mathematical::pi)   -   std::sin((t+dt)*frequency_.component(0)*2*constant::mathematical::pi)  );
        scalar dyt = C[cellI].component(1)  -   amplitude_.component(1)*(   std::sin((t)*frequency_.component(1)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(1)*2*constant::mathematical::pi)  );
        scalar dzt = C[cellI].component(2)  -   amplitude_.component(2)*(   std::sin((t)*frequency_.component(2)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(2)*2*constant::mathematical::pi)  );
        vector tmpDisplacement(dxt,dyt,dzt);  
        displacement.append(tmpDisplacement);
        
    }



    //Calculate new global IDs due to tranlsation
    labelList globalIdsNew;
    forAll(displacement,i)
    {
        globalIdsNew.append(ms.findNearestCell(displacement[i],0,true)); //0 runs faster, -1 slower
    }

    //NOT SURE ABOUT THIS PART; HOW TO SCATTER BACK AND FILL cells_ (fvOPtions!) : cells_ is needed.
    for (label& id : globalIdsNew)
    {
         cells_[id]= globalNumbering.toLocal(id);
    }
    //NOT SURE ABOUT THIS PART; HOW TO SCATTER BACK AND FILL cells_ (fvOPtions!)
    
    volScalarField mask_
    (
        IOobject
        (
            "mask" + ibSource_,
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("mask" + ibSource_,dimensionSet(0,0,0,0,0,0,0),0) 
    );	
    

    forAll(cells_,i)
    {
        mask_[cells_[i]]=1.0;
    }

    if(runTime_.outputTime())
    {
    mask_.write();
    }

[olesen]

The whole point of global data is to translate between the processor-local numbering (eg, each processor has cells starting at zero) and a globally consistent numbering.
To do this you must apply this in the correct place.

// The cells to select - processor-local ids
labelList cellsToSelect(cellSet(mesh_, cellSetName).sortedToc()); 



... do something with them



.... more


// Want a global view. Eg, collecting them on master processor.
// Option1: renumber locally before sending

// Option2: send and renumber later



... option 1: this is better since each processor keeps busy



    const globalIndex globalNumbering(mesh_.nCells());



    const label localCellOffset = globalNumbering.localStart();
     for (label& celli : cellsToSelect)
    {
         celli += localCellOffset;
    }


// Finally collect back into master etc.


List<labelList> allCellIds(Pstream::nProcs);
allCellIds[Pstream::myProcNo()] = cellsToOffset;
...

[hxaxtma]

Hi,

thanks for your help, I tried the last days to find a solution. But still struggling with this case.

In general the concept is now more clear, but as you mentioned, the correct place of placement is important. I tried several ways but failed always. Maybe you can give me an advice here.

    labelList cells_(cellSet(mesh_, ibSource_).sortedToc());
    const globalIndex globalNumbering(mesh_.nCells());
    const label nTotalCells = globalNumbering.size();
    
    meshSearch ms(mesh_); 
    vectorField displacement;
    
    forAll(cells_,i)
    {
        label cellI = cells_[i];                 
        scalar dxt = C[cellI].component(0)  -   amplitude_.component(0)*(   std::sin((t)*frequency_.component(0)*2*constant::mathematical::pi)   -   std::sin((t+dt)*frequency_.component(0)*2*constant::mathematical::pi)  );
        scalar dyt = C[cellI].component(1)  -   amplitude_.component(1)*(   std::sin((t)*frequency_.component(1)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(1)*2*constant::mathematical::pi)  );
        scalar dzt = C[cellI].component(2)  -   amplitude_.component(2)*(   std::sin((t)*frequency_.component(2)*2*constant::mathematical::pi)  -   std::sin((t+dt)*frequency_.component(2)*2*constant::mathematical::pi)  );
        vector tmpDisplacement(dxt,dyt,dzt); 
        displacement.append(tmpDisplacement);  
        
    }


//Calculate newCellLabels by findNearestcell
    labelList cellsTrans_(displacement.size());
    forAll(displacement,i)
    {
        cellsTrans_[i]=ms.findNearestCell(displacement[i],0,true); 
    }

    /*//#############    Same as to.Global so can be missed
    DynamicList<label> localCellOffsetList;
    const label localCellOffset = globalNumbering.localStart();
     for (label& celli : cellsTrans_)
    {
         celli += localCellOffset;
         localCellOffsetList.append(celli);
    }
    */

   //Gather and scatter data
    forAll(cellsTrans_,i)
    {
    label globalCId = globalNumbering.toGlobal(cellsTrans_[i]); //to Master
    label proci = globalNumbering.whichProcID(globalCId); //and Back!
    label cellI=globalNumbering.toLocal(proci, globalCId);
    cells_[i]=cell]; //->
    }

[olesen]

I might have finally figured out what you are trying to do (I think). Since you've ask about parallel, will assume that you need to check things in parallel.
I think this means you need the following

• your starting cells (from cellSet), these are local ids and can stay that way. In fact you don't even care about them after the next step anyhow.
• You can create an nProcs-sized list of vectorFields. The size of each local vectorField is given by the number of cells in each processor-local cellSet. You fill these local fields with the correspond cell centre value - C[celli]
• While we are at it, we also apply your displacement function to obtain the new position (not sure why you stored displacements in your original code)
• now gather/scatter those lists. This ensure that each processor has the same knowledge of the new locations.
• perform your nearest cell search for these new positions. Can transform the nearest cell to either use global numbering, or simply tag your value with an additional processor id. At this stage it would also make sense to add in the absolute distance from this nearest cell and the stored displacement.
• Look for code using index point hit type of storage. You want to reduce these hits (containing cellid, procid, distance) according to the min distance. The cell with the min distance gets kept.
• The final stage is to "unwalk" these vectors and pick off the reduced closest hits that are on the local processor.

[hxaxtma]

Dear Mark,

i tried to follow your instructions. Still I am not able to scatter the right cellLabels back.
For clearification what I am trying to achieve, you can find an image in the attachement.

These are my steps:
1.) Create n-sized proc PointField of mesh_.C() and gather(scatter
2.) Calculate displacement using localIDs of cells_ = ccT_ (translated Coords)
3.) Gather Scatter ccT_ so that every proc has this information
4.) Perform meshSearch over ccT_ go geht new cellLabels : cCelltIDs_
5-) Calculate absDisplacement (at the moment nice to have, not using further at this stage)
6.) globalnumerbing of cCelltIDs_ to local

7.) Filling cells_ with new cCelltIDs_ (Fails in parallel)

The index hit point function is from this point of view not necessary, since from my understanding the findNearestCell funtcion returns already the nearest label, so in a first step this procedure should run as mentioned.

Thanks a lot for your help and time,

       
        //Create n-sized Proc PointField 
        List<Field<point> > allCo(Pstream::nProcs());
        allCo[Pstream::myProcNo()] = mesh_.C();
        Pstream::gatherList(allCo);
        Pstream::scatterList(allCo);

        //Calculate new Cooridnates by using local IDs cells_
        pointField ccT_(cells_.size(),vector::zero);
        forAll(cells_,i)
        {
            label cellI = cells_[i];                 
            ---do sin(2pif)
            ccT_[i]=point(dxt,dyt,dzt);
        }


        //Scatter all these Coordinates so that every proc knows 
        List<Field<point> > translatedCo(Pstream::nProcs());
        translatedCo[Pstream::myProcNo()] = ccT_;
        Pstream::gatherList(translatedCo);
        Pstream::scatterList(translatedCo);
        // //translatedCo = ListListOps::combine<Field<point> >(translatedCo, accessOp<Field<point> >());


meshSearch ms(mesh_);
        globalIndex globalCells(mesh_.nCells());

        (void)mesh_.tetBasePtIs();         // Force calculation of tet-diag decomposition (for use in findCell)
        labelList cellLabels(ccT_.size());
    
        forAll(ccT_,i)
        {
            const point& location=ccT_[i];
            label localCellI = ms.findNearestCell(location,-1,true);
            label globalCellI = -1;
            if (localCellI != -1)
            {
                globalCellI = globalCells.toGlobal(localCellI);
            }
            //reduce(globalCellI, maxOp<label>());

           label procI = globalCells.whichProcID(globalCellI);
           label procCellI = globalCells.toLocal(procI, globalCellI);
   
           Pout<< "Found point " << ccT_[i] << " in cell " << procCellI << " on processor " << procI << endl;
   
        if (globalCells.isLocal(globalCellI))
           {
               cellLabels[i] = localCellI;
           }
           else
           {
               cellLabels[i] = -1;
           }
       }

       cells_=cellLabels;