From 0ac1fb80288cad57bc038e0adc2a8199d372d502 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 2 Feb 2016 19:53:04 +0000 Subject: [PATCH] lagrangian: Move penetration function from KinematicCloud to SprayCloud Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=992 --- .../Templates/KinematicCloud/KinematicCloud.H | 4 - .../KinematicCloud/KinematicCloudI.H | 136 ------------------ .../kinematicCloud/kinematicCloud.H | 8 +- .../clouds/Templates/SprayCloud/SprayCloud.H | 6 + .../clouds/Templates/SprayCloud/SprayCloudI.H | 136 ++++++++++++++++++ 5 files changed, 143 insertions(+), 147 deletions(-) diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H index b25942b5e1..46d12b09bf 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H @@ -400,7 +400,6 @@ public: //- Optional particle forces -// inline const typename parcelType::forceType& forces() const; inline const forceType& forces() const; //- Return the optional particle forces @@ -499,9 +498,6 @@ public: //- Total rotational kinetic energy in the system inline scalar rotationalKineticEnergyOfSystem() const; - //- Penetration for fraction [0-1] of the current total mass - inline scalar penetration(const scalar fraction) const; - //- Mean diameter Dij inline scalar Dij(const label i, const label j) const; diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H index f386091e27..ced4c77f1f 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H @@ -355,142 +355,6 @@ inline Foam::scalar Foam::KinematicCloud::Dmax() const } -template -inline Foam::scalar Foam::KinematicCloud::penetration -( - const scalar fraction -) const -{ - if ((fraction < 0) || (fraction > 1)) - { - FatalErrorInFunction - << "fraction should be in the range 0 < fraction < 1" - << exit(FatalError); - } - - scalar distance = 0.0; - - const label nParcel = this->size(); - globalIndex globalParcels(nParcel); - const label nParcelSum = globalParcels.size(); - - if (nParcelSum == 0) - { - return distance; - } - - // lists of parcels mass and distance from initial injection point - List> procMass(Pstream::nProcs()); - List> procDist(Pstream::nProcs()); - - List& mass = procMass[Pstream::myProcNo()]; - List& dist = procDist[Pstream::myProcNo()]; - - mass.setSize(nParcel); - dist.setSize(nParcel); - - label i = 0; - scalar mSum = 0.0; - forAllConstIter(typename KinematicCloud, *this, iter) - { - const parcelType& p = iter(); - scalar m = p.nParticle()*p.mass(); - scalar d = mag(p.position() - p.position0()); - mSum += m; - - mass[i] = m; - dist[i] = d; - - i++; - } - - // calculate total mass across all processors - reduce(mSum, sumOp()); - Pstream::gatherList(procMass); - Pstream::gatherList(procDist); - - if (Pstream::master()) - { - // flatten the mass lists - List allMass(nParcelSum, 0.0); - SortableList allDist(nParcelSum, 0.0); - for (label procI = 0; procI < Pstream::nProcs(); procI++) - { - SubList - ( - allMass, - globalParcels.localSize(procI), - globalParcels.offset(procI) - ).assign(procMass[procI]); - - // flatten the distance list - SubList - ( - allDist, - globalParcels.localSize(procI), - globalParcels.offset(procI) - ).assign(procDist[procI]); - } - - // sort allDist distances into ascending order - // note: allMass masses are left unsorted - allDist.sort(); - - if (nParcelSum > 1) - { - const scalar mLimit = fraction*mSum; - const labelList& indices = allDist.indices(); - - if (mLimit > (mSum - allMass[indices.last()])) - { - distance = allDist.last(); - } - else - { - // assuming that 'fraction' is generally closer to 1 than 0, - // loop through in reverse distance order - const scalar mThreshold = (1.0 - fraction)*mSum; - scalar mCurrent = 0.0; - label i0 = 0; - - forAllReverse(indices, i) - { - label indI = indices[i]; - - mCurrent += allMass[indI]; - - if (mCurrent > mThreshold) - { - i0 = i; - break; - } - } - - if (i0 == indices.size() - 1) - { - distance = allDist.last(); - } - else - { - // linearly interpolate to determine distance - scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]]; - distance = - allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]); - } - } - } - else - { - distance = allDist.first(); - } - } - - Pstream::scatter(distance); - - return distance; -} - - template inline Foam::cachedRandom& Foam::KinematicCloud::rndGen() { diff --git a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H index 07bcd970b0..b3caba28d2 100644 --- a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,12 +85,6 @@ public: //- Total linear kinetic energy in the system virtual scalar linearKineticEnergyOfSystem() const = 0; - //- Total rotational kinetic energy in the system -// virtual scalar rotationalKineticEnergyOfSystem() const = 0; - - //- Penetration for percentage of the current total mass -// virtual scalar penetration(const scalar& fraction) const = 0; - //- Mean diameter Dij virtual scalar Dij(const label i, const label j) const = 0; diff --git a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H index 68b84d850f..f5f0e887dd 100644 --- a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H +++ b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H @@ -183,6 +183,12 @@ public: inline scalar averageParcelMass() const; + // Check + + //- Penetration for fraction [0-1] of the current total mass + inline scalar penetration(const scalar fraction) const; + + // Sub-models //- Return const-access to the atomization model diff --git a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H index ecabcb71a2..ba5bc78aa7 100644 --- a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H +++ b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H @@ -72,4 +72,140 @@ inline Foam::scalar Foam::SprayCloud::averageParcelMass() const } +template +inline Foam::scalar Foam::SprayCloud::penetration +( + const scalar fraction +) const +{ + if ((fraction < 0) || (fraction > 1)) + { + FatalErrorInFunction + << "fraction should be in the range 0 < fraction < 1" + << exit(FatalError); + } + + scalar distance = 0.0; + + const label nParcel = this->size(); + globalIndex globalParcels(nParcel); + const label nParcelSum = globalParcels.size(); + + if (nParcelSum == 0) + { + return distance; + } + + // lists of parcels mass and distance from initial injection point + List> procMass(Pstream::nProcs()); + List> procDist(Pstream::nProcs()); + + List& mass = procMass[Pstream::myProcNo()]; + List& dist = procDist[Pstream::myProcNo()]; + + mass.setSize(nParcel); + dist.setSize(nParcel); + + label i = 0; + scalar mSum = 0.0; + forAllConstIter(typename SprayCloud, *this, iter) + { + const parcelType& p = iter(); + scalar m = p.nParticle()*p.mass(); + scalar d = mag(p.position() - p.position0()); + mSum += m; + + mass[i] = m; + dist[i] = d; + + i++; + } + + // calculate total mass across all processors + reduce(mSum, sumOp()); + Pstream::gatherList(procMass); + Pstream::gatherList(procDist); + + if (Pstream::master()) + { + // flatten the mass lists + List allMass(nParcelSum, 0.0); + SortableList allDist(nParcelSum, 0.0); + for (label procI = 0; procI < Pstream::nProcs(); procI++) + { + SubList + ( + allMass, + globalParcels.localSize(procI), + globalParcels.offset(procI) + ).assign(procMass[procI]); + + // flatten the distance list + SubList + ( + allDist, + globalParcels.localSize(procI), + globalParcels.offset(procI) + ).assign(procDist[procI]); + } + + // sort allDist distances into ascending order + // note: allMass masses are left unsorted + allDist.sort(); + + if (nParcelSum > 1) + { + const scalar mLimit = fraction*mSum; + const labelList& indices = allDist.indices(); + + if (mLimit > (mSum - allMass[indices.last()])) + { + distance = allDist.last(); + } + else + { + // assuming that 'fraction' is generally closer to 1 than 0, + // loop through in reverse distance order + const scalar mThreshold = (1.0 - fraction)*mSum; + scalar mCurrent = 0.0; + label i0 = 0; + + forAllReverse(indices, i) + { + label indI = indices[i]; + + mCurrent += allMass[indI]; + + if (mCurrent > mThreshold) + { + i0 = i; + break; + } + } + + if (i0 == indices.size() - 1) + { + distance = allDist.last(); + } + else + { + // linearly interpolate to determine distance + scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]]; + distance = + allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]); + } + } + } + else + { + distance = allDist.first(); + } + } + + Pstream::scatter(distance); + + return distance; +} + + // ************************************************************************* //