[hxaxtma]

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:

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

[hxaxtma]

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:

            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:

    // 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:

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:

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?

[hxaxtma]

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.

    // 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>());

		}

[alexeym]

Hi,

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

[hxaxtma]

Hi Alex,
thank you for your advice,
I was able to collect and combine the labels by using the following code snippet.

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.

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;
    
}

[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.

[hxaxtma]

No problem, in the following is the code. My source term is a vectorField, calculated by flowDir * localPressureGradient for each cell.

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;   
}

[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.

[hxaxtma]

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

[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.

[hxaxtma]

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

[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

[alexeym]

Hi,

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

    scalarField gradPlocal_(gradP0_+ dGradPlocal_);

instead of

    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).

[hxaxtma]

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]

Solved , index error:

Instead of

U[cellI] =  flowDir_ * rAU[cellI] * dGradPlocal_[cellI];

it should be

U[cellI] =  flowDir_ * rAU[cellI] * dGradPlocal_[i];