CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Set values of a volVectorField to zero (https://www.cfd-online.com/Forums/openfoam-programming-development/198660-set-values-volvectorfield-zero.html)

GustavGans February 14, 2018 06:26

Set values of a volVectorField to zero
 
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:

Code:

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)

Code:

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

Any ideas?

metalfox February 14, 2018 06:38

You can use either
Code:

U = vector::zero
or
Code:

U = vector(0.0,0.0,0.0)

GustavGans February 14, 2018 07:00

Thanks for your fast answer!

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

Code:

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


metalfox February 14, 2018 07:11

Then try:
Code:

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

GustavGans February 14, 2018 07:18

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 February 14, 2018 07:25

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

U[i] = vector::zero
If not, try the dimensioned version.

GustavGans February 14, 2018 08:11

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.

Code:

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 February 15, 2018 09:58

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 August 30, 2018 06:43

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

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 September 3, 2018 05:27

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 September 4, 2018 09:02

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 March 28, 2019 02:04

Hey,


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


Code:

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 April 9, 2019 00:12

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.


All times are GMT -4. The time now is 02:38.