mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
removed molecular dynamics functionality (being updated)
This commit is contained in:
andy
@ -1,8 +0,0 @@
|
||||
latticeStructures = latticeStructures
|
||||
velocityDistributions = velocityDistributions
|
||||
|
||||
createMolecules.C
|
||||
molConfig.C
|
||||
genMolConfig.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/molConfig
|
||||
@ -1,21 +0,0 @@
|
||||
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;
|
||||
@ -1,253 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,26 +0,0 @@
|
||||
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);
|
||||
}
|
||||
@ -1,13 +0,0 @@
|
||||
vector velocity;
|
||||
|
||||
vector momentumSum = vector::zero;
|
||||
|
||||
if (velocityDistribution == "uniform")
|
||||
{
|
||||
# include "uniform.H"
|
||||
}
|
||||
|
||||
if (velocityDistribution == "maxwellian")
|
||||
{
|
||||
# include "maxwellian.H"
|
||||
}
|
||||
@ -1,129 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,179 +0,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,unitCellLatticePosition);
|
||||
|
||||
partOfLayerInBounds =
|
||||
mesh_.bounds().contains(globalPosition);
|
||||
|
||||
if
|
||||
(
|
||||
findIndex
|
||||
(
|
||||
mesh_.cellZones()[cZ],
|
||||
mesh_.findCell(globalPosition)
|
||||
)
|
||||
!= -1)
|
||||
{
|
||||
molsPlacedThisIteration++;
|
||||
|
||||
initialPositions.append(globalPosition);
|
||||
|
||||
initialCelli.append(mesh_.findCell(globalPosition));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Placing sides
|
||||
|
||||
for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
|
||||
{
|
||||
for (label iR = 0; iR <= 2*n -1; iR++)
|
||||
{
|
||||
latticePosition.x() = (n*gap.x());
|
||||
|
||||
latticePosition.y() = ((-n + (iR + 1))*gap.y());
|
||||
|
||||
latticePosition.z() = (iN.z() * gap.z());
|
||||
|
||||
for (label iK = 0; iK < 4; iK++)
|
||||
{
|
||||
for (label iU = 0; iU < 2; iU++)
|
||||
{
|
||||
vector unitCellLatticePosition = latticePosition;
|
||||
|
||||
if (iU == 1)
|
||||
{
|
||||
unitCellLatticePosition += 0.5 * gap;
|
||||
}
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
unitCellLatticePosition -= 0.25*gap;
|
||||
}
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
|
||||
latticePosition =
|
||||
vector
|
||||
(
|
||||
- latticePosition.y(),
|
||||
latticePosition.x(),
|
||||
latticePosition.z()
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1,217 +0,0 @@
|
||||
labelVector iN(0,0,0);
|
||||
|
||||
vector gap = (vector::one)*pow((numberDensity/4.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 4 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 < 4; iU++)
|
||||
{
|
||||
vector unitCellLatticePosition = latticePosition;
|
||||
|
||||
if (iU != 3)
|
||||
{
|
||||
if (iU != 0)
|
||||
{
|
||||
unitCellLatticePosition.x() += 0.5 * gap.x();
|
||||
}
|
||||
|
||||
if (iU != 1)
|
||||
{
|
||||
unitCellLatticePosition.y() += 0.5 * gap.y();
|
||||
}
|
||||
|
||||
if (iU != 2)
|
||||
{
|
||||
unitCellLatticePosition.z() += 0.5 * gap.z();
|
||||
}
|
||||
}
|
||||
|
||||
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 < 4; iU++)
|
||||
{
|
||||
vector unitCellLatticePosition = latticePosition;
|
||||
|
||||
if (iU != 3)
|
||||
{
|
||||
if (iU != 0)
|
||||
{
|
||||
unitCellLatticePosition.x() += 0.5 * gap.x();
|
||||
}
|
||||
|
||||
if (iU != 1)
|
||||
{
|
||||
unitCellLatticePosition.y() += 0.5 * gap.y();
|
||||
}
|
||||
|
||||
if (iU != 2)
|
||||
{
|
||||
unitCellLatticePosition.z() += 0.5 * gap.z();
|
||||
}
|
||||
}
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
unitCellLatticePosition -= 0.25*gap;
|
||||
}
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Placing sides
|
||||
|
||||
for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
|
||||
{
|
||||
for (label iR = 0; iR <= 2*n -1; iR++)
|
||||
{
|
||||
latticePosition.x() = (n * gap.x());
|
||||
|
||||
latticePosition.y() = ((-n + (iR + 1)) * gap.y());
|
||||
|
||||
latticePosition.z() = (iN.z() * gap.z());
|
||||
|
||||
for (label iK = 0; iK < 4; iK++)
|
||||
{
|
||||
for (label iU = 0; iU < 4; iU++)
|
||||
{
|
||||
vector unitCellLatticePosition = latticePosition;
|
||||
|
||||
if (iU != 3)
|
||||
{
|
||||
if (iU != 0)
|
||||
{
|
||||
unitCellLatticePosition.x() += 0.5 * gap.x();
|
||||
}
|
||||
|
||||
if (iU != 1)
|
||||
{
|
||||
unitCellLatticePosition.y() += 0.5 * gap.y();
|
||||
}
|
||||
|
||||
if (iU != 2)
|
||||
{
|
||||
unitCellLatticePosition.z() += 0.5 * gap.z();
|
||||
}
|
||||
}
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
unitCellLatticePosition -= 0.25*gap;
|
||||
}
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
|
||||
latticePosition =
|
||||
vector
|
||||
(
|
||||
- latticePosition.y(),
|
||||
latticePosition.x(),
|
||||
latticePosition.z()
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1,127 +0,0 @@
|
||||
labelVector iN(0,0,0);
|
||||
|
||||
vector gap = (vector::one)*pow(numberDensity, -(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 = vector::zero;
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
latticePosition += 0.5*gap;
|
||||
}
|
||||
|
||||
globalPosition = origin + transform(latticeToGlobal,latticePosition);
|
||||
|
||||
partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
|
||||
|
||||
if (findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition)) != -1)
|
||||
{
|
||||
molsPlacedThisIteration++;
|
||||
|
||||
initialPositions.append(globalPosition);
|
||||
|
||||
initialCelli.append(mesh_.findCell(globalPosition));
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
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());
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
latticePosition += 0.5*gap;
|
||||
}
|
||||
|
||||
globalPosition =
|
||||
origin + transform(latticeToGlobal,latticePosition);
|
||||
|
||||
partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
|
||||
|
||||
if
|
||||
(
|
||||
findIndex
|
||||
(
|
||||
mesh_.cellZones()[cZ],
|
||||
mesh_.findCell(globalPosition)
|
||||
)
|
||||
!= -1
|
||||
)
|
||||
{
|
||||
molsPlacedThisIteration++;
|
||||
|
||||
initialPositions.append(globalPosition);
|
||||
|
||||
initialCelli.append(mesh_.findCell(globalPosition));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
tensor quarterRotate(EulerCoordinateRotation(-90, 0, 0, true).R());
|
||||
|
||||
iN.x() = n;
|
||||
for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
|
||||
{
|
||||
for (iN.y() = -(n-1); iN.y() <= n; iN.y()++)
|
||||
{
|
||||
latticePosition.x() = (iN.x()*gap.x());
|
||||
|
||||
latticePosition.y() = (iN.y()*gap.y());
|
||||
|
||||
latticePosition.z() = (iN.z()*gap.z());
|
||||
|
||||
for (label iR = 0; iR < 4; iR++)
|
||||
{
|
||||
vector offsetCorrectedLatticePosition = latticePosition;
|
||||
|
||||
if (originSpecifies == "corner")
|
||||
{
|
||||
offsetCorrectedLatticePosition += 0.5*gap;
|
||||
}
|
||||
|
||||
globalPosition =
|
||||
origin
|
||||
+ transform(latticeToGlobal,offsetCorrectedLatticePosition);
|
||||
|
||||
partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
|
||||
|
||||
if
|
||||
(
|
||||
findIndex
|
||||
(
|
||||
mesh_.cellZones()[cZ],
|
||||
mesh_.findCell(globalPosition)
|
||||
)
|
||||
!= -1
|
||||
)
|
||||
{
|
||||
molsPlacedThisIteration++;
|
||||
|
||||
initialPositions.append(globalPosition);
|
||||
|
||||
initialCelli.append(mesh_.findCell(globalPosition));
|
||||
}
|
||||
|
||||
latticePosition = transform(quarterRotate,latticePosition);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,50 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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"
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::molConfig::molConfig
|
||||
(
|
||||
IOdictionary& molConfigDescription,
|
||||
const polyMesh& mesh
|
||||
)
|
||||
:
|
||||
molConfigDescription_(molConfigDescription),
|
||||
mesh_(mesh)
|
||||
{
|
||||
createMolecules();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::molConfig::~molConfig()
|
||||
{}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,147 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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
|
||||
|
||||
Class
|
||||
Foam::molConfig
|
||||
|
||||
Description
|
||||
|
||||
SourceFiles
|
||||
molConfigI.H
|
||||
molConfig.C
|
||||
molConfigIO.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef molConfig_H
|
||||
#define molConfig_H
|
||||
|
||||
#include "labelVector.H"
|
||||
#include "scalar.H"
|
||||
#include "vector.H"
|
||||
#include "labelField.H"
|
||||
#include "scalarField.H"
|
||||
#include "vectorField.H"
|
||||
#include "IOField.H"
|
||||
#include "EulerCoordinateRotation.H"
|
||||
#include "Random.H"
|
||||
|
||||
#include "Time.H"
|
||||
#include "IOdictionary.H"
|
||||
#include "IOstreams.H"
|
||||
#include "moleculeCloud.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class molConfig Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class molConfig
|
||||
{
|
||||
// Private data
|
||||
|
||||
const IOdictionary& molConfigDescription_;
|
||||
|
||||
const polyMesh& mesh_;
|
||||
|
||||
DynamicList<word> idList_;
|
||||
|
||||
labelField id_;
|
||||
|
||||
scalarField mass_;
|
||||
|
||||
vectorField positions_;
|
||||
|
||||
labelField cells_;
|
||||
|
||||
vectorField U_;
|
||||
|
||||
vectorField A_;
|
||||
|
||||
labelField tethered_;
|
||||
|
||||
vectorField tetherPositions_;
|
||||
|
||||
label nMol_;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from IOdictionary and mesh
|
||||
molConfig(IOdictionary&, const polyMesh&);
|
||||
|
||||
|
||||
// Destructor
|
||||
|
||||
~molConfig();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
void createMolecules();
|
||||
|
||||
|
||||
// Access
|
||||
|
||||
inline const List<word>& molIdList() const;
|
||||
|
||||
inline const labelField& id() const;
|
||||
|
||||
inline const scalarField& mass() const;
|
||||
|
||||
inline const vectorField& positions() const;
|
||||
|
||||
inline const labelField& cells() const;
|
||||
|
||||
inline const vectorField& U() const;
|
||||
|
||||
inline const vectorField& A() const;
|
||||
|
||||
inline const labelField& tethered() const;
|
||||
|
||||
inline const vectorField& tetherPositions() const;
|
||||
|
||||
inline label nMol() const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "molConfigI.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,98 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
inline const List<word>& molConfig::molIdList() const
|
||||
{
|
||||
return idList_;
|
||||
}
|
||||
|
||||
|
||||
inline const labelField& molConfig::id() const
|
||||
{
|
||||
return id_;
|
||||
}
|
||||
|
||||
|
||||
inline const scalarField& molConfig::mass() const
|
||||
{
|
||||
return mass_;
|
||||
}
|
||||
|
||||
|
||||
inline const vectorField& molConfig::positions() const
|
||||
{
|
||||
return positions_;
|
||||
}
|
||||
|
||||
|
||||
inline const labelField& molConfig::cells() const
|
||||
{
|
||||
return cells_;
|
||||
}
|
||||
|
||||
|
||||
inline const vectorField& molConfig::U() const
|
||||
{
|
||||
return U_;
|
||||
}
|
||||
|
||||
|
||||
inline const vectorField& molConfig::A() const
|
||||
{
|
||||
return A_;
|
||||
}
|
||||
|
||||
|
||||
inline const labelField& molConfig::tethered() const
|
||||
{
|
||||
return tethered_;
|
||||
}
|
||||
|
||||
|
||||
inline const vectorField& molConfig::tetherPositions() const
|
||||
{
|
||||
return tetherPositions_;
|
||||
}
|
||||
|
||||
|
||||
inline label molConfig::nMol() const
|
||||
{
|
||||
return nMol_;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,49 +0,0 @@
|
||||
// Please refer to notes
|
||||
|
||||
// 1. Determine the unit cell dimensions: xU, yU and zU
|
||||
|
||||
const scalar xU = gap.x();
|
||||
const scalar yU = gap.y();
|
||||
const scalar zU = gap.z();
|
||||
|
||||
// 2. Determine the anchorPoint co-ordinates: xA, yA and zA
|
||||
|
||||
const scalar xA = anchorPoint.x();
|
||||
const scalar yA = anchorPoint.y();
|
||||
const scalar zA = anchorPoint.z();
|
||||
|
||||
// 3. Determine the vector rAB from global co-ordinate system:
|
||||
|
||||
const vector rAB((xMid - xA), (yMid - yA), (zMid - zA));
|
||||
|
||||
// 4. Transform vector rAS into lattice co-ordinate system:
|
||||
|
||||
const vector rASTransf = transform(latticeToGlobal.T(), rAB);
|
||||
|
||||
// Info << "The vector rAS = " << rAS << endl;
|
||||
// Info << "The vector rAStransf = " << rAStransf << endl;
|
||||
|
||||
// 5. Calculate the integer values: ni, nj and nk
|
||||
scalar nIscalar = rASTransf.x()/xU;
|
||||
scalar nJscalar = rASTransf.y()/yU;
|
||||
scalar nKscalar = rASTransf.z()/zU;
|
||||
|
||||
// Info << "The nI, nJ, nK values before are: " << nIscalar <<" "<< nJscalar <<" "<< nKscalar << endl;
|
||||
|
||||
label nI = label(nIscalar + 0.5*sign(nIscalar));
|
||||
label nJ = label(nJscalar + 0.5*sign(nJscalar));
|
||||
label nK = label(nKscalar + 0.5*sign(nKscalar));
|
||||
|
||||
// Info << "The nI, nJ, nK values after are: " << nI <<" "<< nJ <<" "<< nK << endl;
|
||||
|
||||
// 6. Calculate the corrected starting point, rAC (in the lattice co-ordinate system):
|
||||
const vector rAC((nI*xU), (nJ*yU), (nK*zU));
|
||||
|
||||
// 7. Transform the corrected starting point in the global co-ordinate system, rC:
|
||||
const vector rC = anchorPoint + transform(latticeToGlobal, rAC);
|
||||
|
||||
|
||||
const vector& origin = rC;
|
||||
|
||||
// Pout << "The Corrected Starting Point: " << origin << endl;
|
||||
|
||||
@ -1,93 +0,0 @@
|
||||
|
||||
// Info << "Zone description subDict " << cZ <<": " << cellZoneI[cZ].name() << endl;
|
||||
|
||||
const dictionary& subDictI =
|
||||
molConfigDescription_.subDict(cellZoneI[cZ].name());
|
||||
|
||||
const scalar temperature(readScalar(subDictI.lookup("temperature")));
|
||||
|
||||
const word velocityDistribution(subDictI.lookup("velocityDistribution"));
|
||||
|
||||
const vector bulkVelocity(subDictI.lookup("bulkVelocity"));
|
||||
|
||||
const word id(subDictI.lookup("id"));
|
||||
|
||||
const scalar mass(readScalar(subDictI.lookup("mass")));
|
||||
|
||||
scalar numberDensity_read(0.0);
|
||||
|
||||
if (subDictI.found("numberDensity"))
|
||||
{
|
||||
numberDensity_read = readScalar(subDictI.lookup("numberDensity"));
|
||||
}
|
||||
else if (subDictI.found("massDensity"))
|
||||
{
|
||||
numberDensity_read = readScalar(subDictI.lookup("massDensity"))/mass;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn("readZoneSubDict.H\n")
|
||||
<< "massDensity or numberDensity not specified " << nl
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
const scalar numberDensity(numberDensity_read);
|
||||
|
||||
const word latticeStructure(subDictI.lookup("latticeStructure"));
|
||||
|
||||
const vector anchorPoint(subDictI.lookup("anchor"));
|
||||
|
||||
const word originSpecifies(subDictI.lookup("anchorSpecifies"));
|
||||
|
||||
if
|
||||
(
|
||||
originSpecifies != "corner"
|
||||
&& originSpecifies != "molecule"
|
||||
)
|
||||
{
|
||||
FatalErrorIn("readZoneSubDict.H\n")
|
||||
<< "anchorSpecifies must be either 'corner' or 'molecule', found "
|
||||
<< originSpecifies << nl
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
bool tethered = false;
|
||||
|
||||
if (subDictI.found("tethered"))
|
||||
{
|
||||
tethered = Switch(subDictI.lookup("tethered"));
|
||||
}
|
||||
|
||||
const vector orientationAngles(subDictI.lookup("orientationAngles"));
|
||||
|
||||
scalar phi(orientationAngles.x()*mathematicalConstant::pi/180.0);
|
||||
scalar theta(orientationAngles.y()*mathematicalConstant::pi/180.0);
|
||||
scalar psi(orientationAngles.z()*mathematicalConstant::pi/180.0);
|
||||
|
||||
const tensor latticeToGlobal
|
||||
(
|
||||
cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
|
||||
cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
|
||||
sin(psi)*sin(theta),
|
||||
- sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
|
||||
- sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
|
||||
cos(psi)*sin(theta),
|
||||
sin(theta)*sin(phi),
|
||||
- sin(theta)*cos(phi),
|
||||
cos(theta)
|
||||
);
|
||||
|
||||
// Info << "\tcells: " << cellZoneI[cZ].size() << endl;
|
||||
// Info << "\tnumberDensity: " << numberDensity << endl;
|
||||
// Info << "\ttemperature: " << temperature << endl;
|
||||
// Info << "\tvelocityDistribution: " << velocityDistribution << endl;
|
||||
// Info << "\tbulkVelocity: " << bulkVelocity << endl;
|
||||
// Info << "\tid: " << id << endl;
|
||||
// Info << "\tmass: " << mass << endl;
|
||||
// Info << "\tlatticeStructure: " << latticeStructure << endl;
|
||||
// Info << "\tanchor: " << anchorPoint << endl;
|
||||
// Info << "\toriginSpecifies: " << originSpecifies << endl;
|
||||
// Info << "\ttethered: " << tethered << endl;
|
||||
// Info << "\torientationAngles: " << orientationAngles << endl;
|
||||
// Info << "\tlatticeToGlobal: " << latticeToGlobal << endl;
|
||||
|
||||
@ -1,97 +0,0 @@
|
||||
scalar xMax = 0;
|
||||
|
||||
scalar yMax = 0;
|
||||
|
||||
scalar zMax = 0;
|
||||
|
||||
scalar xMin = 0;
|
||||
|
||||
scalar yMin = 0;
|
||||
|
||||
scalar zMin = 0;
|
||||
|
||||
label xMaxPtLabel = 0;
|
||||
|
||||
label yMaxPtLabel = 0;
|
||||
|
||||
label zMaxPtLabel = 0;
|
||||
|
||||
label xMinPtLabel = 0;
|
||||
|
||||
label yMinPtLabel = 0;
|
||||
|
||||
label zMinPtLabel = 0;
|
||||
|
||||
forAll (cellZoneI[cZ], nC)
|
||||
{
|
||||
const labelList& cellPointsJ = mesh_.cellPoints()[cellZoneI[cZ][nC]];
|
||||
|
||||
forAll(cellPointsJ, nP)
|
||||
{
|
||||
const point& ptI = mesh_.points()[cellPointsJ[nP]];
|
||||
|
||||
const label& ptILabel = cellPointsJ[nP];
|
||||
|
||||
if (ptI.x() > xMax || nC == 0)
|
||||
{
|
||||
xMax = ptI.x();
|
||||
xMaxPtLabel = ptILabel;
|
||||
}
|
||||
if (ptI.y() > yMax || nC == 0)
|
||||
{
|
||||
yMax = ptI.y();
|
||||
yMaxPtLabel = ptILabel;
|
||||
}
|
||||
if (ptI.z() > zMax || nC == 0)
|
||||
{
|
||||
zMax = ptI.z();
|
||||
zMaxPtLabel = ptILabel;
|
||||
}
|
||||
if (ptI.x() < xMin || nC == 0)
|
||||
{
|
||||
xMin = ptI.x();
|
||||
xMinPtLabel = ptILabel;
|
||||
}
|
||||
if (ptI.y() < yMin || nC == 0)
|
||||
{
|
||||
yMin = ptI.y();
|
||||
yMinPtLabel = ptILabel;
|
||||
}
|
||||
if (ptI.z() < zMin || nC == 0)
|
||||
{
|
||||
zMin = ptI.z();
|
||||
zMinPtLabel = ptILabel;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Info << "Xmax: label = " << xMaxPtLabel2 << "; vector = " <<mesh_.points()[xMaxPtLabel2]
|
||||
// <<"; x-component = " << mesh_.points()[xMaxPtLabel2].x() << endl;
|
||||
// Info << "Ymax: label = " << yMaxPtLabel2 << "; vector = " <<mesh_.points()[yMaxPtLabel2]
|
||||
// <<"; y-component = " << mesh_.points()[yMaxPtLabel2].y() << endl;
|
||||
// Info << "Zmax: label = " << zMaxPtLabel2 << "; vector = " <<mesh_.points()[zMaxPtLabel2]
|
||||
// <<"; z-component = " << mesh_.points()[zMaxPtLabel2].z() << endl;
|
||||
//
|
||||
// Info << "Xmin: label = " << xMinPtLabel << "; vector = " <<mesh_.points()[xMinPtLabel]
|
||||
// <<"; x-component = " << mesh_.points()[xMinPtLabel].x() << endl;
|
||||
// Info << "Ymin: label = " << yMinPtLabel << "; vector = " <<mesh_.points()[yMinPtLabel]
|
||||
// <<"; y-component = " << mesh_.points()[yMinPtLabel].y() << endl;
|
||||
// Info << "Zmin: label = " << zMinPtLabel << "; vector = " <<mesh_.points()[zMinPtLabel]
|
||||
// <<"; z-component = " << mesh_.points()[zMinPtLabel].z() << endl;
|
||||
|
||||
scalar xMid =
|
||||
(mesh_.points()[xMaxPtLabel].x()
|
||||
+ mesh_.points()[xMinPtLabel].x()) / 2;
|
||||
|
||||
scalar yMid =
|
||||
(mesh_.points()[yMaxPtLabel].y()
|
||||
+ mesh_.points()[yMinPtLabel].y()) / 2;
|
||||
|
||||
scalar zMid =
|
||||
(mesh_.points()[zMaxPtLabel].z()
|
||||
+ mesh_.points()[zMinPtLabel].z()) / 2;
|
||||
|
||||
vector rS(xMid, yMid, zMid);
|
||||
|
||||
// Info << "\t The Estimated Starting Point: " << rS << endl;
|
||||
|
||||
@ -1,26 +0,0 @@
|
||||
scalar velCmptMag = sqrt(moleculeCloud::kb*temperature/mass);
|
||||
|
||||
for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
|
||||
{
|
||||
// Assign velocity: random direction, magnitude determined by desired
|
||||
// maxwellian distribution at temperature
|
||||
|
||||
// Temperature gradients could be created by specifying a gradient in the
|
||||
// zone subDict, or by reading a field from a mesh.
|
||||
|
||||
// The velocities are treated on a zone-by-zone basis for the purposes of
|
||||
// removal of bulk momentum - hence nMols becomes totalZoneMols
|
||||
|
||||
velocity = vector
|
||||
(
|
||||
velCmptMag*rand.GaussNormal(),
|
||||
velCmptMag*rand.GaussNormal(),
|
||||
velCmptMag*rand.GaussNormal()
|
||||
);
|
||||
|
||||
momentumSum += mass*velocity;
|
||||
|
||||
initialVelocities.append(velocity);
|
||||
}
|
||||
|
||||
|
||||
@ -1,27 +0,0 @@
|
||||
scalar initVelMag =
|
||||
sqrt
|
||||
(
|
||||
3.0*(1.0 - 1.0 / totalZoneMols)
|
||||
*moleculeCloud::kb*temperature
|
||||
/mass
|
||||
);
|
||||
|
||||
for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
|
||||
{
|
||||
// Assign velocity: random direction, magnitude determined by desired
|
||||
// temperature
|
||||
|
||||
// Temperature gradients could be created by specifying a gradient in the
|
||||
// zone subDict, or by reading a field from a mesh.
|
||||
|
||||
// The velocities are treated on a zone-by-zone basis for the purposes of
|
||||
// removal of bulk momentum - hence nMols becomes totalZoneMols
|
||||
|
||||
velocity = (2.0*rand.vector01() - vector::one);
|
||||
|
||||
velocity *= initVelMag/mag(velocity);
|
||||
|
||||
momentumSum += mass*velocity;
|
||||
|
||||
initialVelocities.append(velocity);
|
||||
}
|
||||
@ -1,73 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.3 |
|
||||
| \\ / A nd | Web: http://www.openfoam.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
|
||||
root "";
|
||||
case "";
|
||||
instance "";
|
||||
local "";
|
||||
|
||||
class dictionary;
|
||||
object blockMeshDict;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
convertToMeters 2.462491658e-9;
|
||||
|
||||
vertices
|
||||
(
|
||||
(-1 -1 -1)
|
||||
(1 -1 -1)
|
||||
(1 1 -1)
|
||||
(-1 1 -1)
|
||||
(-1 -1 1)
|
||||
(1 -1 1)
|
||||
(1 1 1)
|
||||
(-1 1 1)
|
||||
);
|
||||
|
||||
blocks
|
||||
(
|
||||
hex (0 1 2 3 4 5 6 7) liquid (12 12 12) simpleGrading (1 1 1)
|
||||
);
|
||||
|
||||
patches
|
||||
(
|
||||
cyclic
|
||||
periodicX
|
||||
(
|
||||
(1 2 6 5)
|
||||
(0 4 7 3)
|
||||
)
|
||||
|
||||
cyclic
|
||||
periodicY
|
||||
(
|
||||
(2 3 7 6)
|
||||
(0 1 5 4)
|
||||
)
|
||||
|
||||
cyclic
|
||||
periodicZ
|
||||
(
|
||||
(0 3 2 1)
|
||||
(4 5 6 7)
|
||||
)
|
||||
)
|
||||
|
||||
mergePatchPairs
|
||||
(
|
||||
);
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,67 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.3 |
|
||||
| \\ / A nd | Web: http://www.openfoam.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
|
||||
root "";
|
||||
case "";
|
||||
instance "";
|
||||
local "";
|
||||
|
||||
class dictionary;
|
||||
object fvSchemes;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
ddtSchemes
|
||||
{
|
||||
default Euler;
|
||||
}
|
||||
|
||||
gradSchemes
|
||||
{
|
||||
default Gauss linear;
|
||||
grad(p) Gauss linear;
|
||||
}
|
||||
|
||||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
{
|
||||
default none;
|
||||
laplacian(nu,U) Gauss linear corrected;
|
||||
laplacian(1|A(U),p) Gauss linear corrected;
|
||||
}
|
||||
|
||||
interpolationSchemes
|
||||
{
|
||||
default linear;
|
||||
interpolate(HbyA) linear;
|
||||
}
|
||||
|
||||
snGradSchemes
|
||||
{
|
||||
default corrected;
|
||||
}
|
||||
|
||||
fluxRequired
|
||||
{
|
||||
default no;
|
||||
p;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,40 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.3 |
|
||||
| \\ / A nd | Web: http://www.openfoam.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
|
||||
root "";
|
||||
case "";
|
||||
instance "";
|
||||
local "";
|
||||
|
||||
class dictionary;
|
||||
object fvSolution;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
solvers
|
||||
{
|
||||
p ICCG 1e-06 0;
|
||||
U BICCG 1e-05 0;
|
||||
}
|
||||
|
||||
PISO
|
||||
{
|
||||
nCorrectors 2;
|
||||
nNonOrthogonalCorrectors 0;
|
||||
pRefCell 0;
|
||||
pRefValue 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,27 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.4.1 |
|
||||
| \\ / A nd | Web: http://www.openfoam.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
|
||||
root "";
|
||||
case "";
|
||||
instance "";
|
||||
local "";
|
||||
|
||||
class dictionary;
|
||||
object mdEquilibrationDict;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
equilibrationTargetTemperature 300.0;
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,43 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.3 |
|
||||
| \\ / A nd | Web: http://www.openfoam.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
|
||||
root "";
|
||||
case "";
|
||||
instance "";
|
||||
local "";
|
||||
|
||||
class dictionary;
|
||||
object molConfigDict;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
// Euler angles, expressed in degrees as phi, theta, psi,
|
||||
// see http://mathworld.wolfram.com/EulerAngles.html
|
||||
|
||||
liquid
|
||||
{
|
||||
massDensity 1220.0;
|
||||
temperature 300.0;
|
||||
velocityDistribution maxwellian;
|
||||
bulkVelocity (0.0 0.0 0.0);
|
||||
id Ar;
|
||||
mass 6.63352033e-26;
|
||||
latticeStructure SC;
|
||||
anchor (0.0 0.0 0.0);
|
||||
anchorSpecifies molecule;
|
||||
tethered no;
|
||||
orientationAngles (0 0 0);
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user