[chrisb2244]

Foam::interfaceProperties::calculateK() is causing sigFpe crashes in an interDyMFoam based solver.

the alpha fields are generated by another class (see http://www.cfd-online.com/Forums/ope...ned-faces.html 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();

                mesh.update();

                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;
                    interface.correct();
                    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.

https://github.com/chrisb2244/2D-AMR-OFlib

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,
Christian

[chrisb2244]

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.

Best,
Christian