CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

nan values in gradAlphaf (interfaceProperties.C)

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   April 7, 2014, 04:11
Default nan values in gradAlphaf (interfaceProperties.C)
Christian Butcher
Join Date: Jul 2013
Location: Japan
Posts: 85
Rep Power: 12
chrisb2244 is on a distinguished road
Foam::interfaceProperties::calculateK() is causing sigFpe crashes in an interDyMFoam based solver.

the alpha fields are generated by another class (see last post for details)
and are then copied into an incompressibleTwoPhaseMixture by createFields.H as in the standard createFields.H file of interDyMFoam.C .

On interfaceProperties object is created in the same createFields.H file, named interface. The simulation crashes during a call to interface.correct() within the pimple.firstIter() loop.

 while (pimple.loop())
            if (pimple.firstIter() || moveMeshOuterCorrectors)
                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();


                if (mesh.changing())
                    Info<< "Execution time for mesh.update() = "
                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                        << " s" << endl;

                    gh = g & mesh.C();
                    ghf = g & mesh.Cf();

                if (mesh.changing() && correctPhi)
                    // Calculate absolute flux from the mapped surface velocity
                    phi = mesh.Sf() & Uf;

                    #include "correctPhi.H"

                    // Make the flux relative to the mesh motion
                    fvc::makeRelative(phi, U);
					Pout<< "Before interface.correct()" << endl;
                    Pout<< "After interface.correct();" << endl;

                if (mesh.changing() && checkMeshCourantNo) // This never(?) triggers due to checkMeshCourantNo
                    #include "meshCourantNo.H"

            #include "alphaControls.H"
"Before interface.correct()" is printed, "After interface.correct()" is not.

A modified (for debugging purposes) interfaceProperties.C is used.
The function is below:

void Foam::interfaceProperties::calculateK()
	surfaceInterpolation::debug = 1;

    const fvMesh& mesh = alpha1_.mesh();
    const surfaceVectorField& Sf = mesh.Sf();

    // Cell gradient of alpha
    const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));

    // Interpolated face-gradient of alpha
    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
    //~ gradAlphaf.write();
    Pout<< "gradAlphaf = " << gradAlphaf << endl;
    Pout<< "mag(gradAlphaf) = "; 
    Pout<< mag(gradAlphaf) << endl;
    //~ Pout<< "deltaN_ = " << deltaN_ << endl;

    //gradAlphaf -=
    //    (mesh.Sf()/mesh.magSf())
    //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());

    // Face unit interface normal
    //~ Pout<< "gradAlphaf = " << gradAlphaf << endl; // This is throwing up nan values.

    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
    Pout<< "nHatfv = " << nHatfv << endl;
    // surfaceVectorField nHatfv
    // (
    //     (gradAlphaf + deltaN_*vector(0, 0, 1)
    //    *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
    // );
    Pout<< "About to call correctContactAngle(..) " << endl;
    correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());

    // Face unit interface normal flux
    nHatf_ = nHatfv & Sf;

    // Simple expression for curvature
    K_ = -fvc::div(nHatf_);

    // Complex expression for curvature.
    // Correction is formally zero but numerically non-zero.
    volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
    forAll(nHat.boundaryField(), patchi)
        nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];

    K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
    Pout<< "Done with calculateK()" << endl;
The only changes made are the addition of Pout<< lines, so as to enable the tracking of progress in a simplistic manner.

gradAlphaf shows several (-nan -nan -nan) values.
This then causes OF to crash during the calculation of mag(gradAlphaf) since presumably it is impossible to take the magnitude of 'not a number'...
Depending on the shape of the interface, and the initial mesh, sometimes no nan values are generated in the first timestep (0, before the real first timestep) and so the errors are applied in the first real iteration.

The alpha file generated appears to be fine, and the values of gradAlphaf output by gradAlphaf.write() (commented in the code currently) can be shown in paraview with foamToVTK -surfaceFields as glyphs, and suggest a sensible field output (ie, the glyphs are along the interface, as would be expected)

I have tried the case for both slip and zeroGradient conditions on the side walls. In the case of zeroGradient conditions on alpha, nan values can be found for the side boundary conditions. In the case of a slip boundary condition, I have no observed any in the test cases I have run (although that doesn't mean that they are impossible to generate)

Attached are the case files for the test refinement, and also the solver source files. The library for refinement can be found at the github link added below, and also linked in the earlier forum post referred to at the top of this post.

Any help or suggestions as to how I can avoid getting these nan values, or where they're coming from (more as a why, than a "because you're subtracting small values between adjacent cells in surfaceInterpolation.C or wherever, and then mapping them to faces, and they're smaller than e-3xx machine limit (side note - is e-321, e+multiple digit value sensible in any way?)) would be much appreciated.

Best wishes,
Attached Files
File Type: gz myInterDyMFoam.tar.gz (13.7 KB, 1 views)
File Type: gz refTest.tar.gz (59.6 KB, 4 views)
chrisb2244 is offline   Reply With Quote

Old   April 8, 2014, 22:52
Christian Butcher
Join Date: Jul 2013
Location: Japan
Posts: 85
Rep Power: 12
chrisb2244 is on a distinguished road
So some updates (and another request for thoughts/help)

mag(vector(1e-300 2e-300 3e-300)) = 0, so no problems from small number overflow.

mag(vector(1e-300 2e+250 3e-300)) = inf, so potentially a problem there? Not sure OF will enjoy dividing by (inf+2e-7).

Problem is caused by my mesh refinement - if I set refinementInterval to 30 in dynamicMeshDict, then the simulation will happily run for 29 iterations before crashing.

The mesh.update() function used in the solver has the same structure and implementation within dynamicRefineFvMeshHexRef4.C as for the default dynamicRefineFvMesh.C. As such, my problem is either

a) (probably, I think) that my flux corrections are somehow not properly set up, and a 3D calculation is being performed for 2D refinement, or
b) (really hoping not this) my mesh refinement is still broken (faces don't join, or volumes aren't registered, or something) or
c) (this would be nice, but I think it's unlikely) I'm missing a function call somewhere to update cell volumes or face areas or something (outside of the refinement call). I don't see an equivalent call that I'm missing for 3D refinement, so I doubt it's this.

Any ideas or suggestions would be much appreciated.

chrisb2244 is offline   Reply With Quote


interdymfoam, nan, openfoam 2.3.0

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
TimeVaryingMappedFixedValue irishdave OpenFOAM Running, Solving & CFD 32 June 16, 2021 06:55
cell values vs node values reversemermaid FLUENT 0 March 13, 2014 18:06
[General] NaN color for TecPlot data Eloise ParaView 5 February 15, 2014 14:56
error in grmsch? Adam Beamish FLOW-3D 2 September 2, 2011 04:05
Plotting raw data values Wilesco Siemens 0 January 5, 2006 05:34

All times are GMT -4. The time now is 17:25.