View Issue Details

IDProjectCategoryView StatusLast Update
0002286OpenFOAMBugpublic2016-11-18 14:36
Reportertniemi Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Product Versiondev 
Fixed in Versiondev 
Summary0002286: PatchInjectionBase: inaccurate initialization of tetFaceI and tetPtI
DescriptionIn PatchInjectionBase.C lines 211-212, the tetFaceI and tetPtI of the injected parcel are initialized to arbitrary values. In the comment it is said that the tracking should handle them in the first step. However, in my tests I have found out that this is not true and the wrongly initialized tetFaceI and tetPtI can cause a lot of tracking issues.

I have attached a simple fix, where tetFaceI and tetPtI are initialized to correct values with mesh.findTetFacePt. This seems to greatly reduce the amount of tracking rescues needed. In an earlier bug report (http://bugs.openfoam.org/view.php?id=1436), 30 % of parcels will get stuck in a very simple mesh. After this patch no particles get stuck even if tracking rescues are disabled.

I have not thoroughly tested this yet, but in a few test cases it has worked well.
TagsNo tags attached.

Activities

tniemi

2016-10-07 12:16

reporter  

patchInjectionBase.C (6,572 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2016 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 "patchInjectionBase.H"
#include "polyMesh.H"
#include "SubField.H"
#include "cachedRandom.H"
#include "triPointRef.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::patchInjectionBase::patchInjectionBase
(
    const polyMesh& mesh,
    const word& patchName
)
:
    patchName_(patchName),
    patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
    patchArea_(0.0),
    patchNormal_(),
    cellOwners_(),
    triFace_(),
    triToFace_(),
    triCumulativeMagSf_(),
    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
{
    if (patchId_ < 0)
    {
        FatalErrorInFunction
            << "Requested patch " << patchName_ << " not found" << nl
            << "Available patches are: " << mesh.boundaryMesh().names() << nl
            << exit(FatalError);
    }

    updateMesh(mesh);
}


Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
:
    patchName_(pib.patchName_),
    patchId_(pib.patchId_),
    patchArea_(pib.patchArea_),
    patchNormal_(pib.patchNormal_),
    cellOwners_(pib.cellOwners_),
    triFace_(pib.triFace_),
    triToFace_(pib.triToFace_),
    triCumulativeMagSf_(pib.triCumulativeMagSf_),
    sumTriMagSf_(pib.sumTriMagSf_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::patchInjectionBase::~patchInjectionBase()
{}


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

void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
{
    // Set/cache the injector cells
    const polyPatch& patch = mesh.boundaryMesh()[patchId_];
    const pointField& points = patch.points();

    cellOwners_ = patch.faceCells();

    // Triangulate the patch faces and create addressing
    DynamicList<label> triToFace(2*patch.size());
    DynamicList<scalar> triMagSf(2*patch.size());
    DynamicList<face> triFace(2*patch.size());
    DynamicList<face> tris(5);

    // Set zero value at the start of the tri area list
    triMagSf.append(0.0);

    forAll(patch, facei)
    {
        const face& f = patch[facei];

        tris.clear();
        f.triangles(points, tris);

        forAll(tris, i)
        {
            triToFace.append(facei);
            triFace.append(tris[i]);
            triMagSf.append(tris[i].mag(points));
        }
    }

    forAll(sumTriMagSf_, i)
    {
        sumTriMagSf_[i] = 0.0;
    }

    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);

    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
    Pstream::listCombineScatter(sumTriMagSf_);

    for (label i = 1; i < triMagSf.size(); i++)
    {
        triMagSf[i] += triMagSf[i-1];
    }

    // Transfer to persistent storage
    triFace_.transfer(triFace);
    triToFace_.transfer(triToFace);
    triCumulativeMagSf_.transfer(triMagSf);

    // Convert sumTriMagSf_ into cumulative sum of areas per proc
    for (label i = 1; i < sumTriMagSf_.size(); i++)
    {
        sumTriMagSf_[i] += sumTriMagSf_[i-1];
    }

    const scalarField magSf(mag(patch.faceAreas()));
    patchArea_ = sum(magSf);
    patchNormal_ = patch.faceAreas()/magSf;
    reduce(patchArea_, sumOp<scalar>());
}


void Foam::patchInjectionBase::setPositionAndCell
(
    const polyMesh& mesh,
    cachedRandom& rnd,
    vector& position,
    label& cellOwner,
    label& tetFacei,
    label& tetPti
)
{
    scalar areaFraction = rnd.globalPosition(scalar(0), patchArea_);

    if (cellOwners_.size() > 0)
    {
        // Determine which processor to inject from
        label proci = 0;
        forAllReverse(sumTriMagSf_, i)
        {
            if (areaFraction >= sumTriMagSf_[i])
            {
                proci = i;
                break;
            }
        }

        if (Pstream::myProcNo() == proci)
        {
            // Find corresponding decomposed face triangle
            label triI = 0;
            scalar offset = sumTriMagSf_[proci];
            forAllReverse(triCumulativeMagSf_, i)
            {
                if (areaFraction > triCumulativeMagSf_[i] + offset)
                {
                    triI = i;
                    break;
                }
            }

            // Set cellOwner
            label facei = triToFace_[triI];
            cellOwner = cellOwners_[facei];

            // Find random point in triangle
            const polyPatch& patch = mesh.boundaryMesh()[patchId_];
            const pointField& points = patch.points();
            const face& tf = triFace_[triI];
            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
            const point pf(tri.randomPoint(rnd));

            // Position perturbed away from face (into domain)
            const scalar a = rnd.position(scalar(0.1), scalar(0.5));
            const vector& pc = mesh.cellCentres()[cellOwner];
            const vector d = mag(pf - pc)*patchNormal_[facei];

            position = pf - a*d;

            mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);
        }
        else
        {
            cellOwner = -1;
            tetFacei = -1;
            tetPti = -1;

            // Dummy position
            position = pTraits<vector>::max;
        }
    }
    else
    {
        cellOwner = -1;
        tetFacei = -1;
        tetPti = -1;

        // Dummy position
        position = pTraits<vector>::max;
    }
}


// ************************************************************************* //
patchInjectionBase.C (6,572 bytes)   

tniemi

2016-10-07 14:26

reporter   ~0006979

Hmm, in a more complex test case the patched code crashed, because the cell of the parcel had changed due to the perturbation of its position. Changing the findTetFacePt to findCellFacePt seemed to help.

I guess this needs more testing. In any case the logic in patchInjectionBase seems quite dangerous to me.

tniemi

2016-10-07 16:11

reporter   ~0006980

With findCellFacePt the code seems to work in several more complex test cases, but there is some unnecessary overhead because usually the cell does not change and searching for it is not needed. Also I'm not sure, but I think that it might be even possible that the parcel gets pushed out of the domain when it is perturbed. In such a case even findCellFacePt won't probably work.

Perhaps a following strategy could be used?

1. Try to use findTetFacePt
2. If step 1 fails, the parcel is likely in another cell or out of the domain. Move the parcel to the center of the cell and try findTetFacePt again.
3. As a final fallback we could revert to the current method, ie. choose an arbitrary tet and hope that the tracking algorithm can rescue the parcel

In any case the current situation, where the tetFaceI and tetPtI are random and even the cell can be different than what the parcel thinks it is, is very demanding to the tracking algorithm.

henry

2016-10-07 16:46

manager   ~0006981

I agree that the current approach of hoping that the tracking algorithm will rescue the parcel if it is not in the cell it is supposed to be in is VERY poor and should be replaced. The problem with the "move the parcel to the center of the cell" is that in DEM the particle may overlap with an existing particle and in the first time-step they repel eachother at the speed of light or greater.

tniemi

2016-10-11 13:03

reporter  

patchInjectionBase_rand.H (3,996 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2016 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/>.

Class
    Foam::PatchInjectionBase

Description
    Base class for patch-based injection models.

    Class handles injecting at a random point adjacent to the patch faces to
    provide a more stochastic view of the injection process.  Patch faces are
    triangulated, and area fractions accumulated.  The fractional areas are
    then applied to determine across which face a parcel is to be injected.

SourceFiles
    PatchInjectionBase.C

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

#ifndef patchInjectionBase_H
#define patchInjectionBase_H

#include "word.H"
#include "labelList.H"
#include "scalarList.H"
#include "vectorList.H"
#include "faceList.H"

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

namespace Foam
{

// Forward class declarations
class polyMesh;
class fvMesh;
class cachedRandom;

/*---------------------------------------------------------------------------*\
                     Class patchInjectionBase Declaration
\*---------------------------------------------------------------------------*/

class patchInjectionBase
{
protected:

    // Protected data

        //- Patch name
        const word patchName_;

        //- Patch ID
        const label patchId_;

        //- Patch area - total across all processors
        scalar patchArea_;

        //- Patch face normal directions
        vectorList patchNormal_;

        //- List of cell labels corresponding to injector positions
        labelList cellOwners_;

        //- Decomposed patch faces as a list of triangles
        faceList triFace_;

        //- Addressing from per triangle to patch face
        labelList triToFace_;

        //- Cumulative triangle area per triangle face
        scalarList triCumulativeMagSf_;

        //- Cumulative area fractions per processor
        scalarList sumTriMagSf_;


public:

    // Constructors

        //- Construct from mesh and patch name
        patchInjectionBase(const polyMesh& mesh, const word& patchName);

        //- Copy constructor
        patchInjectionBase(const patchInjectionBase& pib);


    //- Destructor
    virtual ~patchInjectionBase();


    // Member Functions

        //- Update patch geometry and derived info for injection locations
        virtual void updateMesh(const polyMesh& mesh);

        //- Set the injection position and owner cell, tetFace and tetPt
        virtual void setPositionAndCell
        (
            const fvMesh& mesh,
            cachedRandom& rnd,
            vector& position,
            label& cellOwner,
            label& tetFacei,
            label& tetPti
        );
};

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

} // end namespace Foam

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

#endif

// ************************************************************************* //
patchInjectionBase_rand.H (3,996 bytes)   

tniemi

2016-10-11 13:03

reporter  

patchInjectionBase_rand.C (8,001 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2016 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 "patchInjectionBase.H"
#include "polyMesh.H"
#include "SubField.H"
#include "cachedRandom.H"
#include "triPointRef.H"
#include "volFields.H"
#include "polyMeshTetDecomposition.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::patchInjectionBase::patchInjectionBase
(
    const polyMesh& mesh,
    const word& patchName
)
:
    patchName_(patchName),
    patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
    patchArea_(0.0),
    patchNormal_(),
    cellOwners_(),
    triFace_(),
    triToFace_(),
    triCumulativeMagSf_(),
    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
{
    if (patchId_ < 0)
    {
        FatalErrorInFunction
            << "Requested patch " << patchName_ << " not found" << nl
            << "Available patches are: " << mesh.boundaryMesh().names() << nl
            << exit(FatalError);
    }

    updateMesh(mesh);
}


Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
:
    patchName_(pib.patchName_),
    patchId_(pib.patchId_),
    patchArea_(pib.patchArea_),
    patchNormal_(pib.patchNormal_),
    cellOwners_(pib.cellOwners_),
    triFace_(pib.triFace_),
    triToFace_(pib.triToFace_),
    triCumulativeMagSf_(pib.triCumulativeMagSf_),
    sumTriMagSf_(pib.sumTriMagSf_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::patchInjectionBase::~patchInjectionBase()
{}


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

void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
{
    // Set/cache the injector cells
    const polyPatch& patch = mesh.boundaryMesh()[patchId_];
    const pointField& points = patch.points();

    cellOwners_ = patch.faceCells();

    // Triangulate the patch faces and create addressing
    DynamicList<label> triToFace(2*patch.size());
    DynamicList<scalar> triMagSf(2*patch.size());
    DynamicList<face> triFace(2*patch.size());
    DynamicList<face> tris(5);

    // Set zero value at the start of the tri area list
    triMagSf.append(0.0);

    forAll(patch, facei)
    {
        const face& f = patch[facei];

        tris.clear();
        f.triangles(points, tris);

        forAll(tris, i)
        {
            triToFace.append(facei);
            triFace.append(tris[i]);
            triMagSf.append(tris[i].mag(points));
        }
    }

    forAll(sumTriMagSf_, i)
    {
        sumTriMagSf_[i] = 0.0;
    }

    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);

    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
    Pstream::listCombineScatter(sumTriMagSf_);

    for (label i = 1; i < triMagSf.size(); i++)
    {
        triMagSf[i] += triMagSf[i-1];
    }

    // Transfer to persistent storage
    triFace_.transfer(triFace);
    triToFace_.transfer(triToFace);
    triCumulativeMagSf_.transfer(triMagSf);

    // Convert sumTriMagSf_ into cumulative sum of areas per proc
    for (label i = 1; i < sumTriMagSf_.size(); i++)
    {
        sumTriMagSf_[i] += sumTriMagSf_[i-1];
    }

    const scalarField magSf(mag(patch.faceAreas()));
    patchArea_ = sum(magSf);
    patchNormal_ = patch.faceAreas()/magSf;
    reduce(patchArea_, sumOp<scalar>());
}


void Foam::patchInjectionBase::setPositionAndCell
(
    const fvMesh& mesh,
    cachedRandom& rnd,
    vector& position,
    label& cellOwner,
    label& tetFacei,
    label& tetPti
)
{
    scalar areaFraction = rnd.globalPosition(scalar(0), patchArea_);

    if (cellOwners_.size() > 0)
    {
        // Determine which processor to inject from
        label proci = 0;
        forAllReverse(sumTriMagSf_, i)
        {
            if (areaFraction >= sumTriMagSf_[i])
            {
                proci = i;
                break;
            }
        }

        if (Pstream::myProcNo() == proci)
        {
            // Find corresponding decomposed face triangle
            label triI = 0;
            scalar offset = sumTriMagSf_[proci];
            forAllReverse(triCumulativeMagSf_, i)
            {
                if (areaFraction > triCumulativeMagSf_[i] + offset)
                {
                    triI = i;
                    break;
                }
            }

            // Set cellOwner
            label facei = triToFace_[triI];
            cellOwner = cellOwners_[facei];

            // Find random point in triangle
            const polyPatch& patch = mesh.boundaryMesh()[patchId_];
            const pointField& points = patch.points();
            const face& tf = triFace_[triI];
            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
            const point pf(tri.randomPoint(rnd));

            // Position perturbed away from face (into domain)
            const scalar a = rnd.position(scalar(0.1), scalar(0.5));
            const vector& pc = mesh.cellCentres()[cellOwner];
            const vector d = mag(pf - pc)*patchNormal_[facei];

            position = pf - a*d;

            //Try to find tetFacei and tetPti in the current position
            mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);

            //Search failed, choose a random position
            if (tetFacei == -1 ||tetPti == -1)
            {
                const scalarField& V = mesh.V();

                // Construct cell tet indices
                const List<tetIndices> cellTetIs =
                    polyMeshTetDecomposition::cellTetIndices(mesh, cellOwner);

                // Construct cell tet volume fractions
                scalarList cTetVFrac(cellTetIs.size(), 0.0);
                for (label tetI = 1; tetI < cellTetIs.size() - 1; tetI++)
                {
                    cTetVFrac[tetI] =
                        cTetVFrac[tetI-1] + cellTetIs[tetI].tet(mesh).mag()/V[cellOwner];
                }
                cTetVFrac.last() = 1.0;

                // Set new particle position
                const scalar volFrac = rnd.sample01<scalar>();
                label tetI = 0;
                forAll(cTetVFrac, vfI)
                {
                    if (cTetVFrac[vfI] > volFrac)
                    {
                        tetI = vfI;
                        break;
                    }
                }
                position = cellTetIs[tetI].tet(mesh).randomPoint(rnd);
                tetFacei = cellTetIs[tetI].face();
                tetPti = cellTetIs[tetI].tetPt();
            }
        }
        else
        {
            cellOwner = -1;
            tetFacei = -1;
            tetPti = -1;

            // Dummy position
            position = pTraits<vector>::max;
        }
    }
    else
    {
        cellOwner = -1;
        tetFacei = -1;
        tetPti = -1;

        // Dummy position
        position = pTraits<vector>::max;
    }
}


// ************************************************************************* //
patchInjectionBase_rand.C (8,001 bytes)   

tniemi

2016-10-11 13:06

reporter   ~0007000

Yes, with DEM there could be problems... How about randomizing the position instead of using the cell centre in cases where the normal perturbation does not work? Of course this doesn't prevent overlap but it would make it less probable.

I have attached a modified patchInjectionBase with randomizing code taken from CellZoneInjection. I have not yet tested it thoroughly, but in a few test cases (steady state, no DEM) it has worked ok with a clear reduction in the amount of stalled parcels (=parcels requiring more than 1 sequential rescue).

henry

2016-11-17 16:21

manager   ~0007226

Are you satisfied with the tests of new patchInjectionBase with randomized location? If so I will merge the change into OpenFOAM-dev.

tniemi

2016-11-18 06:57

reporter   ~0007236

Yes, I have had the modified code running for a while now and I have not encountered any issues.

henry

2016-11-18 14:36

manager   ~0007242

Resolved in OpenFOAM-dev by commit 2001653352040f4020b598742de4e7efafdc2379

Issue History

Date Modified Username Field Change
2016-10-07 12:16 tniemi New Issue
2016-10-07 12:16 tniemi File Added: patchInjectionBase.C
2016-10-07 14:26 tniemi Note Added: 0006979
2016-10-07 16:11 tniemi Note Added: 0006980
2016-10-07 16:46 henry Note Added: 0006981
2016-10-11 13:03 tniemi File Added: patchInjectionBase_rand.H
2016-10-11 13:03 tniemi File Added: patchInjectionBase_rand.C
2016-10-11 13:06 tniemi Note Added: 0007000
2016-11-17 16:21 henry Note Added: 0007226
2016-11-18 06:57 tniemi Note Added: 0007236
2016-11-18 14:36 henry Assigned To => henry
2016-11-18 14:36 henry Status new => resolved
2016-11-18 14:36 henry Resolution open => fixed
2016-11-18 14:36 henry Fixed in Version => dev
2016-11-18 14:36 henry Note Added: 0007242