Merge branch 'feature-lagrangian' into 'develop'

ENH: Updated lagrangian libraries from internal development line



See merge request !22

src/lagrangian/coalCombustion/Make/options

@@ -18,6 +18,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -43,6 +44,7 @@ LIB_LIBS = \
     -lincompressibleTransportModels \
     -lregionModels \
     -lsurfaceFilmModels \
+    -ldynamicMesh \
     -ldynamicFvMesh \
     -lfiniteVolume \
     -lmeshTools

src/lagrangian/intermediate/Make/files

+submodels/CloudFunctionObjects/CloudToVTK/vtkTools.C
+
 PARCELS=parcels
 BASEPARCELS=$(PARCELS)/baseClasses
 DERIVEDPARCELS=$(PARCELS)/derived

src/lagrangian/intermediate/Make/options

@@ -14,6 +14,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -35,6 +36,7 @@ LIB_LIBS = \
     -lincompressibleTransportModels \
     -lregionModels \
     -lsurfaceFilmModels \
+    -ldynamicMesh \
     -ldynamicFvMesh \
     -lsampling \
     -lfiniteVolume \

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C

@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,11 +28,7 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cloudSolution::cloudSolution
-(
-    const fvMesh& mesh,
-    const dictionary& dict
-)
+Foam::cloudSolution::cloudSolution(const fvMesh& mesh, const dictionary& dict)
 :
     mesh_(mesh),
     dict_(dict),
@@ -42,6 +38,7 @@ Foam::cloudSolution::cloudSolution
     maxCo_(0.3),
     iter_(1),
     trackTime_(0.0),
+    deltaTMax_(GREAT),
     coupled_(false),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
@@ -52,13 +49,25 @@ Foam::cloudSolution::cloudSolution
     {
         read();
     }
+    else
+    {
+        // see if existing source terms should be reset
+        const dictionary sourceTerms(dict_.subOrEmptyDict("sourceTerms"));
+        sourceTerms.readIfPresent("resetOnStartup", resetSourcesOnStartup_);
+
+        if (resetSourcesOnStartup_)
+        {
+            Info<< "Cloud source terms will be reset" << endl;
+        }
+        else
+        {
+            Info<< "Cloud source terms will be held constant" << endl;
+        }
+    }
 }
 
 
-Foam::cloudSolution::cloudSolution
-(
-    const cloudSolution& cs
-)
+Foam::cloudSolution::cloudSolution(const cloudSolution& cs)
 :
     mesh_(cs.mesh_),
     dict_(cs.dict_),
@@ -68,6 +77,7 @@ Foam::cloudSolution::cloudSolution
     maxCo_(cs.maxCo_),
     iter_(cs.iter_),
     trackTime_(cs.trackTime_),
+    deltaTMax_(cs.deltaTMax_),
     coupled_(cs.coupled_),
     cellValueSourceCorrection_(cs.cellValueSourceCorrection_),
     maxTrackTime_(cs.maxTrackTime_),
@@ -76,10 +86,7 @@ Foam::cloudSolution::cloudSolution
 {}
 
 
-Foam::cloudSolution::cloudSolution
-(
-    const fvMesh& mesh
-)
+Foam::cloudSolution::cloudSolution(const fvMesh& mesh)
 :
     mesh_(mesh),
     dict_(dictionary::null),
@@ -89,6 +96,7 @@ Foam::cloudSolution::cloudSolution
     maxCo_(GREAT),
     iter_(0),
     trackTime_(0.0),
+    deltaTMax_(GREAT),
     coupled_(false),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
@@ -111,6 +119,7 @@ void Foam::cloudSolution::read()
     dict_.lookup("coupled") >> coupled_;
     dict_.lookup("cellValueSourceCorrection") >> cellValueSourceCorrection_;
     dict_.readIfPresent("maxCo", maxCo_);
+    dict_.readIfPresent("deltaTMax", deltaTMax_);
 
     if (steadyState())
     {
@@ -229,4 +238,23 @@ bool Foam::cloudSolution::output() const
 }
 
 
+Foam::scalar Foam::cloudSolution::deltaTMax(const scalar trackTime) const
+{
+    if (transient_)
+    {
+        return min(deltaTMax_, maxCo_*trackTime);
+    }
+    else
+    {
+        return min(deltaTMax_, trackTime);
+    }
+}
+
+
+Foam::scalar Foam::cloudSolution::deltaLMax(const scalar lRef) const
+{
+    return maxCo_*lRef;
+}
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,9 @@ class cloudSolution
         //- Particle track time
         scalar trackTime_;
 
+        //- Maximum integration time step (optional)
+        scalar deltaTMax_;
+
 
         // Run-time options
 
@@ -171,6 +174,9 @@ public:
             //- Return the particle track time
             inline scalar trackTime() const;
 
+            //- Return the maximum integation time step
+            inline scalar deltaTMax() const;
+
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
@@ -207,6 +213,12 @@ public:
 
         //- Returns true if writing this step
         bool output() const;
+
+        //- Return the maximum integration time
+        scalar deltaTMax(const scalar trackTime) const;
+
+        //- Return the maximum integration length
+        scalar deltaLMax(const scalar lRef) const;
 };
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,6 +101,12 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
+inline Foam::scalar Foam::cloudSolution::deltaTMax() const
+{
+    return deltaTMax_;
+}
+
+
 inline Foam::Switch& Foam::cloudSolution::coupled()
 {
     return coupled_;

src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,25 +80,25 @@ Foam::MPPICCloud<CloudType>::MPPICCloud
     dampingModel_(NULL),
     isotropyModel_(NULL)
 {
-    if (this->solution().steadyState())
-    {
-        FatalErrorIn
-        (
-            "Foam::MPPICCloud<CloudType>::MPPICCloud"
-            "("
-                "const word&, "
-                "const volScalarField&, "
-                "const volVectorField&, "
-                "const volScalarField&, "
-                "const dimensionedVector&, "
-                "bool"
-            ")"
-        )   << "MPPIC modelling not available for steady state calculations"
-            << exit(FatalError);
-    }
-
     if (this->solution().active())
     {
+        if (this->solution().steadyState())
+        {
+            FatalErrorIn
+            (
+                "Foam::MPPICCloud<CloudType>::MPPICCloud"
+                "("
+                    "const word&, "
+                    "const volScalarField&, "
+                    "const volVectorField&, "
+                    "const volScalarField&, "
+                    "const dimensionedVector&, "
+                    "bool"
+                ")"
+            )   << "MPPIC modelling not available for steady state calculations"
+                << exit(FatalError);
+        }
+
         setModels();
 
         if (readFields)

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,10 +51,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    if (this->solution().coupled())
-    {
-        this->subModelProperties().lookup("radiation") >> radiation_;
-    }
+    this->subModelProperties().lookup("radiation") >> radiation_;
 
     if (radiation_)
     {

src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C

@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,8 @@ bool Foam::CollidingParcel<ParcelType>::move
     typename TrackData::cloudType::parcelType& p =
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
+    td.keepParticle = true;
+
     switch (td.part())
     {
         case TrackData::tpVelocityHalfStep:

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "forceSuSp.H"
 #include "IntegrationScheme.H"
 #include "meshTools.H"
+#include "cloudSolution.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -271,15 +272,11 @@ bool Foam::KinematicParcel<ParcelType>::move
 
     const polyMesh& mesh = td.cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
+    const cloudSolution& solution = td.cloud().solution();
     const scalarField& cellLengthScale = td.cloud().cellLengthScale();
-    const scalar maxCo = td.cloud().solution().maxCo();
 
     scalar tEnd = (1.0 - p.stepFraction())*trackTime;
-    scalar dtMax = trackTime;
-    if (td.cloud().solution().transient())
-    {
-        dtMax *= maxCo;
-    }
+    scalar dtMax = solution.deltaTMax(trackTime);
 
     bool tracking = true;
     label nTrackingStalled = 0;
@@ -301,7 +298,8 @@ bool Foam::KinematicParcel<ParcelType>::move
         if (p.active() && tracking && (magU > ROOTVSMALL))
         {
             const scalar d = dt*magU;
-            const scalar dCorr = min(d, maxCo*cellLengthScale[cellI]);
+            const scalar deltaLMax = solution.deltaLMax(cellLengthScale[cellI]);
+            const scalar dCorr = min(d, deltaLMax);
             dt *=
                 dCorr/d
                *p.trackToFace(p.position() + dCorr*U_/magU, td);
@@ -309,7 +307,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         tEnd -= dt;
 
-        scalar newStepFraction = 1.0 - tEnd/trackTime;
+        const scalar newStepFraction = 1.0 - tEnd/trackTime;
 
         if (tracking)
         {
@@ -335,7 +333,7 @@ bool Foam::KinematicParcel<ParcelType>::move
         p.stepFraction() = newStepFraction;
 
         bool calcParcel = true;
-        if (!tracking && td.cloud().solution().steadyState())
+        if (!tracking && solution.steadyState())
         {
             calcParcel = false;
         }
@@ -346,7 +344,7 @@ bool Foam::KinematicParcel<ParcelType>::move
             // Update cell based properties
             p.setCellValues(td, dt, cellI);
 
-            if (td.cloud().solution().cellValueSourceCorrection())
+            if (solution.cellValueSourceCorrection())
             {
                 p.cellValueSourceCorrection(td, dt, cellI);
             }
@@ -466,6 +464,8 @@ void Foam::KinematicParcel<ParcelType>::hitPatch
 )
 {
     td.keepParticle = false;
+
+    td.cloud().patchInteraction().addToEscapedParcels(nParticle_*mass());
 }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -121,7 +121,13 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 
     YMix[GAS] = massGas/massNew;
     YMix[LIQ] = massLiquid/massNew;
-    YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
+    YMix[SLD] = massSolid/massNew;
+
+    scalar Ytotal = sum(YMix);
+
+    YMix[GAS] /= Ytotal;
+    YMix[LIQ] /= Ytotal;
+    YMix[SLD] /= Ytotal;
 
     return massNew;
 }

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -80,6 +80,9 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
     this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
 
     const scalar CpMean = td.CpInterp().psi()[cellI];
+
+    tetIndices tetIs = this->currentTetIndices();
+    Tc_ = td.TInterp().interpolate(this->position(), tetIs);
     Tc_ += td.cloud().hsTrans()[cellI]/(CpMean*this->massCell(cellI));
 
     if (Tc_ < td.cloud().constProps().TMin())

src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "CloudToVTK.H"
 #include "FacePostProcessing.H"
 #include "ParticleCollector.H"
 #include "ParticleErosion.H"
@@ -42,6 +43,7 @@ License
                                                                               \
     makeCloudFunctionObject(CloudType);                                       \
                                                                               \
+    makeCloudFunctionObjectType(CloudToVTK, CloudType);                       \
     makeCloudFunctionObjectType(FacePostProcessing, CloudType);               \
     makeCloudFunctionObjectType(ParticleCollector, CloudType);                \
     makeCloudFunctionObjectType(ParticleErosion, CloudType);                  \

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudToVTK/CloudToVTK.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 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 "CloudToVTK.H"
+#include "vtkTools.H"
+#include "floatScalar.H"
+
+// * * * * * * * * * * * * * Protectd Member Functions * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::CloudToVTK<CloudType>::writeData
+(
+    std::ostream& vtkOs,
+    const bool binary,
+    const List<floatScalar>& data
+) const
+{
+    const label procI = Pstream::myProcNo();
+
+    List<List<floatScalar> > allProcData(Pstream::nProcs());
+    allProcData[procI] = data;
+    Pstream::gatherList(allProcData);
+    List<floatScalar> allData =
+        ListListOps::combine<List<floatScalar> >
+        (
+            allProcData,
+            accessOp<List<floatScalar> >()
+        );
+
+    vtkTools::write(vtkOs, binary, allData);
+}
+
+
+template<class CloudType>
+template<class Type>
+void Foam::CloudToVTK<CloudType>::writeFieldData
+(
+    std::ostream& vtkOs,
+    const bool binary,
+    const List<floatScalar>& data,
+    const word& title,
+    const label nParcels
+) const
+{
+    vtkOs
+        << title << ' ' << pTraits<Type>::nComponents << ' '
+        << nParcels << " float" << std::endl;
+    writeData(vtkOs, binary, data);
+}
+
+
+template<class CloudType>
+void Foam::CloudToVTK<CloudType>::write()
+{
+    label nParcels = this->owner().size();
+    DynamicList<floatScalar> position(3*nParcels);
+    DynamicList<floatScalar> U(3*nParcels);
+    DynamicList<floatScalar> d(nParcels);
+    DynamicList<floatScalar> age(nParcels);
+    DynamicList<floatScalar> rho(nParcels);
+
+    forAllConstIter(typename CloudType, this->owner(), iter)
+    {
+        vtkTools::insert(iter().position(), position);
+        vtkTools::insert(iter().U(), U);
+        vtkTools::insert(iter().d(), d);
+        vtkTools::insert(iter().age(), age);
+        vtkTools::insert(iter().rho(), rho);
+    }
+
+    reduce(nParcels, sumOp<label>());
+
+
+binary_ = false;
+    if (Pstream::master())
+    {
+        // Create directory if does not exist
+        mkDir(this->outputTimeDir());
+
+        // Open new file at start up
+
+        const fileName fName = this->outputTimeDir()/(type() + ".vtk");
+        this->setModelProperty("file", fName);
+
+        OFstream os(fName, binary_ ? IOstream::BINARY : IOstream::ASCII);
+        std::ostream& vtkOs = os.stdStream();
+
+
+        vtkTools::writeHeader(vtkOs, binary_, this->modelName().c_str());
+        vtkOs
+            << "DATASET POLYDATA" << std::endl
+            << "POINTS " << nParcels << " float" << std::endl;
+
+        writeData(vtkOs, binary_, position);
+
+        vtkOs
+            << "POINT_DATA " << nParcels << std::endl
+            << "FIELD attributes " << 4
+            << std::endl;
+
+        writeFieldData<vector>(vtkOs, binary_, U, "U", nParcels);
+        writeFieldData<scalar>(vtkOs, binary_, d, "d", nParcels);
+        writeFieldData<scalar>(vtkOs, binary_, age, "age", nParcels);
+        writeFieldData<scalar>(vtkOs, binary_, rho, "rho", nParcels);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::CloudToVTK<CloudType>::CloudToVTK
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
+    binary_(dict.lookupOrDefault<bool>("binary", true))
+{}
+
+
+template<class CloudType>
+Foam::CloudToVTK<CloudType>::CloudToVTK
+(
+    const CloudToVTK<CloudType>& c
+)
+:
+    CloudFunctionObject<CloudType>(c),
+    binary_(c.binary_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::CloudToVTK<CloudType>::~CloudToVTK()
+{}
+
+
+// ************************************************************************* //