CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

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

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By strakakl
  • 1 Post By strakakl

Reply
 
LinkBack Thread Tools Display Modes
Old   September 28, 2012, 08:46
Default particle tracking & topo changes -> simulation crashes in 1.6-ext
  #1
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 6
strakakl is on a distinguished road
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
strakakl is offline   Reply With Quote

Old   October 15, 2012, 11:04
Default
  #2
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 6
strakakl is on a distinguished road
Hi,

any ideas about this topic?

Best Klaus
strakakl is offline   Reply With Quote

Old   December 19, 2012, 02:32
Default
  #3
u22
Member
 
Anthony Nitski
Join Date: Aug 2009
Location: Earth
Posts: 34
Rep Power: 7
u22 is on a distinguished road
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?
u22 is offline   Reply With Quote

Old   December 19, 2012, 09:06
Default
  #4
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 6
strakakl is on a distinguished road
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 likes this.
strakakl is offline   Reply With Quote

Old   December 20, 2012, 01:30
Default
  #5
u22
Member
 
Anthony Nitski
Join Date: Aug 2009
Location: Earth
Posts: 34
Rep Power: 7
u22 is on a distinguished road
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?
u22 is offline   Reply With Quote

Old   December 20, 2012, 04:35
Default
  #6
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 6
strakakl is on a distinguished road
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

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

best,

Klaus
u22 likes this.
strakakl is offline   Reply With Quote

Old   December 23, 2012, 22:42
Default
  #7
u22
Member
 
Anthony Nitski
Join Date: Aug 2009
Location: Earth
Posts: 34
Rep Power: 7
u22 is on a distinguished road
Thank you, Klaus.
It really works!
Quote:
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?
u22 is offline   Reply With Quote

Old   January 31, 2013, 16:28
Default
  #8
Member
 
Jace
Join Date: Oct 2012
Posts: 77
Rep Power: 5
zhengzh5 is on a distinguished road
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
zhengzh5 is offline   Reply With Quote

Old   February 3, 2013, 08:14
Default
  #9
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 6
strakakl is on a distinguished road
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
strakakl is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
DPM particle tracking in stirred tank parisa- FLUENT 1 August 7, 2012 12:03
particle tracking time step alireza2475 Main CFD Forum 1 June 21, 2011 00:51
Blood Damage Modelling via Particle Tracking in a Centrifugal Heart Pump scatman CFX 5 May 5, 2011 07:23
massless particle tracking problem Renold FLUENT 0 January 26, 2011 15:23
particle tracking with droplets simone marras CFX 8 September 2, 2006 04:09


All times are GMT -4. The time now is 23:15.