From ec98d7b49b01641eef96c618b1db289adb5f2c6d Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 28 Jun 2012 13:38:06 +0100 Subject: [PATCH] regionSizeDistribution: New functionObject to calculate the size distribution of regions defined by a field e.g. the size distribution of droplets in a VoF spray calculation --- .../functionObjects/field/Make/files | 3 + .../regionSizeDistribution.C | 588 ++++++++++++++++++ .../regionSizeDistribution.H | 161 +++++ .../regionSizeDistributionFunctionObject.C | 46 ++ .../regionSizeDistributionFunctionObject.H | 54 ++ .../incompressible/LES/Make/files | 1 - 6 files changed, 852 insertions(+), 1 deletion(-) create mode 100644 src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C create mode 100644 src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H create mode 100644 src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.C create mode 100644 src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.H diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files index e351784bdc..0abf314f96 100644 --- a/src/postProcessing/functionObjects/field/Make/files +++ b/src/postProcessing/functionObjects/field/Make/files @@ -41,4 +41,7 @@ wallBoundedStreamLine/wallBoundedParticle.C surfaceInterpolateFields/surfaceInterpolateFields.C surfaceInterpolateFields/surfaceInterpolateFieldsFunctionObject.C +regionSizeDistribution/regionSizeDistribution.C +regionSizeDistribution/regionSizeDistributionFunctionObject.C + LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C new file mode 100644 index 0000000000..b3aaac0dbf --- /dev/null +++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C @@ -0,0 +1,588 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 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 . + +\*---------------------------------------------------------------------------*/ + +#include "regionSizeDistribution.H" +#include "volFields.H" +#include "regionSplit.H" +#include "fvcVolumeIntegrate.H" +#include "Histogram.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(regionSizeDistribution, 0); + + //- plus op for FixedList + template + class ListPlusEqOp + { + public: + void operator() + ( + FixedList& x, + const FixedList& y + ) const + { + forAll(x, i) + { + x[i] += y[i]; + } + } + }; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::regionSizeDistribution::regionSizeDistribution +( + const word& name, + const objectRegistry& obr, + const dictionary& dict, + const bool loadFromFiles +) +: + name_(name), + obr_(obr), + active_(true), + alphaName_(dict.lookup("field")), + patchNames_(dict.lookup("patches")) +{ + // Check if the available mesh is an fvMesh, otherwise deactivate + if (!isA(obr_)) + { + active_ = false; + WarningIn + ( + "regionSizeDistribution::regionSizeDistribution" + "(const objectRegistry&, const dictionary&)" + ) << "No fvMesh available, deactivating." << nl + << endl; + } + + read(dict); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::regionSizeDistribution::~regionSizeDistribution() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::regionSizeDistribution::read(const dictionary& dict) +{ + if (active_) + { + dict.lookup("field") >> alphaName_; + dict.lookup("patches") >> patchNames_; + dict.lookup("threshold") >> threshold_; + dict.lookup("volFraction") >> volFraction_; + dict.lookup("nBins") >> nBins_; + + word format(dict.lookup("setFormat")); + formatterPtr_ = writer::New(format); + } +} + + +void Foam::regionSizeDistribution::execute() +{ + // Do nothing - only valid on write +} + + +void Foam::regionSizeDistribution::end() +{ + // Do nothing - only valid on write +} + + +void Foam::regionSizeDistribution::write() +{ + if (active_) + { + const fvMesh& mesh = refCast(obr_); + + autoPtr alphaPtr; + if (obr_.foundObject(alphaName_)) + { + Info<< "Looking up field " << alphaName_ << endl; + } + else + { + Info<< "Reading field " << alphaName_ << endl; + alphaPtr.reset + ( + new volScalarField + ( + IOobject + ( + alphaName_, + mesh.time().timeName(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh + ) + ); + } + + + const volScalarField& alpha = + ( + alphaPtr.valid() + ? alphaPtr() + : obr_.lookupObject(alphaName_) + ); + + Info<< "Volume of alpha = " + << fvc::domainIntegrate(alpha).value() + << endl; + + const scalar meshVol = gSum(mesh.V()); + Info<< "Mesh volume = " << meshVol << endl; + Info<< "Background region volume limit = " << volFraction_*meshVol + << endl; + + + // Determine blocked faces + boolList blockedFace(mesh.nFaces(), false); + label nBlocked = 0; + + { + for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + { + scalar ownVal = alpha[mesh.faceOwner()[faceI]]; + scalar neiVal = alpha[mesh.faceNeighbour()[faceI]]; + + if + ( + (ownVal < threshold_ && neiVal > threshold_) + || (ownVal > threshold_ && neiVal < threshold_) + ) + { + blockedFace[faceI] = true; + nBlocked++; + } + } + + // Block coupled faces + forAll(alpha.boundaryField(), patchI) + { + const fvPatchScalarField& fvp = alpha.boundaryField()[patchI]; + if (fvp.coupled()) + { + tmp townFld(fvp.patchInternalField()); + const scalarField& ownFld = townFld(); + tmp tnbrFld(fvp.patchNeighbourField()); + const scalarField& nbrFld = tnbrFld(); + + label start = fvp.patch().patch().start(); + + forAll(ownFld, i) + { + scalar ownVal = ownFld[i]; + scalar neiVal = nbrFld[i]; + + if + ( + (ownVal < threshold_ && neiVal > threshold_) + || (ownVal > threshold_ && neiVal < threshold_) + ) + { + blockedFace[start+i] = true; + nBlocked++; + } + } + } + } + } + + + regionSplit regions(mesh, blockedFace); + + Info<< "Determined " << regions.nRegions() << " disconnected regions" + << endl; + + + if (debug) + { + volScalarField region + ( + IOobject + ( + "region", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimless, 0) + ); + Info<< "Dumping region as volScalarField to " << region.name() + << endl; + + forAll(regions, cellI) + { + region[cellI] = regions[cellI]; + } + region.correctBoundaryConditions(); + region.write(); + } + + + // Sum all regions + Map > regionVolume(regions.nRegions()/Pstream::nProcs()); + forAll(alpha, cellI) + { + scalar cellVol = mesh.V()[cellI]; + scalar alphaVol = alpha[cellI]*cellVol; + + label regionI = regions[cellI]; + + Map >::iterator fnd = regionVolume.find(regionI); + if (fnd == regionVolume.end()) + { + regionVolume.insert + ( + regionI, + Pair(cellVol, alphaVol) + ); + } + else + { + fnd().first() += cellVol; + fnd().second() += alphaVol; + } + } + Pstream::mapCombineGather(regionVolume, ListPlusEqOp()); + Pstream::mapCombineScatter(regionVolume); + + + if (debug) + { + Info<< token::TAB << "Region" + << token::TAB << "Volume(mesh)" + << token::TAB << "Volume(" << alpha.name() << "):" + << endl; + scalar meshSumVol = 0.0; + scalar alphaSumVol = 0.0; + + forAllConstIter(Map >, regionVolume, iter) + { + Info<< token::TAB << iter.key() + << token::TAB << iter().first() + << token::TAB << iter().second() << endl; + + meshSumVol += iter().first(); + alphaSumVol += iter().second(); + } + Info<< token::TAB << "Total:" + << token::TAB << meshSumVol + << token::TAB << alphaSumVol << endl; + Info<< endl; + } + + + + // Mark all regions starting at patches + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // Count number of patch faces (just for initial sizing) + label nPatchFaces = 0; + forAll(patchNames_, i) + { + const word& pName = patchNames_[i]; + label patchI = mesh.boundaryMesh().findPatchID(pName); + if (patchI == -1) + { + WarningIn("regionSizeDistribution::write()") + << "Cannot find patch " << pName << ". Valid patches are " + << mesh.boundaryMesh().names() + << endl; + } + else + { + nPatchFaces += mesh.boundaryMesh()[patchI].size(); + } + } + + Map