[strakakl]

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

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  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

// 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

[strakakl]

Hi,

any ideas about this topic?

Best Klaus

[u22]

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?

[strakakl]

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

[u22]

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?

[strakakl]

HI,

in the file volPointInterpolation.C

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

i mofified the updateMesh() member function

...

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

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

best,

Klaus

[u22]

Thank you, Klaus.
It really works!

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

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

[zhengzh5]

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

[strakakl]

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