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

Added momentum and energy monitors, angular momentum and energy running away - need to investigate.

Browse Source
This commit is contained in:
graham
2008-11-24 16:57:17 +00:00
parent 0e78f04bd8
commit c54c92dabb
5 changed files with 72 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/solvers/molecularDynamics/mdFoam/mdFoam.C
View File

@ -45,6 +45,8 @@ int main(int argc, char *argv[])
moleculeCloud molecules(mesh, pot);
# include "temperatureAndPressureVariables.H"
Info << "\nStarting time loop\n" << endl;
while (runTime.run())
@ -55,6 +57,8 @@ int main(int argc, char *argv[])
molecules.evolve();
# include "meanMomentumEnergyAndNMols.H"
runTime.write();
Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

1
src/lagrangian/molecularDynamics/Allwmake
View File

@ -3,6 +3,7 @@ cd ${0%/*} || exit 1 # run from this directory
set -x
wmake libso potential
wmake libso molecularMeasurements
wmake libso molecule
# ----------------------------------------------------------------- end-of-file

72
src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
View File

@ -32,13 +32,17 @@ Description
\*---------------------------------------------------------------------------*/
vector singleStepTotalMomentum(vector::zero);
vector singleStepTotalLinearMomentum(vector::zero);
vector singleStepTotalAngularMomentum(vector::zero);
scalar singleStepMaxVelocityMag = 0.0;
scalar singleStepTotalMass = 0.0;
scalar singleStepTotalKE = 0.0;
scalar singleStepTotalLinearKE = 0.0;
scalar singleStepTotalAngularKE = 0.0;
scalar singleStepTotalPE = 0.0;
@ -54,20 +58,30 @@ scalar singleStepTotalrDotf = 0.0;
++mol
)
{
const scalar molM(mol().mass());
label molId = mol().id();
const vector& molU(mol().U());
scalar molMass(molecules.constProps(molId).mass());
singleStepTotalMomentum += molU * molM;
const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
singleStepTotalMass += molM;
const vector& molV(mol().v());
if(mag(molU) > singleStepMaxVelocityMag)
const vector& molOmega(mol().omega());
singleStepTotalLinearMomentum += molV * molMass;
singleStepTotalAngularMomentum += (molMoI & molOmega);
singleStepTotalMass += molMass;
if(mag(molV) > singleStepMaxVelocityMag)
{
singleStepMaxVelocityMag = mag(molU);
singleStepMaxVelocityMag = mag(molV);
}
singleStepTotalKE += 0.5*molM*magSqr(molU);
singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
singleStepTotalPE += mol().potentialEnergy();
@ -79,13 +93,17 @@ label singleStepNMols = molecules.size();
if (Pstream::parRun())
{
reduce(singleStepTotalMomentum, sumOp<vector>());
reduce(singleStepTotalLinearMomentum, sumOp<vector>());
reduce(singleStepTotalAngularMomentum, sumOp<vector>());
reduce(singleStepMaxVelocityMag, maxOp<scalar>());
reduce(singleStepTotalMass, sumOp<scalar>());
reduce(singleStepTotalKE, sumOp<scalar>());
reduce(singleStepTotalLinearKE, sumOp<scalar>());
reduce(singleStepTotalAngularKE, sumOp<scalar>());
reduce(singleStepTotalPE, sumOp<scalar>());
@ -96,22 +114,32 @@ if (Pstream::parRun())
if (singleStepNMols)
{
Info << "Number of mols in system = "
<< singleStepNMols << nl
Info<< "Number of mols in system = "
<< singleStepNMols << nl
<< "Overall number density = "
<< singleStepNMols/meshVolume << " m^-3" << nl
<< singleStepNMols/meshVolume << nl
<< "Overall mass density = "
<< singleStepTotalMass/meshVolume << " kg/m^3" << nl
<< "Average velocity per mol = "
<< singleStepTotalMomentum/singleStepTotalMass << " m/s" << nl
<< singleStepTotalMass/meshVolume << nl
<< "Average linear momentum per mol = "
<< singleStepTotalLinearMomentum/singleStepNMols << nl
<< "Average angular momentum per mol = "
<< singleStepTotalAngularMomentum/singleStepNMols << nl
<< "Maximum |velocity| = "
<< singleStepMaxVelocityMag << " m/s" << nl
<< "Average KE per mol = "
<< singleStepTotalKE/singleStepNMols << " J" << nl
<< singleStepMaxVelocityMag << nl
<< "Average linear KE per mol = "
<< singleStepTotalLinearKE/singleStepNMols << nl
<< "Average angular KE per mol = "
<< singleStepTotalAngularKE/singleStepNMols << nl
<< "Average PE per mol = "
<< singleStepTotalPE/singleStepNMols << " J" << nl
<< singleStepTotalPE/singleStepNMols << nl
<< "Average TE per mol = "
<< (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << " J"
<<
(
singleStepTotalLinearKE
+ singleStepTotalAngularKE
+ singleStepTotalPE
)
/singleStepNMols
<< endl;
// Info << singleStepNMols << " "

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

@ -250,14 +250,16 @@ void Foam::moleculeCloud::calculateExternalForce()
void Foam::moleculeCloud::removeHighEnergyOverlaps()
{
Info << nl << "Removing high energy overlaps, removal order:";
Info<< nl << "Removing high energy overlaps, limit = "
<< pot_.potentialEnergyLimit()
<< nl << "Removal order:";
forAll(pot_.removalOrder(), rO)
{
Info << ' ' << pot_.idList()[pot_.removalOrder()[rO]];
Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
}
Info << nl ;
Info<< nl ;
label initialSize = this->size();
@ -457,7 +459,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
reduce(molsRemoved, sumOp<label>());
}
Info << tab << molsRemoved << " molecules removed" << endl;
Info<< tab << molsRemoved << " molecules removed" << endl;
}
@ -483,17 +485,19 @@ Foam::moleculeCloud::moleculeCloud
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const Cloud<molecule>& cloud = *this;
vector position1 = vector(0,0,1.5e-9);
addParticle
(
new molecule
(
cloud,
vector(0,0,0.5e-9),
13,
position1,
mesh_.findCell(position1),
tensor(1, 0, 0, 0, 1, 0, 0, 0, 1),
vector(234, 102, -365),
vector(34, 12, -65),
vector::zero,
vector(2e10, -3e11, 1e10),
vector(0.2, -0.012, 0.32),
vector::zero,
vector::zero,
constProps(0),
@ -502,17 +506,19 @@ Foam::moleculeCloud::moleculeCloud
)
);
vector position2 = vector(0,0,0);
addParticle
(
new molecule
(
cloud,
vector(0,0,0),
13,
position2,
mesh_.findCell(position2),
tensor(1, 0, 0, 0, 1, 0, 0, 0, 1),
vector::zero,
vector::zero,
vector(4e11, 1.5e10, 2e-10),
vector(-0.4, 0.015, -0.2),
vector::zero,
vector::zero,
constProps(1),

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

@ -66,8 +66,6 @@ inline void Foam::moleculeCloud::evaluatePair
if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
{
Info<< "Pair potential " << idsI << ' ' << idsJ << endl;
vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
@ -103,8 +101,6 @@ inline void Foam::moleculeCloud::evaluatePair
if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
{
Info<< "Electrostatic " << idsI << ' ' << idsJ << endl;
vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
@ -178,10 +174,6 @@ inline void Foam::moleculeCloud::evaluatePair
if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
{
Info<< "Ref Pair potential "
<< idsReal << ' '
<< idsRef << endl;
vector rsRealsRef =
molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
@ -212,10 +204,6 @@ inline void Foam::moleculeCloud::evaluatePair
if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
{
Info<< "Ref Electrostatic "
<< idsReal << ' '
<< idsRef << endl;
vector rsRealsRef =
molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];