|
[Sponsors] |
nan values in gradAlphaf (interfaceProperties.C) |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
April 7, 2014, 05:11 |
nan values in gradAlphaf (interfaceProperties.C)
|
#1 |
Member
Christian Butcher
Join Date: Jul 2013
Location: Japan
Posts: 85
Rep Power: 13 |
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. Code:
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" A modified (for debugging purposes) interfaceProperties.C is used. The function is below: Code:
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; } 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 |
|
April 8, 2014, 23:52 |
|
#2 |
Member
Christian Butcher
Join Date: Jul 2013
Location: Japan
Posts: 85
Rep Power: 13 |
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 |
|
Tags |
interdymfoam, nan, openfoam 2.3.0 |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
TimeVaryingMappedFixedValue | irishdave | OpenFOAM Running, Solving & CFD | 32 | June 16, 2021 07:55 |
cell values vs node values | reversemermaid | FLUENT | 0 | March 13, 2014 19:06 |
[General] NaN color for TecPlot data | Eloise | ParaView | 5 | February 15, 2014 15:56 |
error in grmsch? | Adam Beamish | FLOW-3D | 2 | September 2, 2011 05:05 |
Plotting raw data values | Wilesco | Siemens | 0 | January 5, 2006 06:34 |