diff --git a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C index 6c474208d5..754f078b30 100644 --- a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C +++ b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C @@ -74,7 +74,7 @@ int main(int argc, char *argv[]) Info<< nl << "Total number of molecules added: " << totalMolecules << nl << endl; - IOstream::defaultPrecision(12); + IOstream::defaultPrecision(15); if (!mesh.write()) { diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C index 87bd90ca9d..2d8cd8f00b 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C @@ -210,6 +210,14 @@ void Foam::molecule::setSitePositions(const constantProperties& constProps) } +void Foam::molecule::setSiteSizes(label size) +{ + sitePositions_.setSize(size); + + siteForces_.setSize(size); +} + + void Foam::molecule::hitProcessorPatch ( const processorPolyPatch&, diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H index 8e4c822efc..619fd63ef5 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H @@ -259,6 +259,8 @@ public: void setSitePositions(const constantProperties& constProps); + void setSiteSizes(label size); + // Access inline const tensor& Q() const; diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C index 858437c4f1..c15a629417 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C @@ -100,6 +100,21 @@ void Foam::moleculeCloud::buildConstProps() } +void Foam::moleculeCloud::setSiteSizesAndPositions() +{ + iterator mol(this->begin()); + + for (mol = this->begin(); mol != this->end(); ++mol) + { + const molecule::constantProperties& cP = constProps(mol().id()); + + mol().setSiteSizes(cP.nSites()); + + mol().setSitePositions(cP); + } +} + + void Foam::moleculeCloud::buildCellOccupancy() { forAll(cellOccupancy_, cO) @@ -499,11 +514,6 @@ void Foam::moleculeCloud::initialiseMolecules readScalar(zoneDict.lookup("temperature")) ); - const word velocityDistribution - ( - zoneDict.lookup("velocityDistribution") - ); - const vector bulkVelocity(zoneDict.lookup("bulkVelocity")); List latticeIds @@ -571,7 +581,7 @@ void Foam::moleculeCloud::initialiseMolecules else { FatalErrorIn("Foam::moleculeCloud::initialiseMolecules") - << "massDensity or numberDensity not specified " << nl + << "massDensity or numberDensity not specified " << nl << abort(FatalError); } @@ -697,20 +707,28 @@ void Foam::moleculeCloud::initialiseMolecules bool partOfLayerInBounds = false; // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // start of place molecules + // start of placement of molecules // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - vector unitCellLatticePosition(0,0,0); - - // Special treatment is required for the first position, - // i.e. iteration zero. - if (n == 0) { + // Special treatment is required for the first position, + // i.e. iteration zero. + + labelVector unitCellLatticePosition(0,0,0); + forAll(latticePositions, p) { + label id = findIndex(pot_.idList(), latticeIds[p]); + const vector& latticePosition = - unitCellLatticePosition + latticePositions[p]; + vector + ( + unitCellLatticePosition.x(), + unitCellLatticePosition.y(), + unitCellLatticePosition.z() + ) + + latticePositions[p]; point globalPosition = anchor + (R & (latticeCellShape & latticePosition)); @@ -722,17 +740,163 @@ void Foam::moleculeCloud::initialiseMolecules if (findIndex(zone, cell) != -1) { - Info<< "add molecule at " - << globalPosition - << " to cell " - << cell - << endl; + createMolecule + ( + globalPosition, + cell, + id, + tethered, + temperature, + bulkVelocity + ); + } + } + } + else + { + // Place top and bottom caps. + + labelVector unitCellLatticePosition(0,0,0); + + for + ( + unitCellLatticePosition.z() = -n; + unitCellLatticePosition.z() <= n; + unitCellLatticePosition.z() += 2*n + ) + { + for + ( + unitCellLatticePosition.y() = -n; + unitCellLatticePosition.y() <= n; + unitCellLatticePosition.y()++ + ) + { + for + ( + unitCellLatticePosition.x() = -n; + unitCellLatticePosition.x() <= n; + unitCellLatticePosition.x()++ + ) + { + forAll(latticePositions, p) + { + label id = findIndex + ( + pot_.idList(), + latticeIds[p] + ); + + const vector& latticePosition = + vector + ( + unitCellLatticePosition.x(), + unitCellLatticePosition.y(), + unitCellLatticePosition.z() + ) + + latticePositions[p]; + + point globalPosition = anchor + + (R + &(latticeCellShape & latticePosition)); + + partOfLayerInBounds = + mesh_.bounds().contains(globalPosition); + + label cell = mesh_.findCell + ( + globalPosition + ); + + if (findIndex(zone, cell) != -1) + { + createMolecule + ( + globalPosition, + cell, + id, + tethered, + temperature, + bulkVelocity + ); + } + } + } + } + } + + for + ( + unitCellLatticePosition.z() = -(n-1); + unitCellLatticePosition.z() <= (n-1); + unitCellLatticePosition.z()++ + ) + { + for (label iR = 0; iR <= 2*n -1; iR++) + { + unitCellLatticePosition.x() = n; + + unitCellLatticePosition.y() = -n + (iR + 1); + + for (label iK = 0; iK < 4; iK++) + { + forAll(latticePositions, p) + { + label id = findIndex + ( + pot_.idList(), + latticeIds[p] + ); + + const vector& latticePosition = + vector + ( + unitCellLatticePosition.x(), + unitCellLatticePosition.y(), + unitCellLatticePosition.z() + ) + + latticePositions[p]; + + point globalPosition = anchor + + (R + &(latticeCellShape & latticePosition)); + + partOfLayerInBounds = + mesh_.bounds().contains(globalPosition); + + label cell = mesh_.findCell + ( + globalPosition + ); + + if (findIndex(zone, cell) != -1) + { + createMolecule + ( + globalPosition, + cell, + id, + tethered, + temperature, + bulkVelocity + ); + } + } + + unitCellLatticePosition = + labelVector + ( + - unitCellLatticePosition.y(), + unitCellLatticePosition.x(), + unitCellLatticePosition.z() + ); + } } } } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // end of place molecules + // end of placement of molecules // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if @@ -741,7 +905,7 @@ void Foam::moleculeCloud::initialiseMolecules && !partOfLayerInBounds ) { - WarningIn("molConfig::createMolecules()") + WarningIn("Foam::moleculeCloud::initialiseMolecules()") << "A whole layer of unit cells was placed " << "outside the bounds of the mesh, but no " << "molecules have been placed in zone '" @@ -770,6 +934,75 @@ void Foam::moleculeCloud::initialiseMolecules } } + +void Foam::moleculeCloud::createMolecule +( + const point& position, + label cell, + label id, + bool tethered, + scalar temperature, + const vector& bulkVelocity +) +{ + const Cloud& cloud = *this; + + if (cell == -1) + { + cell = mesh_.findCell(position); + } + + if (cell == -1) + { + FatalErrorIn("Foam::moleculeCloud::createMolecule") + << "Position specified does not correspond to a mesh cell." << nl + << abort(FatalError); + } + + point specialPosition(vector::zero); + + label special = 0; + + if (tethered) + { + specialPosition = position; + + special = molecule::SPECIAL_TETHERED; + } + + const molecule::constantProperties& cP(constProps(id)); + + vector v = equipartitionLinearVelocity(temperature, cP.mass()); + + v += bulkVelocity; + + vector pi = equipartitionAngularMomentum + ( + temperature, + cP.momentOfInertia() + ); + + addParticle + ( + new molecule + ( + cloud, + position, + cell, + I, + v, + vector::zero, + pi, + vector::zero, + specialPosition, + constProps(id), + special, + id + ) + ); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::moleculeCloud::moleculeCloud @@ -783,13 +1016,14 @@ Foam::moleculeCloud::moleculeCloud pot_(pot), cellOccupancy_(mesh_.nCells()), il_(mesh_, pot_.pairPotentials().rCutMaxSqr(), false), - constPropList_() + constPropList_(), + rndGen_(clock::getTime()) { molecule::readFields(*this); buildConstProps(); - buildCellOccupancy(); + setSiteSizesAndPositions(); removeHighEnergyOverlaps(); @@ -808,7 +1042,8 @@ Foam::moleculeCloud::moleculeCloud mesh_(mesh), pot_(pot), il_(mesh_), - constPropList_() + constPropList_(), + rndGen_(clock::getTime()) { molecule::readFields(*this); @@ -818,29 +1053,6 @@ Foam::moleculeCloud::moleculeCloud initialiseMolecules(mdInitialiseDict); } -// int i; - -// const Cloud& cloud = *this; - -// for (i=0; i constPropList_; + Random rndGen_; + // Private Member Functions void buildConstProps(); + void setSiteSizesAndPositions(); + //- Determine which molecules are in which cells void buildCellOccupancy(); @@ -115,6 +121,28 @@ private: const IOdictionary& mdInitialiseDict ); + void createMolecule + ( + const point& position, + label cell, + label id, + bool tethered, + scalar temperature, + const vector& bulkVelocity + ); + + inline vector equipartitionLinearVelocity + ( + scalar temperature, + scalar mass + ); + + inline vector equipartitionAngularMomentum + ( + scalar temperature, + const diagTensor& momentOfInertia + ); + //- Disallow default bitwise copy construct moleculeCloud(const moleculeCloud&); @@ -178,6 +206,8 @@ public: inline const molecule::constantProperties& constProps(label id) const; + inline Random& rndGen(); + // Member Operators //- Write fields diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index e30bd89ac0..7cde56b4ce 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -367,7 +367,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit if ( mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag)) - > pot_.potentialEnergyLimit() + > pot_.potentialEnergyLimit() ) { return true; @@ -510,7 +510,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit chargeReal*chargeRef *electrostatic.energy(rsRealsRefMag) ) - > pot_.potentialEnergyLimit() + > pot_.potentialEnergyLimit() ) { return true; @@ -524,6 +524,35 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit } +inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity +( + scalar temperature, + scalar mass +) +{ + return sqrt(kb*temperature/mass)*vector + ( + rndGen_.GaussNormal(), + rndGen_.GaussNormal(), + rndGen_.GaussNormal() + ); +} + + +inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum +( + scalar temperature, + const diagTensor& momentOfInertia +) +{ + return vector + ( + sqrt(kb*temperature*momentOfInertia.xx())*rndGen_.GaussNormal(), + sqrt(kb*temperature*momentOfInertia.yy())*rndGen_.GaussNormal(), + sqrt(kb*temperature*momentOfInertia.zz())*rndGen_.GaussNormal() + ); +} + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const @@ -566,4 +595,10 @@ inline const Foam::molecule::constantProperties& } +inline Foam::Random& Foam::moleculeCloud::rndGen() +{ + return rndGen_; +} + + // ************************************************************************* //