CFD Online Discussion Forums - particle tracking & topo changes -> simulation crashes in 1.6-ext

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - particle tracking & topo changes -> simulation crashes in 1.6-ext (https://www.cfd-online.com/Forums/openfoam-programming-development/107487-particle-tracking-topo-changes-simulation-crashes-1-6-ext.html)

particle tracking & topo changes -> simulation crashes in 1.6-ext

Hi,

I try to track particles is a domain with changing topology in OpenFOAM-1.6-ext. Therefore i use a small test case where i add a single kinematic parcel in a cylindrical domain with increasing length in axial direction.
Running the case gives

Code:

/*---------------------------------------------------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |

|  \\    /   O peration     | Version:  1.6-ext                               |

|   \\  /    A nd           | Web:      www.extend-project.de                 |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

Build  : 1.6-ext-c04489736fa1

Exec   : ipim1TransientSimpleDyMFoam

Date   : Sep 28 2012

Time   : 13:43:49

Host   : pimpa

PID    : 29408

Case   : /media/data_xfs/straka/OpenFoam/pimpa-1.6-ext/run/GGI/movingTopo3DBlockMesh

nProcs : 1

SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).



// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Create time



Create dynamic mesh for time = 0



Selecting dynamicFvMesh movingBodyTopoFvMesh

Selecting solid-body motion function translation

Adding 1 mesh modifiers

Updating vertex markup

Reading field p



Reading field U



Reading/calculating face flux field phi



Selecting incompressible transport model Newtonian

Selecting turbulence model type laminar



Reading g

Constructing kinematicCloud kinematicCloud1

--> FOAM Warning : 

    From function Cloud<ParticleType>::initCloud(const bool checkClass)

    in file /home/pimpa/OpenFOAM/OpenFOAM-1.6-ext/src/lagrangian/basic/lnInclude/CloudIO.C at line 125

    Cannot read particle positions file 

    "/media/data_xfs/straka/OpenFoam/pimpa-1.6-ext/run/GGI/movingTopo3DBlockMesh/0/lagrangian/kinematicCloud1"

    assuming the initial cloud contains 0 particles.

Selecting DispersionModel none

Selecting DragModel SphereDrag

Selecting InjectionModel ManualInjection

    Constructing 3-D injection

Selecting pdfType RosinRammler

Selecting PatchInteractionModel StandardWallInteraction

Selecting PostProcessingModel PatchPostProcessing

Selecting U IntegrationScheme Euler



Starting time loop



Volume: new = 2.31928e-06 old = 2.31928e-06 change = 0

Courant Number mean: 0 max: 3 velocity magnitude: 0.01

deltaT = 0.005

Time = 0.005



solidBodyMotionFunctions::translation::transformation(): Time = 0.005 velocity: ((0 0 0.001) (0 (0 0 0)))

Executing mesh motion

Evolving kinematicCloud1



--> Cloud: kinematicCloud1

    Added 1 new parcels



Cloud: kinematicCloud1

    Total number of parcels added   = 1

    Total mass introduced           = 0.0001

    Current number of parcels       = 1

    Current mass in system          = 0.0001



Cloud: kinematicCloud1

    Total number of parcels added   = 1

    Total mass introduced           = 0.0001

    Current number of parcels       = 1

    Current mass in system          = 0.0001

volume continuity errors : volume = 2.32044e-06, max error = 1.22968e-12, sum local = 6.07344e-17, global = -5.32805e-18

DILUPBiCG:  Solving for Ux, Initial residual = 0.936055, Final residual = 4.25267e-06, No Iterations 4

DILUPBiCG:  Solving for Uy, Initial residual = 0.936137, Final residual = 4.23359e-06, No Iterations 4

DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 2.15755e-06, No Iterations 5

DICPCG:  Solving for p, Initial residual = 1, Final residual = 9.62349e-07, No Iterations 62

time step continuity errors : sum local = 4.80934e-09, global = -5.37518e-11, cumulative = -5.37518e-11

DILUPBiCG:  Solving for Ux, Initial residual = 0.581895, Final residual = 9.38825e-06, No Iterations 4

DILUPBiCG:  Solving for Uy, Initial residual = 0.581895, Final residual = 9.45249e-06, No Iterations 4

DILUPBiCG:  Solving for Uz, Initial residual = 0.882132, Final residual = 7.38817e-07, No Iterations 5

DICPCG:  Solving for p, Initial residual = 0.0147739, Final residual = 5.86599e-07, No Iterations 48

time step continuity errors : sum local = 6.93377e-08, global = -2.15609e-10, cumulative = -2.69361e-10

DILUPBiCG:  Solving for Ux, Initial residual = 0.260902, Final residual = 3.83719e-06, No Iterations 4

DILUPBiCG:  Solving for Uy, Initial residual = 0.260901, Final residual = 3.87016e-06, No Iterations 4

DILUPBiCG:  Solving for Uz, Initial residual = 0.547004, Final residual = 5.1555e-06, No Iterations 4

DICPCG:  Solving for p, Initial residual = 0.0494603, Final residual = 8.24994e-07, No Iterations 51

time step continuity errors : sum local = 7.54477e-08, global = -9.93111e-10, cumulative = -1.26247e-09

DILUPBiCG:  Solving for Ux, Initial residual = 0.0871615, Final residual = 1.73518e-06, No Iterations 4

DILUPBiCG:  Solving for Uy, Initial residual = 0.0871609, Final residual = 1.73733e-06, No Iterations 4

DILUPBiCG:  Solving for Uz, Initial residual = 0.135341, Final residual = 3.26356e-06, No Iterations 4

DICPCG:  Solving for p, Initial residual = 0.144846, Final residual = 6.42283e-07, No Iterations 56

time step continuity errors : sum local = 1.08982e-08, global = -1.89543e-11, cumulative = -1.28143e-09

ExecutionTime = 1.49 s  ClockTime = 2 s



Volume: new = 2.32044e-06 old = 2.31928e-06 change = 1.15964e-09

Courant Number mean: 0.06297 max: 0.450956 velocity magnitude: 0.0162153

deltaT = 0.00332627

Time = 0.00832627



Floating point exception

the simulation crashes with a floating point exception :(
After some debugging i found that the exception occurs during the topo change in the makeWeights() member function of the class volPointInterpolation

Code:

// Calculate inverse distances between cell centres and points

    // and store in weighting factor array

    forAll(points, pointi)

    {

        scalarList& pw = pointWeights_[pointi];

        const labelList& pcp = pointCells[pointi];

        forAll(pcp, pointCelli)

        {

             pw[pointCelli] =

                1.0/mag(points[pointi] - cellCentres[pcp[pointCelli]]);

             sumWeights[pointi] += pw[pointCelli];

        }

    }

because one of the cell centres is located at position (0 1.00974e+272 -8.07794e+272) which is far out of the bounding box of the domain. The cell with the above mentioned cell centre is added by the topo changer.

I have no idea how to fix this issue, maybe someone can help me with that. I put the solver (based on the transientSimpleDyMFoam, added particle tracking) and the test case in my dropbox

solver:
https://dl.dropbox.com/u/86131187/ip...pleDyMFoam.zip

test case:
https://dl.dropbox.com/u/86131187/mo...DBlockMesh.zip

Have a nice day,

Klaus

Hi,

any ideas about this topic?

Best Klaus

Klaus,

I have same troubles with dieselEngineDyMFoam in 1.6-ext.
Look at here http://www.cfd-online.com/Forums/ope...tml#post398168

Did you solve a problem? Any ideas?

Hi,

as far as i understand the topo change procedure is the following:

1.) Detect if a new cell layer have to be created.
2.) If yes create the layer. At this time the new cell have a volume equal to zero.
3.) perform mesh movement, this means "extrude" the new cell layer in movement direction. After this step the cells of the new layer have a well defined volume > 0.

I think the problem is that after step (2) the volPointInterpolation is performed. But as mentioned above the cells are not well defined at this time (how to calculate the cell centre of a cell with V=0?). Thus volPointInterpolation may crash with a FPE.

My workaround is to call the makeWeights() member function of the volPointInterpolation (volPointInterpolation.C) only if all cell centres are within the mesh bounding box.

This should not be a problem in layer addition/removal topo changes as the volPointInterpolation is called again in step 3 where the cell centres are correctly defined.

If you have comments on that solution please tell me, maybe (or i'm sure) there is a better way to solve the problem.

Best,

Klaus

Thanks for comments, Klaus.
I also suspect that what you're saying. But unfortunately I'm not good in the OpenFOAM code. How your idea will look in C++? What files needed to modify?

HI,

in the file volPointInterpolation.C

/OpenFOAM/OpenFOAM-1.6-ext/src/finiteVolume/interpolation/volPointInterpolation

i mofified the updateMesh() member function

Code:

...



bool volPointInterpolation::updateMesh(const mapPolyMesh&) const

{

    

    // During layer addition topo change we get cell centres

    // far outside of the mesh bounding box, thus makeWeights()

    // will crash due a floating point exception.

    // Hack to overcome this issue: do not execute the makeWeights

    // member function if some cell centres are outside of the

    // mesh bounding box. Should not be a problem in layer addition

    // topo changes because volPointInterolation is calculated again

    // after mesh movement (movePoints()) but may have side effects

    // in other cases!   

    // stra 23/Oct/12

    const boundBox &bb = mesh().bounds();

    const vectorField& cellCentres = mesh().cellCentres();

    bool updateWeights = true;

    forAll(cellCentres,celli)

    {

        if (!bb.contains(cellCentres[celli]))

        {

            updateWeights = false;

        }

    }



    if (updateWeights) //stra ende

        makeWeights();



    // Updated for MeshObject handling

    // HJ, 30/Aug/2010

    const_cast<pointPatchInterpolation&>(boundaryInterpolator_).updateMesh();



    return true;

}



...

but be aware - modifying the src in this way is not the best idea because performing an update of OF-1.6-ext via git will remove your changes!!
This is a quick and dirty hack, but i have not thought about a better solution :eek:

It may also be not very efficient to check all cell centres - as i said, it's a hack ;)

best,

Klaus

Thank you, Klaus.
It really works!

Quote:

With regard to updating 1.6 - perhaps we can put a bug report and expect fix in future versions?

Hi, I am having the same problem in OpenFOAM-2.0.x.

The member function has changed since 1.6-ext, so I try to pull something very similar.

Old code is:
void volPointInterpolation::updateMesh()
{
makeWeights();
makePatchPatchAddressing();
}

changed to:

void volPointInterpolation::updateMesh()
{
const boundBox &bb = mesh().bounds();
const vectorField& cellCentres = mesh().cellCentres();
bool updateWeights = true;
forAll(cellCentres,celli)
{
if (!bb.contains(cellCentres[celli]))
{
updateWeights = false;
}
}
if (updateWeights){
makeWeights();
makePatchPatchAddressing();
}
}

But the problem is still there. I would still run into makeWeight error. Did I understand the idea wrong? any help is appreciated.

Thanks,

J

HI,

sorry no idea, code looks fine. Have you tried to put debug statements into the code to ensure everything runs fine?
By the way, i tried particle tracking and topo changes in OF2.1.x and didn't recognize above mentioned problems.

Best,
Klaus