ENH: Running and improvements.

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

@@ -49,6 +49,8 @@ int main(int argc, char *argv[])
 
         molecules.evolve();
 
+        runTime.write();
+
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -26,8 +26,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "md.H"
 #include "fvCFD.H"
+#include "polyatomicCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,7 +64,9 @@ int main(int argc, char *argv[])
 
     potential pot(mesh, mdInitialiseDict, idListDict);
 
-    moleculeCloud molecules(mesh, pot, mdInitialiseDict);
+    polyatomicCloud molecules(mesh, pot, mdInitialiseDict);
+
+    Info<< "Cloud is called " << molecules.name() << endl;
 
     label totalMolecules = molecules.size();
 
@@ -81,7 +83,7 @@ int main(int argc, char *argv[])
     if (!mesh.write())
     {
         FatalErrorIn(args.executable())
-            << "Failed writing moleculeCloud."
+            << "Failed writing polyatomicCloud."
             << nl << exit(FatalError);
     }
 

src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C

@@ -71,6 +71,49 @@ Foam::constPropSite::~constPropSite()
 {}
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, constPropSite& cPS)
+{
+    is  >> cPS.siteReferencePosition_
+        >> cPS.siteMass_
+        >> cPS.siteCharge_
+        >> cPS.siteId_
+        >> cPS.name_
+        >> cPS.pairPotentialSite_
+        >> cPS.electrostaticSite_;
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>"
+        "(Foam::Istream&, Foam::constPropSite&)"
+    );
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const constPropSite& cPS)
+{
+
+    os  << token::SPACE << cPS.siteReferencePosition()
+        << token::SPACE << cPS.siteMass()
+        << token::SPACE << cPS.siteCharge()
+        << token::SPACE << cPS.siteId()
+        << token::SPACE << cPS.name()
+        << token::SPACE << cPS.pairPotentialSite()
+        << token::SPACE << cPS.electrostaticSite();
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::constPropSite&)"
+    );
+
+    return os;
+}
+
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H

@@ -35,14 +35,23 @@ SourceFiles
 #define constPropSite_H
 
 #include "vector.H"
-#include "IFstream.H"
+#include "IOstreams.H"
-#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class constPropSite;
+Istream& operator>>(Istream&, constPropSite&);
+Ostream& operator<<(Ostream&, const constPropSite&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class constPropSite Declaration
 \*---------------------------------------------------------------------------*/
@@ -144,6 +153,18 @@ public:
 
             //-
             inline bool& electrostaticSite();
+
+
+    // Member Operators
+
+        inline bool operator==(const constPropSite&) const;
+        inline bool operator!=(const constPropSite&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, constPropSite&);
+        friend Ostream& operator<<(Ostream&, const constPropSite&);
 };
 
 

src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H

@@ -114,4 +114,31 @@ inline bool& Foam::constPropSite::electrostaticSite()
 }
 
 
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+bool Foam::constPropSite::operator==
+(
+    const constPropSite& rhs
+) const
+{
+    return
+       siteReferencePosition_ == rhs.siteReferencePosition_
+    && siteMass_ == rhs.siteMass_
+    && siteCharge_ == rhs.siteCharge_
+    && siteId_ == rhs.siteId_
+    && name_ == rhs.name_
+    && pairPotentialSite_ == rhs.pairPotentialSite_
+    && electrostaticSite_ == rhs.electrostaticSite_;
+}
+
+
+bool Foam::constPropSite::operator!=
+(
+    const constPropSite& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H

@@ -57,11 +57,12 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
     List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
 
+    sites_.setSize(siteReferencePositions.size());
+
     if
     (
-        ( siteIdNames.size() != sites_.size() )
+        (siteIdNames.size() != sites_.size())
-     || ( siteReferencePositions.size() != sites_.size() )
+     || (siteCharges.size() != sites_.size())
-     || ( siteCharges.size() != sites_.size() )
     )
     {
             FatalErrorIn
@@ -72,12 +73,13 @@ inline Foam::polyatomic::constantProperties::constantProperties
                 ")"
             )
             << "Sizes of site id, charge and "
-            << "referencePositions are not the same: " << sites_.size()
+            << "referencePositions are not the same: " << nl
-            << nl << abort(FatalError);
+            << siteIdNames  << nl
+            << siteReferencePositions << nl
+            << siteCharges << nl
+            << abort(FatalError);
     }
 
-    sites_.setSize(siteIdNames.size());
-
     electrostaticSites_.setSize(sites_.size(), -1);
     pairPotSites_.setSize(sites_.size(), -1);
 
@@ -107,7 +109,7 @@ inline Foam::polyatomic::constantProperties::constantProperties
             electrostaticSites_[electrostaticI++] = sI;
         }
 
-        if (sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
+        if (!sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
         {
             WarningIn
             (

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C

@@ -40,7 +40,7 @@ namespace Foam
 
 void Foam::polyatomicCloud::buildConstProps()
 {
-    Info<< nl << "Reading polyatomicProperties dictionary." << endl;
+    Info<< nl << "Reading mdProperties dictionary." << endl;
 
     const List<word>& idList(pot_.idList());
 
@@ -48,11 +48,11 @@ void Foam::polyatomicCloud::buildConstProps()
 
     const List<word>& siteIdList(pot_.siteIdList());
 
-    IOdictionary polyatomicPropertiesDict
+    IOdictionary mdPropertiesDict
     (
         IOobject
         (
-            "polyatomicProperties",
+            "mdProperties",
             mesh_.time().constant(),
             mesh_,
             IOobject::MUST_READ_IF_MODIFIED,
@@ -65,7 +65,7 @@ void Foam::polyatomicCloud::buildConstProps()
     {
         const word& id = idList[i];
 
-        const dictionary& molDict(polyatomicPropertiesDict.subDict(id));
+        const dictionary& molDict(mdPropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -1183,6 +1183,8 @@ void Foam::polyatomicCloud::evolve()
 
     polyatomic::trackingData td3(*this, 3);
     Cloud<polyatomic>::move(td3, mesh_.time().deltaTValue());
+
+    info();
 }
 
 
@@ -1237,6 +1239,144 @@ void Foam::polyatomicCloud::applyConstraintsAndThermostats
 }
 
 
+void Foam::polyatomicCloud::info() const
+{
+    // Calculates and prints the mean momentum and energy in the system
+    // and the number of molecules.
+
+    vector totalLinearMomentum(vector::zero);
+
+    vector totalAngularMomentum(vector::zero);
+
+    scalar maxVelocityMag = 0.0;
+
+    scalar totalMass = 0.0;
+
+    scalar totalLinearKE = 0.0;
+
+    scalar totalAngularKE = 0.0;
+
+    scalar totalPE = 0.0;
+
+    scalar totalrDotf = 0.0;
+
+    //vector CentreOfMass(vector::zero);
+
+    label nMols = this->size();
+
+    label dofs = 0;
+
+    {
+        forAllConstIter(polyatomicCloud, *this, mol)
+        {
+            const label molId = mol().id();
+
+            scalar molMass(this->constProps(molId).mass());
+
+            totalMass += molMass;
+
+            //CentreOfMass += mol().position()*molMass;
+        }
+
+        // if (nMols)
+        // {
+        //     CentreOfMass /= totalMass;
+        // }
+
+        forAllConstIter(polyatomicCloud, *this, mol)
+        {
+            const label molId = mol().id();
+
+            const polyatomic::constantProperties cP(this->constProps(molId));
+
+            scalar molMass(cP.mass());
+
+            const diagTensor& molMoI(cP.momentOfInertia());
+
+            const vector& molV(mol().v());
+
+            const vector& molOmega(inv(molMoI) & mol().pi());
+
+            vector molPiGlobal = mol().Q() & mol().pi();
+
+            totalLinearMomentum += molV * molMass;
+
+            totalAngularMomentum += molPiGlobal;
+            //+((mol().position() - CentreOfMass) ^ (molV * molMass));
+
+            if (mag(molV) > maxVelocityMag)
+            {
+                maxVelocityMag = mag(molV);
+            }
+
+            totalLinearKE += 0.5*molMass*magSqr(molV);
+
+            totalAngularKE += 0.5*(molOmega & molMoI & molOmega);
+
+            totalPE += mol().potentialEnergy();
+
+            totalrDotf += tr(mol().rf());
+
+            dofs += cP.degreesOfFreedom();
+        }
+    }
+
+    scalar meshVolume = sum(mesh_.cellVolumes());
+
+    if (Pstream::parRun())
+    {
+        reduce(totalLinearMomentum, sumOp<vector>());
+        reduce(totalAngularMomentum, sumOp<vector>());
+        reduce(maxVelocityMag, maxOp<scalar>());
+        reduce(totalMass, sumOp<scalar>());
+        reduce(totalLinearKE, sumOp<scalar>());
+        reduce(totalAngularKE, sumOp<scalar>());
+        reduce(totalPE, sumOp<scalar>());
+        reduce(totalrDotf, sumOp<scalar>());
+        reduce(nMols, sumOp<label>());
+        reduce(dofs, sumOp<label>());
+        reduce(meshVolume, sumOp<scalar>());
+    }
+
+    if (nMols)
+    {
+        Info<< "Number of molecules in " << this->name() << "  = "
+            << nMols << nl
+            << "Overall number density = "
+            << nMols/meshVolume << nl
+            << "Overall mass density = "
+            << totalMass/meshVolume << nl
+            << "Average linear momentum per molecule = "
+            << totalLinearMomentum/nMols << ' '
+            << mag(totalLinearMomentum)/nMols << nl
+            << "Average angular momentum per molecule = "
+            << totalAngularMomentum << ' '
+            << mag(totalAngularMomentum)/nMols << nl
+            << "maximum |velocity| = "
+            << maxVelocityMag << nl
+            << "Average linear KE per molecule = "
+            << totalLinearKE/nMols << nl
+            << "Average angular KE per molecule = "
+            << totalAngularKE/nMols << nl
+            << "Average PE per molecule = "
+            << totalPE/nMols << nl
+            << "Average TE per molecule = "
+            <<
+            (
+                  totalLinearKE
+                + totalAngularKE
+                + totalPE
+            )
+            /nMols
+            << endl;
+    }
+    else
+    {
+        Info<< "No molecules in " << this->name() << endl;
+    }
+}
+
+
 void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
 {
     OFstream os(fName);

src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.H

@@ -178,6 +178,9 @@ public:
             const scalar measuredTemperature
         );
 
+        //- Print cloud information
+        void info() const;
+
 
         // Access
 

src/lagrangian/molecularDynamics/potential/potential/potential.C

@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
+void Foam::potential::setSiteIdList(const dictionary& mdPropertiesDict)
 {
     DynamicList<word> siteIdList;
     DynamicList<word> pairPotentialSiteIdList;
@@ -36,14 +36,14 @@ void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
     {
         const word& id(idList_[i]);
 
-        if (!moleculePropertiesDict.found(id))
+        if (!mdPropertiesDict.found(id))
         {
             FatalErrorIn("potential::setSiteIdList(const dictionary&)")
                 << id << " molecule subDict not found"
                 << nl << abort(FatalError);
         }
 
-        const dictionary& molDict(moleculePropertiesDict.subDict(id));
+        const dictionary& molDict(mdPropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -117,7 +117,7 @@ void Foam::potential::potential::readPotentialDict()
         (
             IOobject
             (
-                "moleculeProperties",
+                "mdProperties",
                 mesh_.time().constant(),
                 mesh_,
                 IOobject::MUST_READ_IF_MODIFIED,
@@ -247,11 +247,11 @@ void Foam::potential::potential::readMdInitialiseDict
     IOdictionary& idListDict
 )
 {
-    IOdictionary moleculePropertiesDict
+    IOdictionary mdPropertiesDict
     (
         IOobject
         (
-            "moleculeProperties",
+            "mdProperties",
             mesh_.time().constant(),
             mesh_,
             IOobject::MUST_READ_IF_MODIFIED,
@@ -280,7 +280,7 @@ void Foam::potential::potential::readMdInitialiseDict
         {
             const word& id = latticeIds[i];
 
-            if (!moleculePropertiesDict.found(id))
+            if (!mdPropertiesDict.found(id))
             {
                 FatalErrorIn
                 (
@@ -290,7 +290,7 @@ void Foam::potential::potential::readMdInitialiseDict
                         "IOdictionary&"
                     ")"
                 )   << "Molecule type " << id
-                    << " not found in moleculeProperties dictionary." << nl
+                    << " not found in mdProperties dictionary." << nl
                     << abort(FatalError);
             }
 
@@ -322,7 +322,7 @@ void Foam::potential::potential::readMdInitialiseDict
 
                 List<word> siteIds
                 (
-                    moleculePropertiesDict.subDict(id).lookup("siteIds")
+                    mdPropertiesDict.subDict(id).lookup("siteIds")
                 );
 
                 if (findIndex(siteIds, tetherSiteId) != -1)
@@ -359,7 +359,7 @@ void Foam::potential::potential::readMdInitialiseDict
 
     idListDict.add("tetherSiteIdList", tetherSiteIdList);
 
-    setSiteIdList(moleculePropertiesDict);
+    setSiteIdList(mdPropertiesDict);
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/lagrangian/molecularDynamics/potential/potential/potential.H

@@ -76,7 +76,7 @@ class potential
 
     // Private Member Functions
 
-        void setSiteIdList(const dictionary& moleculePropertiesDict);
+        void setSiteIdList(const dictionary& mdPropertiesDict);
 
         void readPotentialDict();