zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Finished lattice generation and molecule creation. Functions added to set the sizes of the sitePositions_ and siteForces_ members after construction from disk.

Browse Source
This commit is contained in:
graham
2009-01-07 15:33:16 +00:00
parent 23e1c04cca
commit d61af129f1
6 changed files with 336 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/utilities/preProcessing/mdInitialise/mdInitialise.C
View File

@ -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())
{

8
src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
View File

@ -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&,

2
src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
View File

@ -259,6 +259,8 @@ public:
void setSitePositions(const constantProperties& constProps);
void setSiteSizes(label size);
// Access
inline const tensor& Q() const;

304
src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
View File

@ -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<word> 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<molecule>& 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<molecule>& cloud = *this;
// for (i=0; i<nMol; i++)
// {
// addParticle
// (
// new molecule
// (
// cloud,
// positions[i],
// cells[i],
// mass[i],
// U[i],
// A[i],
// tetherPositions[i],
// tethered[i],
// id[i]
// )
// );
// }
// };
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

30
src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
View File

@ -42,6 +42,8 @@ SourceFiles
#include "IOdictionary.H"
#include "potential.H"
#include "interactionLists.H"
#include "labelVector.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -71,10 +73,14 @@ private:
List<molecule::constantProperties> 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

39
src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
View File

@ -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_;
}
// ************************************************************************* //