[GustavGans]

Hey there,

maybe you can help me with the following problem I'm currently dealing with:

In my solver, I defined a volScalarField border which only contains the values 0 and 1.
Now I would like to set the values of a volVectorField U to zero, but only in the cells where the value of border is 1.

So far, I'm doing this with a forAll-loop, as follows:

forAll(U,i)
{
    if(border[i] == 1)
    {
        U[i].component(0) = 0;
        U[i].component(1) = 0;
        U[i].component(2) = 0;
    }
}

but since this is an expensive operation, I would like to avoid this. Instead, I would like to use a Matrix-operation, something like (This is only Pseudo-Code, just to show my plan)

if(border == 1)
{
    U = 0;
}

Any ideas?

[metalfox]

You can use either

U = vector::zero

or

U = vector(0.0,0.0,0.0)

[GustavGans]

Thanks for your fast answer!

But using these two options will result in the dimension error, because U has the dimension m/s.

--> FOAM FATAL ERROR:
Different dimensions for =
    dimensions : [0 1 -1 0 0 0 0] = [0 0 0 0 0 0 0]

[metalfox]

Then try:

U = dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)

[GustavGans]

Thanks again, this works for setting the U-values to zero. But what about my if-condition? How do I have to write my if-condition, that it only sets the U-values to zero in these cells, where border = 1?

[metalfox]

I can't test right now, but I would say this should work:

U[i] = vector::zero

If not, try the dimensioned version.

[GustavGans]

I found another (quite simple) way to solve the problem. If anyone is interested: Just multiply the volVectorField U with the function (-1)*(border-1). This function equals 0 for all cells where border=1 and equals 1 for all cells where border=0.

U = (-1)*(border-1)*U;

So I neither need an if-condition nor a forAll-loop. Nevertheless, thanks a lot for your help metalfox!

[metalfox]

I'm glad you have found a nice solution!

Re-reading your question I realized that you stated that a forAll is not efficient. That's not true. This is C++ not matlab.

In general, if you can encompass all your calculations in a single forAll, that's much better that several "hiden" forAll, because of cache misses and redundant memory accesses.

Do not presume that every operation in OpenFOAM is evaluated lazily (i.e. with a single forAll).

Said that, I prefer code readability over performance.

[cuikahlil]

HI!
I'm doing something somewhat related.
I want to set the velocity of a cell to 0 if the value or alpha1 (water phase) is less than a certain value (like <0.01)

I implemented this in Ueqn right after the convective term:

forAll(U, celli)
	{
		if (alpha1[celli] < 0.01) 
		{
			U[celli] == vector(0.0,0.0,0.0);
			
		}
	}

but it doesn't re-assign the values.
Maybe I'm positioning the code wrong? or it the way it is written?
Any help would be very much appreciated. Thank you!

[GustavGans]

Hey Cuikahlil,

in which solver are you implementing this code? For me it seems right and I copied your code into my own pimpleFoam-based solver, here it works.

[cuikahlil]

Hi!

Thanks for your reply.
Actually I realized late that it did work. What I intially expected was that all U values will actually be zero, but I didn't realize that I was only zero-ing the convective contribution and not the contribution from the other side.

Anyways, sorry for the bother and thank you for your time.

[BuzzB]

Hey,


I've used a similar approach like GustavGans in my solver:

U *= (1-field);

field is containing 0 and 1 values, depending on the temperature.


However, I was wondering at which place it is the best to implement this code. I tested it at the end of the UEqn.H call. But the area, where field equals 1 is still seeing some velocity calculation and values that are greater than 0.


Where exactly did you guys implement this approach, in order to get it to work as intended? Thanks in advance!

[BuzzB]

I came up with a solution:


We, of course, have to implement this line after the momentum corrector in pEqn was called, so e.g. at the end of pEqn.C.


However, to me, this doesn't seam to be an elegant solution. So I will try to add a sink term to the PDE. I will post the solution here, once it is satisfactory.