From 9a93918e1dea2ebb1a43afa5ae0977f748e1345c Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Tue, 27 Jan 2009 14:47:49 +0000
Subject: [PATCH] Measuring temperature by the average kinetic energy per
 degree of freedom to allow for a mix of linear, mono-atomic and 6dof
 molecules in a case.

---
 .../mdTools/meanMomentumEnergyAndNMols.H      | 43 +++++++++++--------
 .../molecule/mdTools/temperatureAndPressure.H | 14 +++---
 .../mdTools/temperatureAndPressureVariables.H |  2 +
 .../molecule/molecule/molecule.H              |  2 +
 .../molecule/molecule/moleculeI.H             | 19 +++++++-
 5 files changed, 54 insertions(+), 26 deletions(-)

diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H b/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
index 8688d63669..be97a829c1 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
@@ -52,6 +52,8 @@ scalar singleStepTotalrDotf = 0.0;
 
 label singleStepNMols = molecules.size();
 
+label singleStepDOFs = 0;
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -85,9 +87,11 @@ label singleStepNMols = molecules.size();
     {
         label molId = mol().id();
 
-        scalar molMass(molecules.constProps(molId).mass());
+        const molecule::constantProperties cP(molecules.constProps(molId));
 
-        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+        scalar molMass(cP.mass());
+
+        const diagTensor& molMoI(cP.momentOfInertia());
 
         const vector& molV(mol().v());
 
@@ -112,6 +116,8 @@ label singleStepNMols = molecules.size();
         singleStepTotalPE += mol().potentialEnergy();
 
         singleStepTotalrDotf += tr(mol().rf());
+
+        singleStepDOFs += cP.degreesOfFreedom();
     }
 }
 
@@ -134,46 +140,49 @@ if (Pstream::parRun())
     reduce(singleStepTotalrDotf, sumOp<scalar>());
 
     reduce(singleStepNMols, sumOp<label>());
+
+    reduce(singleStepDOFs, sumOp<label>());
 }
 
 if (singleStepNMols)
 {
-    Info<< "Number of mols in system = "
+    Info<< "Number of molecules in system = "
         << singleStepNMols << nl
         << "Overall number density = "
         << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
         << singleStepTotalMass/meshVolume << nl
-        << "Average linear momentum per mol = "
+        << "Average linear momentum per molecule = "
         << singleStepTotalLinearMomentum/singleStepNMols << ' '
         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
-        << "Average angular momentum per mol = "
+        << "Average angular momentum per molecule = "
         << singleStepTotalAngularMomentum << ' '
         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
         << singleStepMaxVelocityMag << nl
-        << "Average linear KE per mol = "
+        << "Average linear KE per molecule = "
         << singleStepTotalLinearKE/singleStepNMols << nl
-        << "Average angular KE per mol = "
+        << "Average angular KE per molecule = "
         << singleStepTotalAngularKE/singleStepNMols << nl
-        << "Average PE per mol = "
+        << "Average PE per molecule = "
         << singleStepTotalPE/singleStepNMols << nl
-        << "Average TE per mol = "
+        << "Average TE per molecule = "
         <<
         (
             singleStepTotalLinearKE
-            + singleStepTotalAngularKE
-            + singleStepTotalPE
+          + singleStepTotalAngularKE
+          + singleStepTotalPE
         )
         /singleStepNMols
         << endl;
 
-//     Info << singleStepNMols << " "
-// //     << singleStepTotalMomentum/singleStepTotalMass << " "
-//     << singleStepMaxVelocityMag << " "
-//     << singleStepTotalKE/singleStepNMols << " "
-//     << singleStepTotalPE/singleStepNMols << " "
-//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+        // Info << singleStepNMols << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepMaxVelocityMag << " "
+        //     << singleStepTotalKE/singleStepNMols << " "
+        //     << singleStepTotalPE/singleStepNMols << " "
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
 }
 else
 {
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
index 250fa1e2ee..0718750e53 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
@@ -47,17 +47,17 @@ accumulatedTotalrDotfSum += singleStepTotalrDotf;
 
 accumulatedNMols += singleStepNMols;
 
+accumulatedDOFs += singleStepDOFs;
+
 if (runTime.outputTime())
 {
-    // calculate averages
-
     if (accumulatedNMols)
     {
         Info << "calculating averages" << endl;
 
         averageTemperature =
         (
-            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(moleculeCloud::kb * accumulatedDOFs)
             *
             (
                 accumulatedTotalLinearKE + accumulatedTotalAngularKE
@@ -79,8 +79,6 @@ if (runTime.outputTime())
             meshVolume
         );
 
-        // output values
-
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
@@ -91,12 +89,10 @@ if (runTime.outputTime())
     }
     else
     {
-        Info << "Not averaging temperature and pressure: "
+        Info<< "Not averaging temperature and pressure: "
             << "no molecules in system" << endl;
     }
 
-    // reset counters
-
     accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
@@ -110,6 +106,8 @@ if (runTime.outputTime())
     accumulatedTotalrDotfSum = 0.0;
 
     accumulatedNMols = 0;
+
+    accumulatedDOFs = 0;
 }
 
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
index e5e236271f..196299e793 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
@@ -44,6 +44,8 @@ scalar accumulatedTotalrDotfSum = 0.0;
 
 label accumulatedNMols = 0;
 
+label accumulatedDOFs = 0;
+
 scalar averageTemperature = 0.0;
 
 scalar averagePressure = 0.0;
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 7cf42d351b..6734671596 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -141,6 +141,8 @@ public:
 
         inline bool pointMolecule() const;
 
+        inline label degreesOfFreedom() const;
+
         inline scalar mass() const;
 
         inline label nSites() const;
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index c514c336d0..824a8acf80 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -423,7 +423,7 @@ Foam::molecule::constantProperties::momentOfInertia() const
 
 inline bool Foam::molecule::constantProperties::linearMolecule() const
 {
-    return (momentOfInertia_.xx() < 0);
+    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
 }
 
 
@@ -433,6 +433,23 @@ inline bool Foam::molecule::constantProperties::pointMolecule() const
 }
 
 
+inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
+{
+    if (linearMolecule())
+    {
+        return 5;
+    }
+    else if (pointMolecule())
+    {
+        return 3;
+    }
+    else
+    {
+        return 6;
+    }
+}
+
+
 inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;