fv::zeroDimensionalMassSource: Support use of multiple sources

The computation of the current and total accumulated mass is now shared
across all zeroDimensionalMassSource models. So, multiple models can be
used simultaneously.

src/fvModels/zeroDimensional/zeroDimensionalMassSource/zeroDimensionalMassSource.C

@@ -39,49 +39,143 @@ namespace fv
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::fv::zeroDimensionalMassSource::readCoeffs()
+Foam::tmp<Foam::volScalarField>
+Foam::fv::zeroDimensionalMassSource::calcM0D() const
 {
-    const basicThermo& thermo =
+    tmp<volScalarField> tm =
-        mesh().lookupObject<basicThermo>
+        volScalarField::New
         (
-            IOobject::groupName(physicalProperties::typeName, phaseName_)
+            typedName("m0D"),
+            mesh(),
+            dimensionedScalar(dimMass, 0)
         );
 
-    rho0Ptr_.reset
+    HashTable<const basicThermo*> thermos(mesh().lookupClass<basicThermo>());
-    (
+
-        new volScalarField::Internal
+    forAllConstIter(HashTable<const basicThermo*>, thermos, thermoIter)
-        (
+    {
-            IOobject
+        const basicThermo& thermo = *thermoIter();
-            (
+
-                typedName("rho0"),
+        tmp<volScalarField> tRho = thermo.rho();
-                mesh().time().timeName(),
+        const volScalarField& rho = tRho();
-                mesh(),
+
-                IOobject::READ_IF_PRESENT,
+        const word phaseName = thermo.phaseName();
-                IOobject::AUTO_WRITE
+
-            ),
+        if (thermo.phaseName() != word::null)
-            thermo.rho()().internalField()
+        {
-        )
+            const volScalarField alpha =
-    );
+                mesh().lookupObject<volScalarField>
+                (
+                    IOobject::groupName("alpha", phaseName_)
+                );
+
+            tm.ref().ref() += alpha()*rho()*mesh().V();
+        }
+        else
+        {
+            tm.ref().ref() += rho()*mesh().V();
+        }
+    }
+
+    return tm;
 }
 
 
-Foam::scalar Foam::fv::zeroDimensionalMassSource::massFlowRate() const
+Foam::volScalarField& Foam::fv::zeroDimensionalMassSource::initM0D() const
 {
-    const scalar t = mesh().time().userTimeValue();
+    if (!mesh().foundObject<volScalarField>(typedName("m0D")))
-    const scalar t0 = mesh().time().beginTime().value();
+    {
+        volScalarField* mPtr =
+            new volScalarField
+            (
+                calcM0D()
+            );
 
-    const scalar mDot = massFlowRate_->value(t);
+        mPtr->store();
-    const scalar sumMDot = massFlowRate_->integral(t0, t);
+    }
 
-    const basicThermo& thermo =
+    return mesh().lookupObjectRef<volScalarField>(typedName("m0D"));
-        mesh().lookupObject<basicThermo>
+}
+
+
+const Foam::volScalarField& Foam::fv::zeroDimensionalMassSource::m() const
+{
+    // If not registered, then read or create the mass field
+    if (!mesh().foundObject<volScalarField>(typedName("m")))
+    {
+        typeIOobject<volScalarField> mIo
         (
-            IOobject::groupName(physicalProperties::typeName, phaseName_)
+            typedName("m"),
+            mesh().time().timeName(),
+            mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
         );
 
-    return mDot*thermo.rho()()[0]/(rho0Ptr_()[0] + sumMDot/mesh().V()[0]);
+        volScalarField* mPtr =
+            new volScalarField
+            (
+                mIo,
+                mesh(),
+                dimensionedScalar(dimMass, 0)
+            );
+
+        mPtr->store();
+
+        if (!mIo.headerOk())
+        {
+            *mPtr = m0D_;
+        }
+
+        volScalarField* factorPtr =
+            new volScalarField
+            (
+                IOobject
+                (
+                    typedName("factor"),
+                    mesh().time().timeName(),
+                    mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                *mPtr/m0D_
+            );
+
+        factorPtr->store();
+    }
+
+    volScalarField& m =
+        mesh().lookupObjectRef<volScalarField>(typedName("m"));
+
+    volScalarField& factor =
+        mesh().lookupObjectRef<volScalarField>(typedName("factor"));
+
+    // Update the mass if changes are available
+    if (mesh().foundObject<volScalarField>(typedName("deltaM")))
+    {
+        volScalarField& deltaM =
+            mesh().lookupObjectRef<volScalarField>(typedName("deltaM"));
+
+        m = m.oldTime() + deltaM;
+
+        factor = m/m0D_;
+
+        deltaM.checkOut();
+    }
+
+    return m;
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::fv::zeroDimensionalMassSource::massFlowRate() const
+{
+    const scalar mDot = massFlowRate_->value(mesh().time().userTimeValue());
+
+    return mDot*m0D_[0]/m()[0];
 }
 
 
@@ -96,7 +190,7 @@ Foam::fv::zeroDimensionalMassSource::zeroDimensionalMassSource
 )
 :
     massSource(name, modelType, mesh, dict, true),
-    rho0Ptr_()
+    m0D_(initM0D())
 {
     if (mesh.nGeometricD() != 0)
     {
@@ -105,25 +199,41 @@ Foam::fv::zeroDimensionalMassSource::zeroDimensionalMassSource
             << mesh.nGeometricD() << "-dimensional mesh"
             << exit(FatalIOError);
     }
-
-    readCoeffs();
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::fv::zeroDimensionalMassSource::read(const dictionary& dict)
+void Foam::fv::zeroDimensionalMassSource::correct()
 {
-    if (fvModel::read(dict))
+    // Correct the zero-dimensional mass
+    m0D_ = calcM0D();
+
+    // Create the mass change
+    if (!mesh().foundObject<volScalarField>(typedName("deltaM")))
     {
-        massSource::readCoeffs();
+        volScalarField* dMPtr =
-        readCoeffs();
+            new volScalarField
-        return true;
+            (
-    }
+                IOobject
-    else
+                (
-    {
+                    typedName("deltaM"),
-        return false;
+                    mesh().time().timeName(),
+                    mesh()
+                ),
+                mesh(),
+                dimensionedScalar(dimMass, 0)
+            );
+
+        dMPtr->store();
     }
+
+    volScalarField& deltaM =
+        mesh().lookupObjectRef<volScalarField>(typedName("deltaM"));
+
+    const scalar mDot = massFlowRate_->value(mesh().time().userTimeValue());
+
+    deltaM += mesh().time().deltaT()*dimensionedScalar(dimMass/dimTime, mDot);
 }
 
 

src/fvModels/zeroDimensional/zeroDimensionalMassSource/zeroDimensionalMassSource.H

@@ -30,6 +30,19 @@ Description
     for the mass that exits the domain due to expansion in space, so that the
     model correctly applies a total mass flow rate.
 
+    This model will write out additional fields, zeroDimensionalMassSource:m
+    and zeroDimensionalMassSource:factor. The zeroDimensionalMassSource:m field
+    is the total accumulated mass; the sum of the starting mass, plus all mass
+    added by models of this type. This may differ from the actual current mass
+    (= rho*V) within the zero-dimensional cell if pressure or density
+    constraints are being used. The zeroDimensionalMassSource:factor field is
+    the ratio between the current mass or volume and the total accumulated mass
+    or volume.
+
+    If the mass flow rate is positive then user-supplied fixed property values
+    are introduced to the field equations. If the mass flow rate is negative
+    then properties are removed at their current value.
+
 Usage
     Example usage:
     \verbatim
@@ -49,10 +62,10 @@ Usage
     }
     \endverbatim
 
-    Values should be provided for all solved for fields. Warnings will be
+    If the mass flow rate is positive then values should be provided for all
-    issued if values are not provided for fields for which transport equations
+    solved for fields. Warnings will be issued if values are not provided for
-    are solved. Warnings will also be issued if values are provided for fields
+    fields for which transport equations are solved. Warnings will also be
-    which are not solved for.
+    issued if values are provided for fields which are not solved for.
 
 SourceFiles
     zeroDimensionalMassSource.C
@@ -79,19 +92,30 @@ class zeroDimensionalMassSource
 :
     public massSource
 {
+private:
+
+    // Private Data
+
+        //- Reference to the zero-dimensional mass
+        volScalarField& m0D_;
+
+
+    // Private Member Functions
+
+        //- Calculate and return the zero-dimensional mass
+        tmp<volScalarField> calcM0D() const;
+
+        //- Initialise and return a reference to the zero-dimensional mass
+        volScalarField& initM0D() const;
+
+        //- Get a reference to the mass
+        const volScalarField& m() const;
+
+
 protected:
 
-    // Protected Data
-
-        //- Cached initial density field
-        autoPtr<volScalarField::Internal> rho0Ptr_;
-
-
     // Protected Member Functions
 
-        //- Non-virtual read
-        void readCoeffs();
-
         //- Return the mass flow rate
         virtual scalar massFlowRate() const;
 
@@ -116,10 +140,8 @@ public:
 
     // Member Functions
 
-        // IO
+        //- Update the model
-
+        virtual void correct();
-            //- Read source dictionary
-            virtual bool read(const dictionary& dict);
 };