[Gennaro]

Dear all,

I need to solve the following equation in OpenFOAM as a part of an algorithm:

Keff = (Foam:ow(fInv/magSqr(fvc::grad(T)),scalar(1.)/scalar(3.)))/(rho*Cp);

The problem is that when grad(T) = (0 0 0), e.g on the first iteration, the ratio becomes infinite and the solver understandably crashes.

Therefore I tried to to something like:

if (magSqr(fvc::grad(T)) != 0)
{
Keff = (Foam:ow(fInv/magSqr(fvc::grad(T)),scalar(1.)/scalar(3.)))/(rho*Cp);
}
else
{
Keff = (Foam:ow(fInv/(dimensionedScalar("small", dimensionSet(0,-2,0,2,0), 1e-4)),scalar(1.)/scalar(3.)))/(rho*Cp);
}

Unfortunately this won't compile returning the error

error: no match for ‘operator==’ in ‘Foam::magSqr(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]() == 0’

Please find the whole solver in attachment and a case to test it. Please note that I'm not a programmer and I'm just starting to learn C++. Forgive me if the answer is too obvious.

Thanks and best regards

Gennaro

[tomislav_maric]

error: no match for ‘operator==’ in ‘Foam::magSqr(const  Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>  >&) [with Type = Foam::Vector<double>, PatchField =  Foam::fvPatchField, GeoMesh = Foam::volMesh]() == 0’

Even though you're not a programmer, you'll have to be come one to some extent, if you intend to extend OpenFOAM(-extend). Pun intended. The first step is reading compiler errors.

That template error says that there is no 'operator==' in GeometricField that takes '0' scalar value as the second argument.

If you want to stabilize your gradient, don't use if-then-else, you can stabilize it using SMALL (1e-15 for a double scalar). For example, check out

$WM_PROJECT_DIR/src/finiteVolume/cfdTools/general/include/setDeltaT.H

and take a look at the line:

    scalar maxDeltaTFact = maxCo/(CoNum + SMALL);

To prevent division by zero that results with a floating point exception, the denominator is summed with SMALL. Try the same for the magnitude of your gradient and loose the conditional statements, this should help.

[Gennaro]

Thanks Tomi,

Actually both != and == were attempts returning the same error message: no match for ‘operator==’or no match for ‘operator!='. So, the operator itself was not the problem.

Thanks a lot for your tip, it's very interesting. However, while SMALL is non-dimensional, I need a SMALL scalar with the right dimensions. Indeed, while with SMALL the solver would still compile, then when I run the case I get:

--> FOAM FATAL ERROR:
LHS and RHS of + have different dimensions
dimensions : [0 -2 0 2 0 0 0] + [0 0 0 0 0 0 0]

Any tips on how to apply dimensions to SMALL?

Thanks

Gennaro

[tomislav_maric]

Actually both != and == were attempts returning the same error message: no match for ‘operator==’or no match for ‘operator!='. So, the operator itself was not the problem.

Nope, this means that the GeometricField does not implement any of the two requested operators at the call sites, which results in an undefined reference error.

To give a scalar a dimension, include "dimensionedType.H" in the library code. Do not insert it into 'Make/files' for compilation - it's a class template.

Then you can instantiate a dimensionedScalar:

dimensionedScalar smallGrad ("small", gradT.dimensions(), SMALL);

Also, don't call 'fvc::grad' all over the place, it will execute the same gradient calculations over and over, unless you have set the flag for caching gradients, which is not set by default. Compute the gradient field and store it under gradT, then you can get its dimensions as shown above.

[Gennaro]

Thanks Tomi,

your help is incredibly valuable. Congratulations!

I gave it a try, but it doesn't work.

With reference to the case in attachment, I added

#include "dimensionedType.H

in laplacianFoamCustom.C, trying both inside and outside int main.

Then I added

dimensionedScalar smallGrad ("small", gradT.dimensions(), SMALL);

in createFields.H.

Finally, I changed the formula in Eeqn.H to

Keff = (Foam::pow(fInv/magSqr(fvc::grad(T) + small),scalar(1.)/scalar(3.)))/(rho*Cp);

, trying both small and SMALL.

Can you please advise?

Thanks

Genn

[Gennaro]

Hi Tomi, any news?

The error is:

createFields.H:123:40: error: ‘gradT’ was not declared in this scope

meaning that

#include "dimensionedType.H"

is not enough?

I look forward to hearing from you

Thanks

Genn

Originally Posted by tomislav_maric

Nope, this means that the GeometricField does not implement any of the two requested operators at the call sites, which results in an undefined reference error.

To give a scalar a dimension, include "dimensionedType.H" in the library code. Do not insert it into 'Make/files' for compilation - it's a class template.

Then you can instantiate a dimensionedScalar:

dimensionedScalar smallGrad ("small", gradT.dimensions(), SMALL);

Also, don't call 'fvc::grad' all over the place, it will execute the same gradient calculations over and over, unless you have set the flag for caching gradients, which is not set by default. Compute the gradient field and store it under gradT, then you can get its dimensions as shown above.

[j-avdeev]

Originally Posted by Gennaro

createFields.H:123:40: error: ‘gradT’ was not declared in this scope

Hi, Genn.

I am not sure, that it optimal way to get gradT, but you can try

surfaceVectorField gradT = fvc::interpolate(fvc::grad(T));
surfaceScalarField sumGradT = gradT.component(vector::X)+gradT.component(vector::Y)+gradT.component(vector::Z);