View Issue Details

IDProjectCategoryView StatusLast Update
0001310OpenFOAMBugpublic2015-03-07 15:29
ReportercraigWhite Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityrandom
Status resolvedResolutionunable to reproduce 
PlatformLinuxOSUbuntuOS Version12.04
Summary0001310: Lagrangian Particle Tracking Hangs
DescriptionDSMC solver sometimes hangs during the lagrangian particle tracking. The issue is only really seen when using large numbers of particles and happens at random intervals - I have changed the seeding for the random number generator so I wouldn't expect it to hang at the same place every time. I'm not sure what function exactly causes the issue. There is no error output, the process continues to run, but the simulation cannot proceed.

Sometimes it happens soon after a write interval, and then this write interval cannot be reconstructed as there is a particle with a position that does not match any cell - it is actually quite a relatively large distance outside of the mesh when comparing the particle position to the mesh bounding box.
Steps To ReproduceLots of particles (30 million), run in parallel, include wall surfaces and freestream boundaries.
Additional InformationAlthough I using 2.1.x, I have tried adding the fixes from previous reports of what seem to be the same problem: http://www.openfoam.org/mantisbt/view.php?id=1059

My issue still remains.
TagsNo tags attached.

Activities

wyldckat

2014-06-28 22:38

updater   ~0003143

Last edited: 2014-06-28 22:42

If you could provide an example case where this can be reproduced, or at least instructions on how to configure one of OpenFOAM's tutorial cases to reproduce the same issue, it would make it a lot easier/faster to diagnose this.

In addition, knowing how many cores to use and the decomposition method used for the parallel run, would also help.

user924

2014-09-21 18:21

 

solidParticle.C (6,910 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "solidParticleCloud.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTemplateTypeNameAndDebug(Cloud<solidParticle>, 0);
}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //


bool Foam::solidParticle::move( trackingData& td,
                                const scalar trackTime
                              )
{
    td.switchProcessor = false;
    td.keepParticle = true;

    const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();

    scalar tEnd = (1.0 - stepFraction())*trackTime;
    scalar dtMax = tEnd;

    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
    {
        if (debug)
        {
            Info<< "Time = " << mesh_.time().timeName()
                << " trackTime = " << trackTime
                << " tEnd = " << tEnd
                << " steptFraction() = " << stepFraction() << endl;
        }

        // set the lagrangian time-step
        scalar dt = min(dtMax, tEnd);

        // remember which cell the parcel is in
        // since this will change if a face is hit
        label cellI = cell();
  
        //- Track particle to a given position and returns 1.0 if the
        //  trajectory is completed without hitting a face otherwise
        //  stops at the face and returns the fraction of the trajectory
        //  completed.
        //  on entry 'stepFraction()' should be set to the fraction of the
        //  time-step at which the tracking starts.
//        Pout << ".";
        dt *= trackToFace(position() + dt*U_, td);

/*        if(dt<=0) 
        {
            Pout << "S";
        }*/

        tEnd -= dt;
        stepFraction() = 1.0 - tEnd/trackTime;
  
        cellPointWeight cpw(mesh_, position(), cellI, face());
        vector Uc = td.UInterp().interpolate(cpw);
      
        // Add in here the interpolation for p, p_rgh, alpha, which we use in ParaFoam to colour the dots
        p_ = td.pInterp().interpolate(cpw);
        p_rgh_ = td.p_rghInterp().interpolate(cpw);
        alpha_ = td.alphaInterp().interpolate(cpw);
        // -----------
        scalar rhoc = td.rhoInterp().interpolate(cpw);    // density of fluid surrounding 
        scalar rhop  = td.cloud().rhop();
   
        // Lets assume that if the particle are smaller than zero, then we merely want to track the flow
        // So we'll just skip the following
        if(d_>0)
        {
        	if (debug) Info << "Real Particles" <<endl;
            scalar nuc  = td.nuInterp().interpolate(cpw);

            // Seems to be based on: http://www.freestudy.co.uk/fluid%20mechanics/saet3a.pdf
            scalar magUr = mag(Uc - U_);
            
            scalar ReFunc = 1.0;
            scalar Re = magUr*d_/nuc;   // where nu is kinematic viscosity

            if (Re > 0.01)
            {
                // This looks like a coefficient from "Clift et al, 1978", where they give
                // Cd = 24 * (1+0.15 Re^0.687) / Re
                ReFunc = 1.0 + 0.15*pow(Re, 0.687);
            }

            // Compute the drag coefficient
            scalar Dc = (24.0*nuc/d_)*ReFunc*(3.0/4.0)*(rhoc/(d_*rhop));

            // Compute new velocity of particle. Initial velocity + acceleration from drag
            // This is a rearrangement to get U1, for
            //   U1 = U0 + dt*Dc*(Uc-U1) + dt*(1-rhoc/rhop)*g

            U_ = (U_ + dt*(Dc*Uc + (1.0 - rhoc/rhop)*td.g()))/(1.0 + dt*Dc);
        }
        else
        {
        	if (debug) Info << "Trace particles" <<endl;

            // As we are merely tracking the stream flow, but we want to retain acceleration due to 
            // gravity on particle as they go into the air. But I'm assuming these particle reach terminal 
            // velocity immediately, i.e. there is no retained acceleration from a previous time step. 
            U_ = Uc + (1.0 - rhoc/rhop)*td.g()*dt;
        }

        if (onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[patch(face())]))
            {
                td.switchProcessor = true;
            }
            else 
            {
                // Particle has reached the edge of the domain, so it must vanish
                td.keepParticle=false;
            }
        }
    }

    return td.keepParticle;
}


bool Foam::solidParticle::hitPatch
(
    const polyPatch&,
    trackingData&,
    const label,
    const scalar,
    const tetIndices&
)
{
    return false;
}


void Foam::solidParticle::hitProcessorPatch
(
    const processorPolyPatch&,
    trackingData& td
)
{
    td.switchProcessor = true;
}


void Foam::solidParticle::hitWallPatch(
                                        const wallPolyPatch& wpp,
                                        trackingData& td,
                                        const tetIndices& tetIs
                                      )
{
    vector nw = tetIs.faceTri(mesh_).normal();
    nw /= mag(nw);

    scalar Un = U_ & nw;
    vector Ut = U_ - Un*nw;

    if (Un > 0)
    {
        U_ -= (1.0 + td.cloud().e())*Un*nw;
    }

    U_ -= td.cloud().mu()*Ut;
}


void Foam::solidParticle::hitPatch( const polyPatch&,
                                    trackingData& td )
{
    td.keepParticle = false;
}


void Foam::solidParticle::transformProperties (const tensor& T)
{
    particle::transformProperties(T);
    U_ = transform(T, U_);
}


void Foam::solidParticle::transformProperties(const vector& separation)
{
    particle::transformProperties(separation);
}


Foam::scalar Foam::solidParticle::wallImpactDistance(const vector&) const
{
    if(d_<=0) return 0.0005;  // Some small number
    return 0.5*d_;
}


// ************************************************************************* //
solidParticle.C (6,910 bytes)   

user924

2014-09-21 18:22

  ~0003236

We experience the same problem with LPT within waveFoam. For this integration we copied the files from src/lagrangian/solidParticle. Because we wanted to use the particles as tracer for the fluid (rather than physical particles) we modified solidParticle.C amongst other bits. Everything worked fine, except the model would randomly hang.

In debugged this, I believe that the error is in the last part of the function solidParticle.move, where it detects if you've reached a processor boundary, but doesn't address the case when the particle has left the computational domain completely.

I therefore THINK (I would stress I'm still trying to figure out how everything ties together) that the lines:

        if (onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[patch(face())]))
            {
                td.switchProcessor = true;
            }
        }

should really be

        if (onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[patch(face())]))
            {
                td.switchProcessor = true;
            }
            else
            {
                // Particle has reached the edge of the domain, so it must vanish
                td.keepParticle=false;
            }
        }


If tracking this down, we also monitored the line
        dt *= trackToFace(position() + dt*U_, td);

What still worries me is that often dt<=0 at the end of this. Which doesn't all together seems wise. So we've added

if(dt<=0) dt=SMALL

---

The final change we've made, given that I know a lot of people are using LPT simply to trace the movement of the underlying fluid, we have added the option for a particle size of 0, however we left in gravitational accelerations (originally to get around the above bug), but also because we don't want to track particles within the air - so this approach gradually cause particles that leave our water surface to drop back again. At least that the hope.

I'd be interested in feedback on all these alterations. In particular are they appropriate.


Thanks Will

user924

2014-09-25 15:03

  ~0003237

Not sure if this is a bug or a limitation. But LGT doesn't seem to work with a Dynamic Mesh. I assume the particles find that the mesh to which they are associated has changed and are somewhat confused by this.

henry

2015-03-07 15:23

manager   ~0003999

Most but not all the cloud types supported dynamic mesh changes. I added this support to DSMC is OpenFOAM-dev:

commit 218875bf401389262da828a896f7e7633df1480d
Author: Henry <Henry>
Date: Fri Feb 20 17:24:14 2015 +0000

    DSMC: Rationalization and addition of mapping support

Could you check and let me know if this resolves the problem.

henry

2015-03-07 15:28

manager   ~0004000

The issue with dynamic mesh is separate to the original report so if you have some feedback on it please open a separate as I will close this orphaned report.

Note that your proposal to remove particles when they hit boundaries may not be appropriate depending on application, for example if you want to accumulate the particles on walls.

Issue History

Date Modified Username Field Change
2014-06-02 11:23 craigWhite New Issue
2014-06-28 22:38 wyldckat Note Added: 0003143
2014-06-28 22:42 wyldckat Note Edited: 0003143
2014-09-21 18:21 user924 File Added: solidParticle.C
2014-09-21 18:22 user924 Note Added: 0003236
2014-09-25 15:03 user924 Note Added: 0003237
2015-03-07 15:23 henry Note Added: 0003999
2015-03-07 15:28 henry Note Added: 0004000
2015-03-07 15:28 henry Status new => resolved
2015-03-07 15:28 henry Resolution open => unable to reproduce
2015-03-07 15:28 henry Assigned To => henry