CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   raising volScalarField to an exponent (http://www.cfd-online.com/Forums/openfoam/83003-raising-volscalarfield-exponent.html)

hyperion December 10, 2010 18:50

raising volScalarField to an exponent
 
Hello - I'm trying to modify a functionObject to output an additional field using the following chunk of code:


Info<< " Calculating mean free path field." << endl;
volScalarField mfp
(
IOobject
(
"mfp",
obr_.time().timeName(),
obr_,
IOobject::NO_READ
),
1/(8.9968e-19*rhoNMean*pow(overallT/300, .42))
);


Where overallT was defined as a volScalarField earlier in the code. The code doesn't like when I use overallT as an argument to the pow(base, exponent) function. Can I convert overallT ( a volScalarField) to something else like a scalar field? Where can I found out if a particular kind of argument is allowed?

I would greatly appreciate any advice on how to solve this.

Many Thanks

MartinB December 12, 2010 06:38

Hi Daniel,

the pow() function works for me (OpenFOAM-1.6.x) with this dummy code snippet:
Code:

        volScalarField overallT
        (
                        IOobject
                        (
                                        "dummy",
                                        runTime.timeName(),
                                        mesh,
                                        IOobject::NO_READ,
                                        IOobject::AUTO_WRITE
                        ),
                        multiPhaseProperties.mu()
        );

        volScalarField mfp
        (
        IOobject
        (
        "mfp",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ
        ),
        1/(8.9968e-19*overallT*pow(overallT/300, .42))
        );
        mfp.write();

It compiles and it writes results (although these results are nonsense, of course)... and overallT is a volScalarField...

Do you have a compilation problem or a runtime problem?

Martin

l_r_mcglashan December 12, 2010 12:47

If it's a runtime error, then it's most likely that somewhere in your overallT field you have a negative number if the error is from within the Foam::pow() function.

hyperion December 13, 2010 01:59

Thank you both for your responses. It is a runtime error and, as Laurence suggested, it is caused by division by zero. The zero values of overallT only occur on the boundary field and not in the internal field. Is there a way I can ignore boundary field values and only calculate values for the internal field?

Thanks Again

MartinB December 13, 2010 02:15

You can use
Code:

pow(max(overallT,VSMALL)/300, .42)
to avoid division by zero.

Martin

hyperion December 13, 2010 02:32

Thanks Martin - Now I get the following error message at runtime:

--> FOAM FATAL ERROR:
Arguments of max have different dimensions
dimensions : [1 2 -2 0 0 0 0] and [0 0 0 0 0 0 0]


From function max(const dimensionSet& ds1, const dimensionSet& ds2)

in file dimensionSet/dimensionSet.C at line 199.


(There's also a bunch of print stack and FPE errors following the above statement. These errors were there before adding the max statement as well.)

Could you explain what code you suggested is supposed to do? What is VSMALL? What happens when you have two arguments in the max() statement?

Many Thanks

MartinB December 13, 2010 02:49

VSMALL is a constant with a very small positive value.
max(VSMALL, otherValue) therefore guarantees a nonzero positive value to be used in a function that will fail otherwise.

To add the correct dimensions to this constant, try this:
Code:

pow(max(overallT, dimensionedScalar ("VSMALL", dimensionSet(1,2,-2,0,0,0,0), VSMALL))/300, .42)
Since the small value is further divided by 300, it might be necessary to increase its value by a factor of 300 before... just try without first and see what's happening...

Martin

l_r_mcglashan December 13, 2010 05:47

The real issue you have here is that you need to either initialise the boundary conditions, or remove them from overallT. Does mfp need boundary conditions? If not, just make the variable a scalarField. Otherwise, initialise the boundary fields.

You can access the boundaries using .boundaryField() and the cells using internalField()

hyperion December 14, 2010 14:56

Thanks Martin for the snippet of code and for your explanation of its meaning. I used that code and the application runs well now. I do obtain values for the mean free path that approach infinity at a few boundary cell faces, but I can ignore those values for now. I am treating it as a temporary fix to a larger problem, which is that the solver (dsmcFoam) attributes extreme values to particular patches in certain simulation runs. I haven't yet figured out why it does this in for certain geometry and boundary conditions.

I want to make sure I understand Laurence's remarks, so I have a few more questions. If overallT is computed from a volScalarField that has to be specified in the 0 directory, then does that mean that it is initialized and defined with boundary conditions? Can overallT be a scalarField if it's computed from a volScalarField? What do you mean I can "access" the boundary and internal fields? Do you mean as a means to read or write new values to those fields? What would be the purpose of accessing those fields?

Once Again, Many Thanks

l_r_mcglashan December 14, 2010 15:41

scalarField is just like a list of scalar values.
volScalarField is more complicated. It has a list of scalar values that represent the variable's values within the cells (the internalField()) and lists of values that represent the variable's values at each boundary patch. The files in the folder 0/ are just printouts of these volScalarFields/volVectorFields.

My point is, are you solving for mfp/overallT? If you are, you need boundary conditions. If not, then unless mfp/overallT are supplying source terms at the boundary, then they don't need boundary values and can just be defined as scalarFields.


All times are GMT -4. The time now is 19:37.