CFD Online Discussion Forums - Findcell in parallel

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Findcell in parallel (https://www.cfd-online.com/Forums/openfoam-programming-development/60958-findcell-parallel.html)

Findcell in parallel

Hi,

I had some problems defining a momentum source. the first step was to locate the cells that should have nonzero source. If I make the call

label point=mesh.findCell(r) // r is vector located in the domain

if ( point >= 0)
{
Info << "cell found";
}

It works fine in serial, but when I decompose the domain and run in parallel the cell cannot be found. Did I miss something here? Do I need to explicit loop over the processors or what?

I guess this is trivial for most of you so remember; a short answer is infinately more appreciated than no answer http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

Thanks!

Every processor looks for that

Every processor looks for that vector in it's own sub-Grid. Only the one that has the right cell enters the if-statemen. As I understand it (havn't looked at THAT source) the Info-stream is only output for the first processor, all others are mute. So if the point is NOT on the first processor nothing is output (allthough he was found).

My fix (but I'm not very good with that) would be:

label point=mesh.findCell(r)
{
label tmpPoint=myPoint;

reduce(tmpPoint,maxOp<label>());
// is point on ANY processor bigger than 0
if ( tmpPoint >= 0)
{
Info << "cell found";
}
}
if(point>=0) {
// do it only on the processor that HAS this point
}

Bernhard is right, if you want

Bernhard is right, if you want all processors to be able to report the message, use

Pout << "cell found"

Instead of Info. This will prefix the statement with the originating processor number as well (say thanks to Mattijs).

Thanks alot guys! That rea

Thanks alot guys!

That really made me alot clearer on whats going on. Now it works like a charm (Although it should have worked before also but the output was not what i believed it to be).

/M

Hi Foamers,
it's possible to use findCell to find a cell on a processor that's different from the local one? If not, is there any workaround to do it?
Thanks, Ivan

If I understand your question, a possible solution is to use findCell on all the processors, then reduce the result so that each processor knows where the requested point is.

Hope this helps,
Francesco

Quote:

Originally Posted by fra76 (Post 236222)

Hi Francesco,
unfortunately your suggestion can't work, as in my code each cell have to search points that can be on different processors, and findCell finds cells that belong to the local processor...

If it's only one point (I'm writing by heart, so please forgive errors):

Code:

labelList localCellId(Pstream::nProcs(),0);

localCellId[Pstream::myProcNo()] = mesh.findCell(YOUR_POINT);

reduce(localCellId, maxOp<labeList>());

After that, the list localCellId should be != -1 only in the entry that corresponds to the cell owner, with the value that is the owner's cell id.

If you have to look for a point in the mesh, I think that you have to spread this information among all the processors before searching. Each processor can only look for things in its domain, and then can communicate the result of the search operation to the others, if needed.

Francesco

Hi,
I have a somewhat related problem.
I want to perform a simple post-processing task (in parallel, using meshsearch).
Basically:

Quote:

- each processor loops over its interior field

-> current_pt=mesh.C()[cellI];
-> U_current[cellI]=U[cellI];

- if the value of field X (here mesh.C().component(2)) is above a threshold, it looks for an offset point:

-> sup_pt=current_pt+point(0.0 , 0.0, offset);
-> idx_sup = searchEngine.findCell(sup_pt)
-> U_sup[cellI]= Uint.interpolate(sup_pt,idx_sup)

- the new field Y gets the value Y[cellI]=U_curr[cellI]-U_sup[cellI]

Here is my implementation:

Quote:

meshSearch searchEngine(mesh);
interpolationCellPoint<vector> UInt(U);

vector U_sup=vector::zero;
vector U_curr=vector::zero;
point current_pt=point::zero;
point pt_sup=point::zero;

scalar offset=30.0;

forAll(mesh.C(), cellI)
{

current_pt=mesh.C()[cellI];
U_curr=U[cellI];

if ( current_pt[2]>offset )
{
pt_sup=current_pt+point(0.0,0.0,offset);

Pout << " pt :x "<< current_pt << endl;
Pout << " pt sup: "<< pt_sup << endl;

List<label> idx_sup(Pstream::nProcs(),0);
List<scalar> dist_sup(Pstream::nProcs(),GREAT);
List<vector> Usup(Pstream::nProcs(),vector::zero);

idx_sup[Pstream::myProcNo()] = searchEngine.findNearestCell(pt_sup,cellI);
reduce(idx_sup,maxOp<List<label> >());

searcherSup=1000000.0;

for (int i=0; i<Pstream::nProcs() ;i++)
{
dist_sup[i]=mag(mesh.C()[idx_sup[i]]-pt_sup);
Pout << " pt sup found: "<< i << " "<< mesh.C()[idx_sup[i]] << " " << dist_sup[i] <<endl;
if (dist_sup[i]<searcherSup)
{
searcherSup=dist_sup[i];
found_sup=i;

}

}
Pout << " pt sup kept: "<< mesh.C()[idx_sup[found_sup]] << endl << endl;

if (Pstream::myProcNo()==found_sup)
{
Usup[Pstream::myProcNo()]=UInt.interpolate(pt_sup,idx_sup[found_sup]);

}

reduce(Usup,maxOp<List<vector> >());

U_sup = Usup[found_sup];
Y[cellI] = U_curr-U_sup;

Info << " U : "<< U_curr << endl;
Info << " U sup: "<< U_sup << endl;
}

}
Pout << "Finished! << endl;

It seems that the meshSearch finds the right points, but the problem is that processor2 (due to my decomposition and the offset value) finishes its task first, and the other processors seem to get stuck has soon as this happens.

(By the way, it is a simplified version of my problem. I am not trying to find offset points outside my mesh, as this sample code would suggest)

Do you have an idea why? I have been struggling for quiet a while!!!

Thanks,

Fabien

findCell in parallel and reduce

Fabien & Foamers,

Did you solve this problem?

I suspect it may be related to what I am seeing when I loop over cells in a decomposed mesh. Unfortunately, I don't have an answer only another question...

Code:

forAll(alpha.internalField(),cellI)

{

              vector centerR = ptsR[cellI];

               labelList targetCellId(Pstream::nProcs(),0);



               targetCellId[Pstream::myProcNo()] = mesh.findCell(centerR);



        reduce(targetCellId, maxOp<List<label> >());

        Pout << "targetCellId: " << targetCellId << endl;

}

Without reduce it returns:

Code:

[3] targetCellId: 4(0 0 0 -1)

[1] targetCellId: 4(0 -1 0 0)

[2] targetCellId: 4(0 0 -1 0)

[1] targetCellId: 4(0 288 0 0)

[0] targetCellId: 4{0}

[0] targetCellId: 4(1 0 0 0)

[3] targetCellId: 4(0 0 0 -1)

[2] targetCellId: 4(0 0 -1 0)

[2] targetCellId: 4(0 0 -1 0)

[0] targetCellId: 4(34 0 0 0)

[1] targetCellId: 4(0 289 0 0)

[3] targetCellId: 4(0 0 0 -1)

[0] targetCellId: 4(35 0 0 0)

[1] targetCellId: 4(0 290 0 0)

With reduce it returns:

Code:

[0] targetCellId: 4{0}

[1] targetCellId: 4{0}

[2] targetCellId: 4{0}

[3] targetCellId: 4{0}

[0] targetCellId: 4(1 288 0 0)

[1] targetCellId: 4(1 288 0 0)

[2] targetCellId: 4(1 288 0 0)

[3] targetCellId: 4(1 288 0 0)

[3] targetCellId: 4(34 289 0 0)

[2] targetCellId: 4(34 289 0 0)

[0] targetCellId: 4(34 289 0 0)

[1] targetCellId: 4(34 289 0 0)

[3] targetCellId: 4(35 290 0 0)

[2] targetCellId: 4(35 290 0 0)

[0] targetCellId: 4(35 290 0 0)

[1] targetCellId: 4(35 290 0 0)

I am expecting reduce to return one cell ID and the rest zeros. However, you can see that reduce seems to be combining the results of multiple finds from different CPU findCell requests! Each point only exists on one of the four decomposed mesh domains. It seems the CPUs can race ahead of each other and reduce then gives a muddled answer.

How can I match each CPU's unique findCell request with a unique found cell ID?

thanks and regards,

James

Has anybody been able to use findCell() in parallel or got a viable work-around? I'm really pulling my hair out over this one.

Thanks

Reply to Thomas Flint on findCells() in parallel

Hi Thomas,

findCells() tries to find the cell index of the given point in each decomposed mesh. If it doesn't find the point, it returns a cell index of -1. You can find what the function does in src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C. If you are looking to get the cell index of a single point, you can write an if statement so that your code only runs if the cell index is not -1. For example:

Code:

vector position(0.5,0.5,0.5);

label cellI = mesh.findCell(position);

if (cellI != -1) //if the point is in the sub-mesh 

{

   //do stuff

}

Hope this helps!

Quote:

Originally Posted by tom_flint2012 (Post 646643)

Has anybody been able to use findCell() in parallel or got a viable work-around? I'm really pulling my hair out over this one.

Thanks

Lets suppose that you want to store the velocity and density of a point of interest (myPoint). But at the same time you want your code to be able to handle parallel mode.

A possible solution would be:

Code:

  



    point myPoint=point(myX,myY,myZ);//initialize the point of interest

    label myCellId=mesh.findCell(myPoint);//search myPoint in the mesh and return the Cell id(if myCellId!=-1 then the cell is found) 

    



    vector myU = vector(VGREAT, VGREAT, VGREAT);//initialize velocity vector with a very high value (this will help later on when you use reduce command)

    scalar myRho = VGREAT;

    if (myCellId != -1)  //if cell is found 

    {

      myU =  U[myCellId];//store values in myU



      myRho = rho[myCellId];

    }

     reduce(myRho, minOp<scalar>());

     reduce(myU, minOp<vector>());//reduce vector

Hope that this example will help any future readers that may come in front of a similar problem.

Best,

Konstantinos