diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C index 660c12a5b8..65157dd2b4 100644 --- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C +++ b/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; diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/files similarity index 100% rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/files diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/options similarity index 100% rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/options diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/mdEquilibrationFoam.C similarity index 100% rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/mdEquilibrationFoam.C diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/readmdEquilibrationDict.H similarity index 100% rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/readmdEquilibrationDict.H diff --git a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C index 3551e46133..5d06f35fb6 100644 --- a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C +++ b/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); } diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C index 431f37aeb4..5647869233 100644 --- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C +++ b/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; +} + // ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H index 0ffb63b0a8..c53eba2ee8 100644 --- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H +++ b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H @@ -35,14 +35,23 @@ SourceFiles #define constPropSite_H #include "vector.H" -#include "IFstream.H" -#include "OFstream.H" +#include "IOstreams.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&); }; diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H index dbeb395920..b4f1ad9177 100644 --- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H +++ b/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); +} + + // ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H b/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H index faf24d1f1c..efb9cf52f4 100644 --- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H +++ b/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H @@ -57,11 +57,12 @@ inline Foam::polyatomic::constantProperties::constantProperties List pairPotentialIds(dict.lookup("pairPotentialSiteIds")); + sites_.setSize(siteReferencePositions.size()); + if ( - ( siteIdNames.size() != sites_.size() ) - || ( siteReferencePositions.size() != sites_.size() ) - || ( siteCharges.size() != sites_.size() ) + (siteIdNames.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() - << nl << abort(FatalError); + << "referencePositions are not the same: " << nl + << 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 ( diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C b/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C index 0d2a7c4f9a..0ad96b90c8 100644 --- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C +++ b/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& idList(pot_.idList()); @@ -48,11 +48,11 @@ void Foam::polyatomicCloud::buildConstProps() const List& 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 siteIdNames = molDict.lookup("siteIds"); @@ -1183,6 +1183,8 @@ void Foam::polyatomicCloud::evolve() polyatomic::trackingData td3(*this, 3); Cloud::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()); + reduce(totalAngularMomentum, sumOp()); + reduce(maxVelocityMag, maxOp()); + reduce(totalMass, sumOp()); + reduce(totalLinearKE, sumOp()); + reduce(totalAngularKE, sumOp()); + reduce(totalPE, sumOp()); + reduce(totalrDotf, sumOp()); + reduce(nMols, sumOp