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

Making potential energy limit absolute value. Created all lattice generation variables.

Browse Source
This commit is contained in:
graham
2009-01-06 17:38:48 +00:00
parent e5642a9152
commit 23e1c04cca
4 changed files with 235 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -88,8 +88,6 @@ bool Foam::molecule::move(molecule::trackData& td)
const constantProperties& constProps(td.molCloud().constProps(id_));
const diagTensor& momentOfInertia(constProps.momentOfInertia());
scalar deltaT = cloud().pMesh().time().deltaT().value();
if (td.part() == 0)
@ -125,6 +123,8 @@ bool Foam::molecule::move(molecule::trackData& td)
// but after tracking stage, i.e. rotation carried once linear motion
// complete.
const diagTensor& momentOfInertia(constProps.momentOfInertia());
tensor R;
R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());

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

@ -478,25 +478,21 @@ void Foam::moleculeCloud::initialiseMolecules
<< abort(FatalError);
}
forAll (cellZones, cZ)
forAll (cellZones, z)
{
if (cellZones[cZ].size())
const cellZone& zone(cellZones[z]);
if (zone.size())
{
if (!mdInitialiseDict.found(cellZones[cZ].name()))
if (!mdInitialiseDict.found(zone.name()))
{
Info<< "No specification subDictionary for zone "
<< cellZones[cZ].name() << " found, skipping." << endl;
<< zone.name() << " found, skipping." << endl;
}
else
{
label n = 0;
label totalZoneMols = 0;
label molsPlacedThisIteration = 0;
const dictionary& zoneDict =
mdInitialiseDict.subDict(cellZones[cZ].name());
mdInitialiseDict.subDict(zone.name());
const scalar temperature
(
@ -510,43 +506,88 @@ void Foam::moleculeCloud::initialiseMolecules
const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
// const word id(zoneDict.lookup("id"));
List<word> latticeIds
(
zoneDict.lookup("latticeIds")
);
// const scalar mass(readScalar(zoneDict.lookup("mass")));
List<vector> latticePositions
(
zoneDict.lookup("latticePositions")
);
// scalar numberDensity_read(0.0);
if(latticeIds.size() != latticePositions.size())
{
FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
<< "latticeIds and latticePositions must be the same "
<< " size." << nl
<< abort(FatalError);
}
// if (zoneDict.found("numberDensity"))
// {
// numberDensity_read = readScalar
// (
// zoneDict.lookup("numberDensity")
// );
// }
// else if (zoneDict.found("massDensity"))
// {
// numberDensity_read = readScalar
// (
// zoneDict.lookup("massDensity")
// )/mass;
// }
// else
// {
// FatalErrorIn("readZoneSubDict.H\n")
// << "massDensity or numberDensity not specified " << nl
// << abort(FatalError);
// }
diagTensor latticeCellShape
(
zoneDict.lookup("latticeCellShape")
);
// const scalar numberDensity(numberDensity_read);
scalar latticeCellScale = 0.0;
const vector anchorPoint(zoneDict.lookup("anchor"));
if (zoneDict.found("numberDensity"))
{
scalar numberDensity = readScalar
(
zoneDict.lookup("numberDensity")
);
// bool tethered = false;
latticeCellScale = pow
(
latticeIds.size()/(det(latticeCellShape)*numberDensity),
(1.0/3.0)
);
}
else if (zoneDict.found("massDensity"))
{
scalar unitCellMass = 0.0;
// if (zoneDict.found("tethered"))
// {
// tethered = Switch(zoneDict.lookup("tethered"));
// }
forAll(latticeIds, i)
{
label id = findIndex(pot_.idList(), latticeIds[i]);
const molecule::constantProperties& cP(constProps(id));
unitCellMass += cP.mass();
}
scalar massDensity = readScalar
(
zoneDict.lookup("massDensity")
);
latticeCellScale = pow
(
unitCellMass /(det(latticeCellShape)*massDensity),
(1.0/3.0)
);
}
else
{
FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
<< "massDensity or numberDensity not specified " << nl
<< abort(FatalError);
}
latticeCellShape *= latticeCellScale;
vector anchor(zoneDict.lookup("anchor"));
bool tethered = false;
if (zoneDict.found("tetherSiteIds"))
{
tethered = bool
(
List<word>(zoneDict.lookup("tetherSiteIds")).size()
);
}
const vector orientationAngles
(
@ -568,7 +609,7 @@ void Foam::moleculeCloud::initialiseMolecules
orientationAngles.z()*mathematicalConstant::pi/180.0
);
const tensor latticeToGlobal
const tensor R
(
cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
@ -581,12 +622,149 @@ void Foam::moleculeCloud::initialiseMolecules
cos(theta)
);
// Find the optimal anchor position. Finding the approximate
// mid-point of the zone of cells and snapping to the nearest
// lattice location.
vector zoneMin = VGREAT*vector::one;
vector zoneMax = -VGREAT*vector::one;
forAll (zone, cell)
{
const point cellCentre = mesh_.cellCentres()[zone[cell]];
if (cellCentre.x() > zoneMax.x())
{
zoneMax.x() = cellCentre.x();
}
if (cellCentre.x() < zoneMin.x())
{
zoneMin.x() = cellCentre.x();
}
if (cellCentre.y() > zoneMax.y())
{
zoneMax.y() = cellCentre.y();
}
if (cellCentre.y() < zoneMin.y())
{
zoneMin.y() = cellCentre.y();
}
if (cellCentre.z() > zoneMax.z())
{
zoneMax.z() = cellCentre.z();
}
if (cellCentre.z() < zoneMin.z())
{
zoneMin.z() = cellCentre.z();
}
}
point zoneMid = 0.5*(zoneMin + zoneMax);
point latticeMid = inv(latticeCellShape) & (R.T()
& (zoneMid - anchor));
point latticeAnchor
(
label(latticeMid.x() + 0.5*sign(latticeMid.x())),
label(latticeMid.y() + 0.5*sign(latticeMid.y())),
label(latticeMid.z() + 0.5*sign(latticeMid.z()))
);
anchor += (R & (latticeCellShape & latticeAnchor));
// 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.
label n = 0;
label totalZoneMols = 0;
label molsPlacedThisIteration = 0;
while
(
molsPlacedThisIteration != 0
|| totalZoneMols == 0
)
{
label sizeBeforeIteration = this->size();
bool partOfLayerInBounds = false;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// start of place molecules
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vector unitCellLatticePosition(0,0,0);
// Special treatment is required for the first position,
// i.e. iteration zero.
if (n == 0)
{
forAll(latticePositions, p)
{
const vector& latticePosition =
unitCellLatticePosition + latticePositions[p];
point globalPosition =
anchor + (R & (latticeCellShape & latticePosition));
partOfLayerInBounds =
mesh_.bounds().contains(globalPosition);
label cell = mesh_.findCell(globalPosition);
if (findIndex(zone, cell) != -1)
{
Info<< "add molecule at "
<< globalPosition
<< " to cell "
<< cell
<< endl;
}
}
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// end of place molecules
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 '"
<< zone.name()
<< "'. This is likely to be because the zone "
<< "has few cells ("
<< zone.size()
<< " in this case) and no lattice position "
<< "fell inside them. "
<< "Aborting filling this zone."
<< endl;
break;
}
molsPlacedThisIteration =
this->size() - sizeBeforeIteration;
totalZoneMols += molsPlacedThisIteration;
n++;
}
}
}
}

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

@ -321,7 +321,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if
(
pairPot.energy(idsI, idsJ, rsIsJMag)
mag(pairPot.energy(idsI, idsJ, rsIsJMag))
> pot_.potentialEnergyLimit()
)
{
@ -366,7 +366,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if
(
chargeI*chargeJ*electrostatic.energy(rsIsJMag)
mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
> pot_.potentialEnergyLimit()
)
{
@ -460,7 +460,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if
(
pairPot.energy(idsReal, idsRef, rsRealsRefMag)
mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
> pot_.potentialEnergyLimit()
)
{
@ -505,7 +505,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if
(
chargeReal*chargeRef*electrostatic.energy(rsRealsRefMag)
mag
(
chargeReal*chargeRef
*electrostatic.energy(rsRealsRefMag)
)
> pot_.potentialEnergyLimit()
)
{

10
src/lagrangian/molecularDynamics/potential/potential/potential.C
View File

@ -278,14 +278,14 @@ void Foam::potential::potential::readMdInitialiseDict
mdInitialiseDict.toc()[zone]
);
List<word> ids
List<word> latticeIds
(
zoneDict.lookup("latticeIds")
);
forAll(ids, i)
forAll(latticeIds, i)
{
const word& id = ids[i];
const word& id = latticeIds[i];
if(!moleculePropertiesDict.found(id))
{
@ -314,14 +314,14 @@ void Foam::potential::potential::readMdInitialiseDict
bool idFound = false;
forAll(ids, i)
forAll(latticeIds, i)
{
if (idFound)
{
break;
}
const word& id = ids[i];
const word& id = latticeIds[i];
List<word> siteIds
(