nan values in gradAlphaf (interfaceProperties.C)
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.
A modified (for debugging purposes) interfaceProperties.C is used.
The function is below:
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.
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.
|All times are GMT -4. The time now is 20:15.|