Commit af88d9be authored by Andrew Heather's avatar Andrew Heather
Browse files

updating molecularDynamics functionality with Graham's latest changes

parent 7b62e070
EXE_INC = \
-I$(LIB_SRC)/lagrangian/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
......@@ -8,4 +9,5 @@ EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \
-llagrangian \
-lmolecule
-lmolecule \
-lpotential
EXE_INC = \
-I$(LIB_SRC)/lagrangian/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
......@@ -8,4 +9,5 @@ EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \
-llagrangian \
-lmolecule
-lmolecule \
-lpotential
latticeStructures = latticeStructures
velocityDistributions = velocityDistributions
createMolecules.C
molConfig.C
genMolConfig.C
EXE = $(FOAM_APPBIN)/molConfig
EXE_INC = \
-I$(latticeStructures) \
-I$(velocityDistributions) \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lmeshTools \
-ldynamicMesh \
-lfiniteVolume \
-llagrangian \
-lmolecule \
-lpotential
for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
{
// Remove bulk momentum introduced by random numbers and add
// desired bulk velocity
// For systems with molecules of significantly differing masses, this may
// need to be an iterative process or employ a better algorithm for
// removing an appropriate share of the excess momentum from each molecule.
initialVelocities(molN) += bulkVelocity - momentumSum/totalZoneMols/mass;
}
// momentumSum = vector::zero;
//
// for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
// {
// momentumSum += mass*initialVelocities(molN);
// }
//
// Info << "Check momentum adjustment: " << momentumSum << endl;
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*----------------------------------------------------------------------------*/
#include "molConfig.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::molConfig::createMolecules()
{
Info<< nl << "Creating molecules from zone specifications\n" << endl;
DynamicList<vector> initialPositions(0);
DynamicList<label> initialIds(0);
DynamicList<scalar> initialMasses(0);
DynamicList<label> initialCelli(0);
DynamicList<vector> initialVelocities(0);
DynamicList<vector> initialAccelerations(0);
DynamicList<label> initialTethered(0);
DynamicList<vector> initialTetherPositions(0);
label totalMols = 0;
label idAssign;
Random rand(clock::getTime());
// * * * * * * * * Building the IdList * * * * * * * * * //
//Notes: - each processor will have an identical idList_.
// - The order of id's inside the idList_ depends on the order
// of subDicts inside the molConigDict.
Info<< "Building the idList: " ;
forAll(molConfigDescription_.toc(), cZs)
{
word subDictName (molConfigDescription_.toc()[cZs]);
word iD (molConfigDescription_.subDict(subDictName).lookup("id"));
if (findIndex(idList_,iD) == -1)
{
idList_.append(iD);
}
}
forAll(idList_, i)
{
Info << " " << idList_[i];
}
Info << nl << endl;
// * * * * * * * * Filling the Mesh * * * * * * * * * //
const cellZoneMesh& cellZoneI = mesh_.cellZones();
if (cellZoneI.size())
{
Info<< "Filling the zones with molecules." << nl << endl;
}
else
{
FatalErrorIn("void createMolecules()\n")
<< "No cellZones found in mesh description."
<< abort(FatalError);
}
forAll (cellZoneI, cZ)
{
if (cellZoneI[cZ].size())
{
if (!molConfigDescription_.found(cellZoneI[cZ].name()))
{
Info << "Zone specification subDictionary: "
<< cellZoneI[cZ].name() << " not found." << endl;
}
else
{
label n = 0;
label totalZoneMols = 0;
label molsPlacedThisIteration;
# include "readZoneSubDict.H"
idAssign = findIndex(idList_,id);
# include "startingPoint.H"
// Continue trying to place molecules as long as at
// least one molecule is placed in each iteration.
// The "|| totalZoneMols == 0" condition means that the
// algorithm will continue if the origin is outside the
// zone - it will cause an infinite loop if no molecules
// are ever placed by the algorithm.
if (latticeStructure != "empty")
{
while
(
molsPlacedThisIteration != 0
|| totalZoneMols == 0
)
{
molsPlacedThisIteration = 0;
bool partOfLayerInBounds = false;
# include "createPositions.H"
if
(
totalZoneMols == 0
&& !partOfLayerInBounds
)
{
WarningIn("molConfig::createMolecules()")
<< "A whole layer of unit cells was placed "
<< "outside the bounds of the mesh, but no "
<< "molecules have been placed in zone '"
<< cellZoneI[cZ].name()
<< "'. This is likely to be because the zone "
<< "has few cells ("
<< cellZoneI[cZ].size()
<< " in this case) and no lattice position "
<< "fell inside them. "
<< "Aborting filling this zone."
<< endl;
break;
}
totalZoneMols += molsPlacedThisIteration;
n++;
}
label molN;
for
(
molN = totalMols;
molN < totalMols + totalZoneMols;
molN++
)
{
initialIds.append(idAssign);
initialMasses.append(mass);
initialAccelerations.append(vector::zero);
if (tethered)
{
initialTethered.append(1);
initialTetherPositions.append
(
initialPositions[molN]
);
}
else
{
initialTethered.append(0);
initialTetherPositions.append(vector::zero);
}
}
# include "createVelocities.H"
# include "correctVelocities.H"
}
totalMols += totalZoneMols;
}
}
}
idList_.shrink();
positions_ = initialPositions;
positions_.setSize(initialPositions.size());
id_ = initialIds;
id_.setSize(initialIds.size());
mass_ = initialMasses;
mass_.setSize(initialMasses.size());
cells_ = initialCelli;
cells_.setSize(initialCelli.size());
U_ = initialVelocities;
U_.setSize(initialVelocities.size());
A_ = initialAccelerations;
A_.setSize(initialAccelerations.size());
tethered_ = initialTethered;
tethered_.setSize(initialTethered.size());
tetherPositions_ = initialTetherPositions;
tetherPositions_.setSize(initialTetherPositions.size());
nMol_ = totalMols;
}
// ************************************************************************* //
vector latticePosition;
vector globalPosition;
if (latticeStructure == "SC")
{
# include "SC.H"
}
else if (latticeStructure == "FCC")
{
# include "FCC.H"
}
else if (latticeStructure == "BCC")
{
# include "BCC.H"
}
else
{
FatalErrorIn("createPositions.H\n")
<< "latticeStructure " << latticeStructure
<< " not supported."
<< abort(FatalError);
}
vector velocity;
vector momentumSum = vector::zero;
if (velocityDistribution == "uniform")
{
# include "uniform.H"
}
if (velocityDistribution == "maxwellian")
{
# include "maxwellian.H"
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "molConfig.H"
#include "fvCFD.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
Info<< nl << "Reading molecular configuration description dictionary"
<< endl;
IOobject molConfigDescriptionIOobject
(
"molConfigDict",
runTime.system(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
if (!molConfigDescriptionIOobject.headerOk())
{
FatalErrorIn(args.executable())
<< "Cannot find molConfig description file " << nl
<< args.caseName()/runTime.system()/"molConfig"/"molConfigDict"
<< nl << exit(FatalError);
}
IOdictionary molConfigDescription(molConfigDescriptionIOobject);
// Create molCloud, registering object with mesh
Info<< nl << "Creating molecular configuration" << endl;
molConfig molecules(molConfigDescription, mesh);
label totalMolecules = molecules.nMol();
if (Pstream::parRun())
{
reduce(totalMolecules, sumOp<label>());
}
Info<< nl << "Total number of molecules added: " << totalMolecules
<< nl << endl;
moleculeCloud molCloud
(
mesh,
molecules.nMol(),
molecules.id(),
molecules.mass(),
molecules.positions(),
molecules.cells(),
molecules.U(),
molecules.A(),
molecules.tethered(),
molecules.tetherPositions()
);
IOdictionary idListDict
(
IOobject
(
"idList",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
)
);
idListDict.add("idList", molecules.molIdList());
IOstream::defaultPrecision(12);
Info << nl << "Writing molecular configuration" << endl;
if (!mesh.write())
{
FatalErrorIn(args.executable())
<< "Failed writing moleculeCloud."
<< nl << exit(FatalError);
}
Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
Info << nl << "End\n" << endl;
return 0;
}
// ************************************************************************* //
labelVector iN(0,0,0);
vector gap = (vector::one)*pow((numberDensity/2.0),-(1.0/3.0));
#include "origin.H"
// Info<< "gap = " << gap << endl;
// Special treatment is required for the first position, i.e. iteration zero.
if (n == 0)
{
latticePosition.x() = (iN.x()*gap.x());
latticePosition.y() = (iN.y()*gap.y());
latticePosition.z() = (iN.z()*gap.z());
// Placing 2 molecules in each unit cell, using the algorithm from
// D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
for (label iU = 0; iU < 2; iU++)
{
vector unitCellLatticePosition = latticePosition;
if (iU == 1)
{
unitCellLatticePosition += 0.5 * gap;
}
if (originSpecifies == "corner")
{
unitCellLatticePosition -= 0.25*gap;
}
// Info << nl << n << ", " << unitCellLatticePosition;
globalPosition =
origin + transform(latticeToGlobal,unitCellLatticePosition);
partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
if
(
findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition))
!= -1
)
{
molsPlacedThisIteration++;
initialPositions.append(globalPosition);
initialCelli.append(mesh_.findCell(globalPosition));
}
}
}
else
{
// Place top and bottom caps.
for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
{
for (iN.y() = -n; iN.y() <= n; iN.y()++)
{
for (iN.x() = -n; iN.x() <= n; iN.x()++)
{
latticePosition.x() = (iN.x() * gap.x());
latticePosition.y() = (iN.y() * gap.y());
latticePosition.z() = (iN.z() * gap.z());
for (label iU = 0; iU < 2; iU++)
{
vector unitCellLatticePosition = latticePosition;
if (iU == 1)
{
unitCellLatticePosition += 0.5*gap;
}
if(originSpecifies == "corner")
{
unitCellLatticePosition -= 0.25*gap;
}
// Info << nl << iN << ", " << unitCellLatticePosition;
globalPosition =
origin
+ transform(latticeToGlobal,