[IvyLong]

Hello,

I'm modifying a solver and one of my fields of interest is calculated by dividing one scalar by another, say A/B, with the new field C used in later calculations. However, I found that if, in my case folder, I fix both A and B at the boundaries with non-zero values, the calculated result of C is always zero at that boundary, which is obviously wrong.

I'm hoping there's a way to declare the C field differently, or possibly alter the way the b.c. are set in the case folder. Any suggestions?

Thanks

[boger]

I don't observe this. I read in two volScalarFields (A and B), each with fixedValue boundary conditions on all boundaries, then created a new volScalarField C as A/B, and wrote out C. The output file for C gives the boundaries as "type calculated" with the "value" field reflecting the correct (non-zero) value.

How do you observe the values are zero? Can you give a short representative code that demonstrates your problem?

[IvyLong]

Thanks David for your reply!

I've done a bit more tinkering around and it looks like when I divide A and B as entire fields, it does calculate the boundaries correctly. However, my denominator can and does go to zero often, so I'm dividing element by element, checking the denominator each time. In this case, the boundaries are calculated as zero (observed by checking the boundary fields in the output files).

Also, I've noticed that it happens regardless of the boundary type; in my case, an inlet, inletOutlet, and wall.

I hope this sheds a bit more light on the problem, but if you'd still like to see my solver and test case, I'd be happy to share.

Thanks for your help!

[boger]

If C is a volScalarField, for example, then it contains both a collection of the cells (cell-centered values) and the values at the boundary faces. Normally, if you calculate C = A / B, OpenFOAM will take care of updating both the cell values and boundary values for you. In this case, it sounds like you are explicitly calculating the cell values yourself and are neglecting to do anything to update the boundary values. There are lots of ways to get around this -- you could loop over each of the boundaries of C and each of the faces of each boundary, and perform your division operation (usually not recommended), you could set appropriate boundary conditions on C (e.g. by defining a 0/C file) and make sure they are updated with a call to C.correctBoundaryConditions() after you calculate the cell values, or you could do something like

dimensionedScalar BMin("BMin", B.dimensions(), SMALL);
C = A / (B + BMin)

which you'll see sometimes in OpenFOAM, for example, in the k-epsilon model calculation of the eddy viscosity.

When B is zero, how do you define C?

[marupio]

Use the stabilise function on your denominator. Then you can calculate the entire field at a time. Field operations are optimised and can be as much as 100 times faster than going element by element. Your boundary field is likely zero when you loop element by element because you neglected to perform the calculations on the boundary field. A.field() only gives the internal field. You must also use A.boundaryField[patchIndex][cellIndex].

[IvyLong]

Thanks both of you for your help! I did neglect to calculate the boundaries, but it's fixed now.

David Boger, I decided to use an alteration of the last method you suggested, adding in an additional term to prevent dividing by zero while still allowing the full field operation. In the event that B goes to zero, C does as well, so I had to get a bit creative with some of OFs scalar operations.

David Gaden, my method isn't nearly as elegant as stabilise, but it gets the job done.

Thanks again!