CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Loop through processors and collect cellLabels of celLZone

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   March 14, 2017, 14:30
Default Loop through processors and collect cellLabels of celLZone
  #1
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Hi guys,
maybe you have an advice.
I want to loop through a a cellZone collect data on master node, do something and sent data back on a parallelized case. Any advice how to code this properly?

This is what I have so far, unfortunately the lower part seems not to work.
Code:
    // get the number of processors

	labelList fieldSizePerProc( Pstream::nProcs(), 0.0 );
	label thisProcNb = Pstream::myProcNo();
	fieldSizePerProc[thisProcNb] = thisPntFieldSize;

	reduce(fieldSizePerProc, sumOp<labelList>() );

	int totalSize = sum(fieldSizePerProc );
	labelList gobalCells_(totalSize,0.0);


	for( int i = 0; i < Pstream::nProcs(); i++)
		{
			if( thisPntFieldSize > 0 )
			{
				for( int j = 0; j < cells_.size(); j++)
				{
					gobalCells_[i][j]=cells_[j];
				}
			}
		reduce(gobalCells_, sumOp<labelList>());

		}
hxaxtma is offline   Reply With Quote

Old   March 15, 2017, 05:34
Default
  #2
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,590
Rep Power: 26
alexeym will become famous soon enoughalexeym will become famous soon enough
Hi,

There is gatherList function (https://cpp.openfoam.org/v4/a02062.h...6456ae1b323211) for collecting lists from processors.
alexeym is offline   Reply With Quote

Old   March 15, 2017, 16:25
Default
  #3
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Hi Alex,
thank you for your advice,
I was able to collect and combine the labels by using the following code snippet.

Code:
void Foam::fv::localVelocityForce::correct(volVectorField& U)
{
    const scalarField& rAU = rAPtr_().internalField(); //1.0/UEqn.A() Diagonal part of matrix U 
...
List<List<label>> gatheredcells_(Pstream::nProcs());
gatheredcells_[Pstream::myProcNo()] = cells_;
Pstream::gatherList(gatheredcells_);

List<label> cells_  = ListListOps::combine<List<label> >(gatheredcells_, accessOp<List<label> >());

.... // Calculate Source Term.
I was also able to calculate my source, which I want to impose the fluid field.

Since I am using the fvOption framework, the next problem seems to appear in the addSup function:

Unfortunately I am nut sure where to put the Pstream::gatherList function for the source, since the addsup function below is using UIndirectList assignment.


Code:
void Foam::fv::localVelocityForce::addSup
(
    const vectorField& ibSource_,
    fvMatrix<vector>& eqn,
    const label fieldI
)
{   
       
	DimensionedField<vector, volMesh> Su
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "_Su",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedVector("zero", eqn.dimensions()/dimVolume, vector::zero)
    );
 

	UIndirectList<vector>(Su,cells_) = source;

    eqn += Su;
    
}
hxaxtma is offline   Reply With Quote

Old   March 15, 2017, 16:46
Default
  #4
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,590
Rep Power: 26
alexeym will become famous soon enoughalexeym will become famous soon enough
Hi hxaxt,

Since it is YOUR source term, it is YOUR, who decides, where to put gather-scatter operation.

If you have DOUBTS, you can describe your source term in more details, so other people can propose their solutions without necessity to use mind reading powers.
alexeym is offline   Reply With Quote

Old   March 15, 2017, 16:49
Default
  #5
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
No problem, in the following is the code. My source term is a vectorField, calculated by flowDir * localPressureGradient for each cell.
Code:
void Foam::fv::localVelocityForce::correct(volVectorField& U)
{
    const scalarField& rAU = rAPtr_().internalField(); //1.0/UEqn.A() Diagonal part of matrix U [s/m]




List<List<label>> gatheredcells_(Pstream::nProcs());
gatheredcells_[Pstream::myProcNo()] = cells_;
Pstream::gatherList(gatheredcells_);

List<label> cells_  = ListListOps::combine<List<label> >(gatheredcells_, accessOp<List<label> >());

//for (int i=0; i<Pstream::nProcs() ;i++)
 //if (pntFieldSizePerProc[i]>0)
 //{
	//Info << "Processor Nr "<< i << "with cellSize :"<< pntFieldSizePerProc[i] << endl;
	// Integrate flow variables over cell set
	scalarField rAUlocal_(cells_.size(),0.0);//_(cells_.size(),0);
	scalarField magUbarlocal_(cells_.size(),0.0);

	forAll(cells_, i)
	{
		label cellI = cells_[i];
		rAUlocal_[i] = rAU[cellI];
		magUbarlocal_[i] = (flowDir_ & U[cellI]);
	}

	// Calculate the pressure gradient increment needed to adjust the locale
	// flow-rate to the desired value

	dGradPlocal_ = relaxation_*(mag(Ubar_) - magUbarlocal_)/(1e-12+rAUlocal_);

	// Apply correction to velocity field          
	forAll(cells_, i)
	{
		label cellI = cells_[i];
		U[cellI] =  flowDir_ * rAU[cellI] * dGradPlocal_[cellI];
	}


	scalarField gradPlocal_ = gradP0_+ dGradPlocal_; 

	vectorField ibSource_(cells_.size(),vector::zero);
	ibSource_ = (flowDir_*gradPlocal_);

	writePropslocal(gradPlocal_);



	if(debug==1)
	{
		Info << "magUbarlocal" << magUbarlocal_ << endl;
		Info << "rAUlocal_" << rAUlocal_ << endl;
		Info << "gradPlocal_" << gradPlocal_ << endl;
		Info << "ibSource_" << ibSource_ << endl;	
		Info<< "ibSource Size:  = " << ibSource_.size() << endl;
	}
}


void Foam::fv::localVelocityForce::addSup
(
    const vectorField& ibSource_,
    fvMatrix<vector>& eqn,
    const label fieldI
)
{   
       
	DimensionedField<vector, volMesh> Su
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "_Su",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedVector("zero", eqn.dimensions()/dimVolume, vector::zero)
    );
	UIndirectList<vector>(Su,cells_) = ibSource_;

    eqn += Su;   
}
hxaxtma is offline   Reply With Quote

Old   March 15, 2017, 17:04
Default
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,590
Rep Power: 26
alexeym will become famous soon enoughalexeym will become famous soon enough
Mmmm... did I ask for the code?

Right now it looks rather strange, you collect cell labels from all processors and then use these cells in a separate process (so maybe current processor does not have cell with given ID).

OK. Now we do not need mind reading powers to guess your fvOptions aim, now we have to decipher your code. Nice. Though, maybe it is just me, maybe to others everything is already clear.
alexeym is offline   Reply With Quote

Old   March 15, 2017, 17:56
Default
  #7
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Hi,
The problem is I am not sure how to handel the addSup void function within the fvOptions framework. Think this functon is a need. So my goal was to collect the data in the correct() function, and return the source term to the addSup function, keeping it like it was

Maybe this is totally wrong here. Do you have any suggestion?
Thanks a lot for your time
hxaxtma is offline   Reply With Quote

Old   March 19, 2017, 08:10
Default
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,590
Rep Power: 26
alexeym will become famous soon enoughalexeym will become famous soon enough
Well, it is still not clear what your are trying to achieve. Your code is just strange modification of meanVelocityForce fvOption. For some strange and unknown reason, you collect cell IDs from all processors, then use cell IDs from other processors for address your local velocity field and do other bizarre things.
alexeym is offline   Reply With Quote

Old   March 19, 2017, 09:40
Default
  #9
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Hi Alexeym,
well my goal is the following:

Instead of only calculating meanvalues over a cellSet, as it is in meanVelocityForce given, I want to calculate local values in order to get the influence of a local velocity profile in front of the cellSet. Imagine a wall mounted cylinder in a blasius boundary layer.

Here the main goal is to calculate the pressure gradient up to zero velocity is achieved. Let's say we call it mimicking a immersed boundary by explicit pressure gradient.

The code here so far works very well for serial runs and does exactly what I wanted to do.
For parallel runs it results in strange behaviour. For some reason I cannot find are the cells of the cellZone not collected in parallel run.
That's why I tried to combine all labels to one master node, calcualte data and spread it again.

After the last days and due to your comments I am sure that this is no necessary. Cause forAll loops are acting independently on every processor.

Therefore I also started debugging the code by using OF in debug mode.

In parallel (for a minimal example with overall 480 cells in the cellSet, divided by 4 processors in x-direction results in (0 240 240 0) cells per processor), I got always this error

Code:
[0] 
[0] 
[0] --> FOAM FATAL ERROR: 
[0]  Field<scalar> f1(240), Field<scalar> f2(240) and Field<scalar> f3(0)
    for operation f1 = f2 + f3
[0] 
[0]     From function void Foam::checkFields(const Foam::UList<T>&, const Foam::UList<Key>&, const Foam::UList<Type3>&, const char*) [with Type1 = double; Type2 = double; Type3 = double]
[0]     in file /opt/openfoam4_dbg/src/OpenFOAM/lnInclude/FieldM.H at line 75.
and this cannot run for sure. I a just not able to figure out which line really causes the problem.

If you have time and you are willing to help, I uploaded the library + testcase on GitHub under: https://github.com/gaberchen/localVelocityForce

And last but not least, the overall goal is to combine this method within an actuator line method in order to model the tower effect (bluff body) more precisely.

Thanks for your time
hxaxtma is offline   Reply With Quote

Old   March 19, 2017, 15:32
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,590
Rep Power: 26
alexeym will become famous soon enoughalexeym will become famous soon enough
Hi,

Unfortunately I can not reproduce error from your last post (repository commit fd812810a5689dfa05160270ea6055bb3a409c46). I have made only one modification (due to clang):

Code:
    scalarField gradPlocal_(gradP0_+ dGradPlocal_);
instead of

Code:
    scalarField gradPlocal_ = gradP0_+ dGradPlocal_;
And, when I try to run test case (blockMesh -> topoSet -> decomposePar -> mpiexec -np 4 simpleFoam -parallel), I got no errors (usually during second run, first run usually has certain issues, but they are not related to fvOption).
alexeym is offline   Reply With Quote

Old   March 20, 2017, 07:51
Default
  #11
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Hi Alex,
yes, thanks to Pete we found the bug on Sunday, so you looked at the actual version which runs in parallel.

Unfortunately, parallel run only works for certain combinations of decomposition. For higher cpu counts it fails again. Need to figure this out why this happens.

Next, as you already mentioned, sometimes it only runs on second try. Would be also interesting to have a closer look where this comes from. Any idea from your experience?


I appreciate your help,
hxaxtma is offline   Reply With Quote

Old   March 22, 2017, 15:08
Default
  #12
Senior Member
 
Join Date: Jan 2014
Posts: 127
Rep Power: 5
hxaxtma is on a distinguished road
Solved , index error:

Instead of
Code:
U[cellI] =  flowDir_ * rAU[cellI] * dGradPlocal_[cellI];

it should be

Code:
U[cellI] =  flowDir_ * rAU[cellI] * dGradPlocal_[i];
hxaxtma is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 21:05.