View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0002286 | OpenFOAM | Bug | public | 2016-10-07 12:16 | 2016-11-18 14:36 |
| Reporter | tniemi | Assigned To | henry | ||
| Priority | normal | Severity | minor | Reproducibility | always |
| Status | resolved | Resolution | fixed | ||
| Product Version | dev | ||||
| Fixed in Version | dev | ||||
| Summary | 0002286: PatchInjectionBase: inaccurate initialization of tetFaceI and tetPtI | ||||
| Description | In 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. | ||||
| Tags | No tags attached. | ||||
|
|
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;
}
}
// ************************************************************************* //
|
|
|
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. |
|
|
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. |
|
|
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. |
|
|
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.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;
}
}
// ************************************************************************* //
|
|
|
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). |
|
|
Are you satisfied with the tests of new patchInjectionBase with randomized location? If so I will merge the change into OpenFOAM-dev. |
|
|
Yes, I have had the modified code running for a while now and I have not encountered any issues. |
|
|
Resolved in OpenFOAM-dev by commit 2001653352040f4020b598742de4e7efafdc2379 |
| 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 |
