CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User 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 7, 2017, 07:10
Default Loop through processors and collect cellLabels of celLZone
  #1
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
hxaxtma is on a distinguished road
Hi guys,
did someone of you tried to rewrite the meanVelocityForce fvOptions files to a localVelocityForce one?
Which means I want to model a non-slip wall by imposing velocity to (0 0 0) with respect to local velocity distribution at each cell of the cellset.

So far I managed to rewrite the code successfully up to here:
Code:
void Foam::fv::localVelocityForce::addSup
(
    fvMatrix<vector>& eqn,
    const label fieldI
)
{   
       
	scalar gradP = gradP0_ + dGradP_;
	scalarField gradPlocal_ = gradP0_ + dGradPlocal_;

	DimensionedField<vector, volMesh> SuI
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "_SuI",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedVector("zero", eqn.dimensions()/dimVolume, vector::zero)
    );
/*
    forAll(cells_, i)
    {
		label cellI = cells_[i];
		SuI[cellI] += flowDir_*gradPlocal_[cellI];
	}
*/


		UIndirectList<vector>(SuI,cells_) = flowDir_*gradPlocal_; ->This line is the problem!!!!

    eqn += SuI;
  
}
The problem is caused by gradPlocal_ which is a scalarField multiplied by flowDir_ equals in a volVectorField.

So my question is, how can I use UIndirectList as a volVectorField.
Or is there a way to impose Sui as volVectorField directly within the last line eqn +=SuI?
Thanks in advance

Last edited by hxaxtma; March 8, 2017 at 19:39.
hxaxtma is offline   Reply With Quote

Old   March 10, 2017, 09:37
Default Problems parallel run with cellZones while serial is running
  #2
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
hxaxtma is on a distinguished road
Hi guys,
following problem occurs (using openfoam 4.1) in my case. Calculating some data over a cellZone or cellSet runs in serial absolutely fine, in parallel the solver crashes.

In serial magUbarlocal is filled with correct values.
In parallel no values are inserted -> Solver crashes

Here is a small code snippet causing the problem (implementing within the fvOptions framework). where cells_ contain the label list of the cellZone:
Code how cell labels is obtained:
Code:
            label zoneID = mesh_.cellZones().findZoneID(cellSetName_);
            if (zoneID == -1)
            {
                FatalErrorInFunction
                    << "Cannot find cellZone " << cellSetName_ << endl
                    << "Valid cellZones are " << mesh_.cellZones().names()
                    << exit(FatalError);
            }
            cells_ = mesh_.cellZones()[zoneID];

Code which causes the problem:

Code:
    // Integrate flow variables over cell set
    ...
    const volVectorField& U;
    ...
    scalarField magUbarlocal_(cells_,size(),0.0);


    const scalarField& cv = mesh_.V();
    forAll(cells_, i)
    {
        label cellI = cells_[i];
        magUbarlocal_[i] = (flowDir_ & U[cellI]);

    }

    // Collect across all processors
    reduce(magUbarlocal_, sumOp<scalarField>());
Error is here:
Quote:
void Foam::Pstream::gather<Foam::Field<double>, Foam::sumOp<Foam::Field<double> > >(Foam::List<Foam::UPstream::commsStruct> const&, Foam::Field<double>&, Foam::sumOp<Foam::Field<double> > const&, int, int) at ??:?
Running OPenFOAm in Debug Mode gives following hint:
Quote:
Time = 1

smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 0.0842275
526312, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.971361820323, Final residua
l = 0.0643989023786, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 0.0662730
562513, No Iterations 2
[0]
[0]
[0] --> FOAM FATAL ERROR:
[0] Field<scalar> f1(0), Field<scalar> f2(0) and Field<scalar> f3(240)
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; Typ
e2 = double; Type3 = double]
[0] in file /opt/openfoam4_dbg/src/OpenFOAM/lnInclude/FieldM.H at line 75.
[0]
FOAM parallel run aborting
....
:
Next, in FieldM.H , line 75 the command if (f1.size() != f2.size() || f1.size() != f3.size()) results in the above error.
So again, the error is caused by non inserting the local values into the array.

Any advice what I am doing wrong here?

Last edited by hxaxtma; March 10, 2017 at 11:09.
hxaxtma is offline   Reply With Quote

Old   March 14, 2017, 14:30
Default Loop through processors and collect cellLabels of celLZone
  #3
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #4
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,816
Rep Power: 31
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
  #5
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,816
Rep Power: 31
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
  #7
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,816
Rep Power: 31
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
  #9
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,816
Rep Power: 31
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
  #11
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #12
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,816
Rep Power: 31
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
  #13
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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
  #14
Senior Member
 
Join Date: Jan 2014
Posts: 140
Rep Power: 6
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 23:13.