RemoveParcels.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) 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 "RemoveParcels.H"
#include "fvMesh.H"
#include "faceZone.H"
#include "OFstream.H"
#include "surfaceFields.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template<class CloudType>
void Foam::RemoveParcels<CloudType>::makeLogFile
(
    const word& zoneName,
    const label zoneI,
    const label nFaces,
    const scalar totArea
)
{
    // Create the output file if not already created
    if (log_)
    {
        DebugInfo<< "Creating output file." << endl;
 
        if (Pstream::master())
        {
            // Create directory if does not exist
            mkDir(this->writeTimeDir());
 
            // Open new file at start up
            outputFilePtr_.set
            (
                zoneI,
                new OFstream
                (
                    this->writeTimeDir()/(type() + '_' + zoneName + ".dat")
                )
            );
 
            outputFilePtr_[zoneI]
                << "# Source    : " << type() << nl
                << "# Face zone : " << zoneName << nl
                << "# Faces     : " << nFaces << nl
                << "# Area      : " << totArea << nl
                << "# Time" << tab << "nParcels" << tab << "mass" << endl;
        }
    }
}
 
 
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
template<class CloudType>
void Foam::RemoveParcels<CloudType>::postEvolve
(
    const typename parcelType::trackingData& td
)
{
    Info<< this->modelName() << " output:" << nl;
 
    const fvMesh& mesh = this->owner().mesh();
    const faceZoneMesh& fzm = mesh.faceZones();
 
    forAll(faceZoneIDs_, i)
    {
        const word& zoneName = fzm[faceZoneIDs_[i]].name();
 
        scalar zoneMass = returnReduce(mass_[i], sumOp<scalar>());
        label zoneNParcels = returnReduce(nParcels_[i], sumOp<label>());
 
        Info<< "    faceZone " << zoneName
            << ": removed " << zoneNParcels
            << " parcels with mass " << zoneMass
            << nl;
    }
 
    CloudFunctionObject<CloudType>::postEvolve(td);
}
 
 
template<class CloudType>
void Foam::RemoveParcels<CloudType>::write()
{
    const fvMesh& mesh = this->owner().mesh();
    const Time& time = mesh.time();
 
 
    List<scalar> allZoneMass(faceZoneIDs_.size(), 0.0);
    List<label> allZoneNParcels(faceZoneIDs_.size(), 0);
 
    forAll(faceZoneIDs_, i)
    {
        allZoneMass[i] = returnReduce(mass_[i], sumOp<scalar>());
        allZoneNParcels[i] = returnReduce(nParcels_[i], sumOp<label>());
 
        if (outputFilePtr_.set(i))
        {
            OFstream& os = outputFilePtr_[i];
            os  << time.timeName() << token::TAB
                << allZoneNParcels[i] << token::TAB
                << allZoneMass[i] << endl;
        }
    }
 
    Info<< endl;
 
    if (resetOnWrite_)
    {
        forAll(mass_, i)
        {
            mass_[i] = 0.0;
            nParcels_[i] = 0.0;
        }
    }
 
    this->setModelProperty("mass", allZoneMass);
    this->setModelProperty("nParcels", allZoneNParcels);
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::RemoveParcels<CloudType>::RemoveParcels
(
    const dictionary& dict,
    CloudType& owner,
    const word& modelName
)
:
    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
    faceZoneIDs_(),
    nParcels_(),
    mass_(),
    typeId_(this->coeffDict().template getOrDefault<label>("parcelType", -1)),
    log_(this->coeffDict().getBool("log")),
    resetOnWrite_(this->coeffDict().getBool("resetOnWrite")),
    resetOnStart_(this->coeffDict().getBool("resetOnStart")),
    outputFilePtr_()
{
    const wordList faceZoneNames(this->coeffDict().lookup("faceZones"));
 
    nParcels_.setSize(faceZoneNames.size(), 0);
    mass_.setSize(faceZoneNames.size(), 0.0);
 
    if (!resetOnStart_ && Pstream::master())
    {
        this->getModelProperty("mass", mass_);
        this->getModelProperty("nParcels", nParcels_);
    }
 
    outputFilePtr_.setSize(faceZoneNames.size());
 
    DynamicList<label> zoneIDs;
    const faceZoneMesh& fzm = owner.mesh().faceZones();
    const surfaceScalarField& magSf = owner.mesh().magSf();
    const polyBoundaryMesh& pbm = owner.mesh().boundaryMesh();
 
    forAll(faceZoneNames, i)
    {
        const word& zoneName = faceZoneNames[i];
        label zonei = fzm.findZoneID(zoneName);
 
        if (zonei != -1)
        {
            zoneIDs.append(zonei);
            const faceZone& fz = fzm[zonei];
 
            label nFaces = returnReduce(fz.size(), sumOp<label>());
            Info<< "        " << zoneName << " faces: " << nFaces << nl;
 
            scalar totArea = 0.0;
            for (label facei : fz)
            {
                if (facei < owner.mesh().nInternalFaces())
                {
                    totArea += magSf[facei];
                }
                else
                {
                    label bFacei = facei - owner.mesh().nInternalFaces();
                    label patchi = pbm.patchID()[bFacei];
                    const polyPatch& pp = pbm[patchi];
 
                    if
                    (
                        !magSf.boundaryField()[patchi].coupled()
                     || refCast<const coupledPolyPatch>(pp).owner()
                    )
                    {
                        label localFacei = pp.whichFace(facei);
                        totArea += magSf.boundaryField()[patchi][localFacei];
                    }
                }
            }
            totArea = returnReduce(totArea, sumOp<scalar>());
 
            makeLogFile(zoneName, i, nFaces, totArea);
        }
    }
 
    faceZoneIDs_.transfer(zoneIDs);
}
 
 
template<class CloudType>
Foam::RemoveParcels<CloudType>::RemoveParcels
(
    const RemoveParcels<CloudType>& rpf
)
:
    CloudFunctionObject<CloudType>(rpf),
    faceZoneIDs_(rpf.faceZoneIDs_),
    nParcels_(rpf.nParcels_),
    mass_(rpf.mass_),
    typeId_(rpf.typeId_),
    log_(rpf.log_),
    resetOnWrite_(rpf.resetOnWrite_),
    resetOnStart_(rpf.resetOnStart_),
    outputFilePtr_()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
void Foam::RemoveParcels<CloudType>::postFace
(
    const parcelType& p,
    bool& keepParticle
)
{
    if ((typeId_ >= 0) && (p.typeId() != typeId_))
    {
        // Not processing this particle type
        return;
    }
 
    if
    (
        this->owner().solution().output()
     || this->owner().solution().transient()
    )
    {
        const faceZoneMesh& fzm = this->owner().mesh().faceZones();
 
        forAll(faceZoneIDs_, i)
        {
            const faceZone& fz = fzm[faceZoneIDs_[i]];
            if (fz.found(p.face()))
            {
                nParcels_[i]++;
                mass_[i] += p.mass()*p.nParticle();
                keepParticle = false;
                break;
            }
        }
    }
}
 
 
// ************************************************************************* //