PairCollision.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
    Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "PairCollision.H"
#include "PairModel.H"
#include "WallModel.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::scalar Foam::PairCollision<CloudType>::cosPhiMinFlatWall = 1 - SMALL;
 
template<class CloudType>
Foam::scalar Foam::PairCollision<CloudType>::flatWallDuplicateExclusion =
    sqrt(3*SMALL);
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template<class CloudType>
void Foam::PairCollision<CloudType>::preInteraction()
{
    // Set accumulated quantities to zero
    for (typename CloudType::parcelType& p : this->owner())
    {
        p.f() = Zero;
 
        p.torque() = Zero;
    }
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::parcelInteraction()
{
    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
    label startOfRequests = Pstream::nRequests();
 
    il_.sendReferredData(this->owner().cellOccupancy(), pBufs);
 
    realRealInteraction();
 
    il_.receiveReferredData(pBufs, startOfRequests);
 
    realReferredInteraction();
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::realRealInteraction()
{
    // Direct interaction list (dil)
    const labelListList& dil = il_.dil();
 
    typename CloudType::parcelType* pA_ptr = nullptr;
    typename CloudType::parcelType* pB_ptr = nullptr;
 
    List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
        this->owner().cellOccupancy();
 
    forAll(dil, realCelli)
    {
        // Loop over all Parcels in cell A (a)
        forAll(cellOccupancy[realCelli], a)
        {
            pA_ptr = cellOccupancy[realCelli][a];
 
            forAll(dil[realCelli], interactingCells)
            {
                List<typename CloudType::parcelType*> cellBParcels =
                    cellOccupancy[dil[realCelli][interactingCells]];
 
                // Loop over all Parcels in cell B (b)
                forAll(cellBParcels, b)
                {
                    pB_ptr = cellBParcels[b];
 
                    evaluatePair(*pA_ptr, *pB_ptr);
                }
            }
 
            // Loop over the other Parcels in cell A (aO)
            forAll(cellOccupancy[realCelli], aO)
            {
                pB_ptr = cellOccupancy[realCelli][aO];
 
                // Do not double-evaluate, compare pointers, arbitrary
                // order
                if (pB_ptr > pA_ptr)
                {
                    evaluatePair(*pA_ptr, *pB_ptr);
                }
            }
        }
    }
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::realReferredInteraction()
{
    // Referred interaction list (ril)
    const labelListList& ril = il_.ril();
 
    List<IDLList<typename CloudType::parcelType>>& referredParticles =
        il_.referredParticles();
 
    List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
        this->owner().cellOccupancy();
 
    // Loop over all referred cells
    forAll(ril, refCelli)
    {
        IDLList<typename CloudType::parcelType>& refCellRefParticles =
            referredParticles[refCelli];
 
        const labelList& realCells = ril[refCelli];
 
        // Loop over all referred parcels in the referred cell
 
        for
        (
            typename CloudType::parcelType& referredParcel
          : refCellRefParticles
        )
        {
            // Loop over all real cells in that the referred cell is
            // to supply interactions to
 
            forAll(realCells, realCelli)
            {
                List<typename CloudType::parcelType*> realCellParcels =
                    cellOccupancy[realCells[realCelli]];
 
                forAll(realCellParcels, realParcelI)
                {
                    evaluatePair
                    (
                        *realCellParcels[realParcelI],
                        referredParcel
                    );
                }
            }
        }
    }
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::wallInteraction()
{
    const polyMesh& mesh = this->owner().mesh();
 
    const labelListList& dil = il_.dil();
 
    const labelListList& directWallFaces = il_.dwfil();
 
    const labelList& patchID = mesh.boundaryMesh().patchID();
 
    const volVectorField& U = mesh.lookupObject<volVectorField>(il_.UName());
 
    List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
        this->owner().cellOccupancy();
 
    // Storage for the wall interaction sites
    DynamicList<point> flatSitePoints;
    DynamicList<scalar> flatSiteExclusionDistancesSqr;
    DynamicList<WallSiteData<vector>> flatSiteData;
    DynamicList<point> otherSitePoints;
    DynamicList<scalar> otherSiteDistances;
    DynamicList<WallSiteData<vector>> otherSiteData;
    DynamicList<point> sharpSitePoints;
    DynamicList<scalar> sharpSiteExclusionDistancesSqr;
    DynamicList<WallSiteData<vector>> sharpSiteData;
 
    forAll(dil, realCelli)
    {
        // The real wall faces in range of this real cell
        const labelList& realWallFaces = directWallFaces[realCelli];
 
        // Loop over all Parcels in cell
        forAll(cellOccupancy[realCelli], cellParticleI)
        {
            flatSitePoints.clear();
            flatSiteExclusionDistancesSqr.clear();
            flatSiteData.clear();
            otherSitePoints.clear();
            otherSiteDistances.clear();
            otherSiteData.clear();
            sharpSitePoints.clear();
            sharpSiteExclusionDistancesSqr.clear();
            sharpSiteData.clear();
 
            typename CloudType::parcelType& p =
                *cellOccupancy[realCelli][cellParticleI];
 
            const point pos(p.position());
 
            scalar r = wallModel_->pREff(p);
 
            // real wallFace interactions
 
            forAll(realWallFaces, realWallFacei)
            {
                label realFacei = realWallFaces[realWallFacei];
 
                pointHit nearest = mesh.faces()[realFacei].nearestPoint
                (
                    pos,
                    mesh.points()
                );
 
                if (nearest.distance() < r)
                {
                    const vector normal =
                        normalised(mesh.faceAreas()[realFacei]);
 
                    const vector& nearPt = nearest.rawPoint();
 
                    const vector pW = normalised(nearPt - pos);
 
                    const scalar normalAlignment = normal & pW;
 
                    // Find the patchIndex and wallData for WallSiteData object
                    label patchi = patchID[realFacei - mesh.nInternalFaces()];
 
                    label patchFacei =
                        realFacei - mesh.boundaryMesh()[patchi].start();
 
                    WallSiteData<vector> wSD
                    (
                        patchi,
                        U.boundaryField()[patchi][patchFacei]
                    );
 
                    if (normalAlignment > cosPhiMinFlatWall)
                    {
                        // Guard against a flat interaction being
                        // present on the boundary of two or more
                        // faces, which would create duplicate contact
                        // points. Duplicates are discarded.
                        if
                        (
                            !duplicatePointInList
                            (
                                flatSitePoints,
                                nearPt,
                                sqr(r*flatWallDuplicateExclusion)
                            )
                        )
                        {
                            flatSitePoints.append(nearPt);
 
                            flatSiteExclusionDistancesSqr.append
                            (
                                sqr(r) - sqr(nearest.distance())
                            );
 
                            flatSiteData.append(wSD);
                        }
                    }
                    else
                    {
                        otherSitePoints.append(nearPt);
 
                        otherSiteDistances.append(nearest.distance());
 
                        otherSiteData.append(wSD);
                    }
                }
            }
 
            // referred wallFace interactions
 
            // The labels of referred wall faces in range of this real cell
            const labelList& cellRefWallFaces = il_.rwfilInverse()[realCelli];
 
            forAll(cellRefWallFaces, rWFI)
            {
                label refWallFacei = cellRefWallFaces[rWFI];
 
                const referredWallFace& rwf =
                    il_.referredWallFaces()[refWallFacei];
 
                const pointField& pts = rwf.points();
 
                pointHit nearest = rwf.nearestPoint(pos, pts);
 
                if (nearest.distance() < r)
                {
                    const vector normal = rwf.unitNormal(pts);
 
                    const vector& nearPt = nearest.rawPoint();
 
                    const vector pW = normalised(nearPt - pos);
 
                    const scalar normalAlignment = normal & pW;
 
                    // Find the patchIndex and wallData for WallSiteData object
 
                    WallSiteData<vector> wSD
                    (
                        rwf.patchIndex(),
                        il_.referredWallData()[refWallFacei]
                    );
 
                    if (normalAlignment > cosPhiMinFlatWall)
                    {
                        // Guard against a flat interaction being
                        // present on the boundary of two or more
                        // faces, which would create duplicate contact
                        // points. Duplicates are discarded.
                        if
                        (
                            !duplicatePointInList
                            (
                                flatSitePoints,
                                nearPt,
                                sqr(r*flatWallDuplicateExclusion)
                            )
                        )
                        {
                            flatSitePoints.append(nearPt);
 
                            flatSiteExclusionDistancesSqr.append
                            (
                                sqr(r) - sqr(nearest.distance())
                            );
 
                            flatSiteData.append(wSD);
                        }
                    }
                    else
                    {
                        otherSitePoints.append(nearPt);
 
                        otherSiteDistances.append(nearest.distance());
 
                        otherSiteData.append(wSD);
                    }
                }
            }
 
            // All flat interaction sites found, now classify the
            // other sites as being in range of a flat interaction, or
            // a sharp interaction, being aware of not duplicating the
            // sharp interaction sites.
 
            // The "other" sites need to evaluated in order of
            // ascending distance to their nearest point so that
            // grouping occurs around the closest in any group
 
            labelList sortedOtherSiteIndices
            (
                sortedOrder(otherSiteDistances)
            );
 
            for (const label orderedIndex : sortedOtherSiteIndices)
            {
                const point& otherPt = otherSitePoints[orderedIndex];
 
                if
                (
                    !duplicatePointInList
                    (
                        flatSitePoints,
                        otherPt,
                        flatSiteExclusionDistancesSqr
                    )
                )
                {
                    // Not in range of a flat interaction, must be a
                    // sharp interaction.
 
                    if
                    (
                        !duplicatePointInList
                        (
                            sharpSitePoints,
                            otherPt,
                            sharpSiteExclusionDistancesSqr
                        )
                    )
                    {
                        sharpSitePoints.append(otherPt);
 
                        sharpSiteExclusionDistancesSqr.append
                        (
                            sqr(r) - sqr(otherSiteDistances[orderedIndex])
                        );
 
                        sharpSiteData.append(otherSiteData[orderedIndex]);
                    }
                }
            }
 
            evaluateWall
            (
                p,
                flatSitePoints,
                flatSiteData,
                sharpSitePoints,
                sharpSiteData
            );
        }
    }
}
 
 
template<class CloudType>
bool Foam::PairCollision<CloudType>::duplicatePointInList
(
    const DynamicList<point>& existingPoints,
    const point& pointToTest,
    scalar duplicateRangeSqr
) const
{
    forAll(existingPoints, i)
    {
        if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr)
        {
            return true;
        }
    }
 
    return false;
}
 
 
template<class CloudType>
bool Foam::PairCollision<CloudType>::duplicatePointInList
(
    const DynamicList<point>& existingPoints,
    const point& pointToTest,
    const scalarList& duplicateRangeSqr
) const
{
    forAll(existingPoints, i)
    {
        if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr[i])
        {
            return true;
        }
    }
 
    return false;
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::postInteraction()
{
    // Delete any collision records where no collision occurred this step
 
    for (typename CloudType::parcelType& p : this->owner())
    {
        p.collisionRecords().update();
    }
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::evaluatePair
(
    typename CloudType::parcelType& pA,
    typename CloudType::parcelType& pB
) const
{
    pairModel_->evaluatePair(pA, pB);
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::evaluateWall
(
    typename CloudType::parcelType& p,
    const List<point>& flatSitePoints,
    const List<WallSiteData<vector>>& flatSiteData,
    const List<point>& sharpSitePoints,
    const List<WallSiteData<vector>>& sharpSiteData
) const
{
    wallModel_->evaluateWall
    (
        p,
        flatSitePoints,
        flatSiteData,
        sharpSitePoints,
        sharpSiteData
    );
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::PairCollision<CloudType>::PairCollision
(
    const dictionary& dict,
    CloudType& owner
)
:
    CollisionModel<CloudType>(dict, owner, typeName),
    pairModel_
    (
        PairModel<CloudType>::New
        (
            this->coeffDict(),
            this->owner()
        )
    ),
    wallModel_
    (
        WallModel<CloudType>::New
        (
            this->coeffDict(),
            this->owner()
        )
    ),
    il_
    (
        owner.mesh(),
        this->coeffDict().getScalar("maxInteractionDistance"),
        this->coeffDict().getOrDefault
        (
            "writeReferredParticleCloud",
            false
        ),
        this->coeffDict().template getOrDefault<word>("U", "U")
    )
{}
 
 
template<class CloudType>
Foam::PairCollision<CloudType>::PairCollision
(
    const PairCollision<CloudType>& cm
)
:
    CollisionModel<CloudType>(cm),
    pairModel_(nullptr),
    wallModel_(nullptr),
    il_(cm.owner().mesh())
{
    // Need to clone to PairModel and WallModel
    NotImplemented;
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::PairCollision<CloudType>::~PairCollision()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::label Foam::PairCollision<CloudType>::nSubCycles() const
{
    label nSubCycles = 1;
 
    if (pairModel_->controlsTimestep())
    {
        label nPairSubCycles = returnReduce
        (
            pairModel_->nSubCycles(), maxOp<label>()
        );
 
        nSubCycles = max(nSubCycles, nPairSubCycles);
    }
 
    if (wallModel_->controlsTimestep())
    {
        label nWallSubCycles = returnReduce
        (
            wallModel_->nSubCycles(), maxOp<label>()
        );
 
        nSubCycles = max(nSubCycles, nWallSubCycles);
    }
 
    return nSubCycles;
}
 
 
template<class CloudType>
void Foam::PairCollision<CloudType>::collide()
{
    preInteraction();
 
    parcelInteraction();
 
    wallInteraction();
 
    postInteraction();
}
 
 
// ************************************************************************* //