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

Functions from old moleculeCloud transferred in and restructured.

Browse Source
This commit is contained in:
graham
2008-09-21 15:45:09 +01:00
parent 5d51971a49
commit 8138fcf8d7
7 changed files with 846 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
View File

@ -185,9 +185,11 @@ public:
inline const referredCellList& referredInteractionList() const;
inline const List<sendingReferralList>& cellSendingReferralLists() const;
inline const List<sendingReferralList>&
cellSendingReferralLists() const;
inline const List<receivingReferralList>& cellReceivingReferralLists() const;
inline const List<receivingReferralList>&
cellReceivingReferralLists() const;
inline label nInteractingProcs() const;
};

3
src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
View File

@ -40,7 +40,8 @@ const Foam::directInteractionList& Foam::interactionLists::dil() const
}
inline const Foam::referredCellList& Foam::interactionLists::referredInteractionList() const
inline const Foam::referredCellList&
Foam::interactionLists::referredInteractionList() const
{
return ril_;
}

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

@ -85,10 +85,11 @@ bool Foam::molecule::move(molecule::trackData& td)
if (td.part() == 0)
{
// Leapfrog velocity adjust part, required before and after
// tracking+force part
// First leapfrog velocity adjust part, required before tracking+force part
U_ += 0.5*deltaT*A_;
v_ += 0.5*deltaT*a_;
omega1 += 0.5*deltaT*alpha_;
}
else if (td.part() == 1)
{
@ -99,7 +100,7 @@ bool Foam::molecule::move(molecule::trackData& td)
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
dt *= trackToFace(position() + dt*U_, td);
dt *= trackToFace(position() + dt*v_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/deltaT;
@ -115,6 +116,33 @@ bool Foam::molecule::move(molecule::trackData& td)
setSitePositions(td.molCloud().constProps(id_));
}
else if (td.part() == 3)
{
// Second leapfrog velocity adjust part, required after tracking+force part
const diagTensor& I; // TAKE REFERENCE TO CONST PROP TO GET THIS
scalar m; // TAKE REFERENCE TO CONST PROP TO GET THIS
a_ = vector::zero;
vector tau(vector::zero);
forAll(siteForces_, s)
{
const vector& f = siteForces_[s];
a_ += f/m;
tau += (sitePositions_[s] - position_) ^ f;
}
alpha_ = R_ & inv(I) & R_.T() & tau;
v_ += 0.5*deltaT*a_;
omega_ += 0.5*deltaT*alpha_;
}
else
{
FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
@ -173,8 +201,8 @@ void Foam::molecule::hitWallPatch
vector nw = wpp.faceAreas()[wpp.whichFace(face())];
nw /= mag(nw);
scalar Un = U_ & nw;
// vector Ut = U_ - Un*nw;
scalar vn = v_ & nw;
// vector vt = v_ - vn*nw;
// Random rand(clock::getTime());
@ -186,11 +214,11 @@ void Foam::molecule::hitWallPatch
// {
// // Diffuse reflection
//
// vector tw1 = Ut/mag(Ut);
// vector tw1 = vt/mag(vt);
//
// vector tw2 = nw ^ tw1;
//
// U_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
// V_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
// + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2
// - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw;
// }
@ -199,9 +227,9 @@ void Foam::molecule::hitWallPatch
// {
// Specular reflection
if (Un > 0)
if (vn > 0)
{
U_ -= 2*Un*nw;
v_ -= 2*vn*nw;
}
// }

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

@ -42,6 +42,401 @@ Foam::scalar Foam::moleculeCloud::elementaryCharge = 1.602176487e-19;
Foam::scalar Foam::moleculeCloud::vacuumPermittivity = 8.854187817e-12;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::moleculeCloud::buildCellOccupancy()
{
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].clear();
}
iterator mol(this->begin());
for
(
mol = this->begin();
mol != this->end();
++mol
)
{
cellOccupancy_[mol().cell()].append(&mol());
}
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].shrink();
}
il_ril().referMolecules(cellOccupancy_);
}
void Foam::moleculeCloud::calculatePairForce()
{
{
iterator mol(this->begin());
// Real-Real interactions
vector rIJ;
scalar rIJMag;
scalar rIJMagSq;
vector fIJ;
label idI;
label idJ;
mol = this->begin();
molecule* molI = &mol();
molecule* molJ = &mol();
forAll(directInteractionList_, dIL)
{
forAll(cellOccupancy_[dIL],cellIMols)
{
molI = cellOccupancy_[dIL][cellIMols];
forAll(directInteractionList_[dIL], interactingCells)
{
List< molecule* > cellJ =
cellOccupancy_[directInteractionList_[dIL][interactingCells]];
forAll(cellJ, cellJMols)
{
molJ = cellJ[cellJMols];
evaluatelPair(molI, molJ);
}
}
forAll(cellOccupancy_[dIL],cellIOtherMols)
{
molJ = cellOccupancy_[dIL][cellIOtherMols];
if (molJ > molI)
{
evaluatePair(molI, molJ);
}
}
}
}
}
{
// Real-Referred interactions
vector rKL;
scalar rKLMag;
scalar rKLMagSq;
vector fKL;
label idK;
label idL;
molecule* molK = &mol();
forAll(referredInteractionList_, rIL)
{
const List<label>& realCells =
referredInteractionList_[rIL].realCellsForInteraction();
forAll(referredInteractionList_[rIL], refMols)
{
referredMolecule* molL = &(referredInteractionList_[rIL][refMols]);
forAll(realCells, rC)
{
List<molecule*> cellK = cellOccupancy_[realCells[rC]];
forAll(cellK, cellKMols)
{
molK = cellK[cellKMols];
evaluatePair(molK, molL);
}
}
}
}
}
}
void Foam::moleculeCloud::calculateTetherForce()
{
iterator mol(this->begin());
vector rIT;
vector fIT;
for (mol = this->begin(); mol != this->end(); ++mol)
{
if (mol().tethered())
{
rIT = mol().position() - mol().tetherPosition();
fIT = tetherPotentials_.force
(
mol().id(),
rIT
);
mol().A() += fIT/(mol().mass());
mol().potentialEnergy() += tetherPotentials_.energy
(
mol().id(),
rIT
);
mol().rf() += rIT*fIT;
}
}
}
void Foam::moleculeCloud::calculateExternalForce()
{
iterator mol(this->begin());
// Info<< "Warning! Includes dissipation term!" << endl;
for (mol = this->begin(); mol != this->end(); ++mol)
{
mol().A() += gravity_;
// mol().A() += -1.0 * mol().U() /mol().mass();
}
}
void Foam::moleculeCloud::removeHighEnergyOverlaps()
{
Info << nl << "Removing high energy overlaps, removal order:";
forAll(removalOrder_, rO)
{
Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
}
Info << nl ;
label initialSize = this->size();
buildCellOccupancy();
// Real-Real interaction
iterator mol(this->begin());
{
vector rIJ;
scalar rIJMag;
scalar rIJMagSq;
label idI;
label idJ;
mol = this->begin();
molecule* molI = &mol();
molecule* molJ = &mol();
DynamicList<molecule*> molsToDelete;
forAll(directInteractionList_, dIL)
{
forAll(cellOccupancy_[dIL],cellIMols)
{
molI = cellOccupancy_[dIL][cellIMols];
forAll(directInteractionList_[dIL], interactingCells)
{
List< molecule* > cellJ =
cellOccupancy_[directInteractionList_[dIL][interactingCells]];
forAll(cellJ, cellJMols)
{
molJ = cellJ[cellJMols];
if (evaluatePotentialLimit(molI, molJ))
{
if
(
idI == idJ
|| findIndex(removalOrder_, idJ)
< findIndex(removalOrder_, idI)
)
{
if (findIndex(molsToDelete, molJ) == -1)
{
molsToDelete.append(molJ);
}
}
else if (findIndex(molsToDelete, molI) == -1)
{
molsToDelete.append(molI);
}
}
}
}
}
forAll(cellOccupancy_[dIL],cellIOtherMols)
{
molJ = cellOccupancy_[dIL][cellIOtherMols];
if (molJ > molI)
{
if (evaluatePotentialLimit(molI, molJ))
{
if
(
idI == idJ
|| findIndex(removalOrder_, idJ)
< findIndex(removalOrder_, idI)
)
{
if (findIndex(molsToDelete, molJ) == -1)
{
molsToDelete.append(molJ);
}
}
else if (findIndex(molsToDelete, molI) == -1)
{
molsToDelete.append(molI);
}
}
}
}
}
forAll (molsToDelete, mTD)
{
deleteParticle(*(molsToDelete[mTD]));
}
}
buildCellOccupancy();
// Real-Referred interaction
{
vector rKL;
scalar rKLMag;
scalar rKLMagSq;
label idK;
label idL;
molecule* molK = &mol();
DynamicList<molecule*> molsToDelete;
forAll(referredInteractionList_, rIL)
{
referredCell& refCell = referredInteractionList_[rIL];
forAll(refCell, refMols)
{
referredMolecule* molL = &(refCell[refMols]);
List <label> realCells = refCell.realCellsForInteraction();
forAll(realCells, rC)
{
label cellK = realCells[rC];
List<molecule*> cellKMols = cellOccupancy_[cellK];
forAll(cellKMols, cKM)
{
molK = cellKMols[cKM];
if (evaluatePotentialLimit(molK, molL))
{
if
(
findIndex(removalOrder_, idK)
< findIndex(removalOrder_, idL)
)
{
if (findIndex(molsToDelete, molK) == -1)
{
molsToDelete.append(molK);
}
}
else if
(
findIndex(removalOrder_, idK)
== findIndex(removalOrder_, idL)
)
{
if
(
Pstream::myProcNo() == refCell.sourceProc()
&& cellK <= refCell.sourceCell()
)
{
if (findIndex(molsToDelete, molK) == -1)
{
molsToDelete.append(molK);
}
}
else if
(
Pstream::myProcNo() < refCell.sourceProc()
)
{
if (findIndex(molsToDelete, molK) == -1)
{
molsToDelete.append(molK);
}
}
}
}
}
}
}
}
forAll (molsToDelete, mTD)
{
deleteParticle(*(molsToDelete[mTD]));
}
}
buildCellOccupancy();
label molsRemoved = initialSize - this->size();
if (Pstream::parRun())
{
reduce(molsRemoved, sumOp<label>());
}
Info << tab << molsRemoved << " molecules removed" << endl;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
@ -115,48 +510,74 @@ Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::moleculeCloud::buildCellOccupancy()
{
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].clear();
}
iterator mol(this->begin());
for
(
mol = this->begin();
mol != this->end();
++mol
)
{
cellOccupancy_[mol().cell()].append(&mol());
}
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].shrink();
}
il_ril().referMolecules(cellOccupancy_);
}
void Foam::moleculeCloud::evolve()
{
molecule::trackData td0(*this, 0);
molecule::trackData td1(*this, 1);
molecule::trackData td2(*this, 2);
Cloud<molecule>::move(td0);
molecule::trackData td1(*this, 1);
Cloud<molecule>::move(td1);
molecule::trackData td2(*this, 2);
Cloud<molecule>::move(td2);
calculateForce();
Cloud<molecule>::move(td1);
molecule::trackData td2(*this, 3);
Cloud<molecule>::move(td3);
}
Cloud<molecule>::move(td2);
Cloud<molecule>::move(td0);
void Foam::moleculeCloud::calculateForce()
{
buildCellOccupancy();
iterator mol(this->begin());
// Set accumulated quantities to zero
for (mol = this->begin(); mol != this->end(); ++mol)
{
mol().siteForces() = vector::zero;
mol().potentialEnergy() = 0.0;
mol().rf() = tensor::zero;
}
calculatePairForce();
calculateTetherForce();
calculateExternalForce();
}
void Foam::moleculeCloud::applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
)
{
scalar temperatureCorrectionFactor =
sqrt(targetTemperature/measuredTemperature);
Info<< "----------------------------------------" << nl
<< "Temperature equilibration" << nl
<< "Target temperature = "
<< targetTemperature << nl
<< "Measured temperature = "
<< measuredTemperature << nl
<< "Temperature correction factor = "
<< temperatureCorrectionFactor << nl
<< "----------------------------------------"
<< endl;
iterator mol(this->begin());
for (mol = this->begin(); mol != this->end(); ++mol)
{
mol().u() *= temperatureCorrectionFactor;
}
}

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

@ -63,13 +63,13 @@ private:
const polyMesh& mesh_;
// MD solution parameters
// MD solution parameters
List<DynamicList<molecule*> > cellOccupancy_;
List<DynamicList<molecule*> > cellOccupancy_;
interactionLists il_;
interactionLists il_;
List<molecule::constantProperties> constPropList_;;
List<molecule::constantProperties> constPropList_;
// Move to potential class and store reference to it ->
scalar potentialEnergyLimit_;
@ -83,29 +83,55 @@ private:
vector gravity_;
// <- Move to potential class and store reference to it
// Private Member Functions
// Private Member Functions
void applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
);
//- Determine which molecules are in which cells
void buildCellOccupancy();
void calculateForce();
void calculatePairForce();
void calculatePairForce();
inline void evaluatePair
(
molecule* molI,
molecule* molJ
);
void calculateTetherForce();
inline void evaluatePair
(
molecule* molReal,
referredMolecule* molRef
)
void calculateExternalForce();
inline bool evaluatePotentialLimit
(
label idI,
label idJ,
scalar rIJMagSq
) const;
void removeHighEnergyOverlaps();
inline bool evaluatePotentialLimit
(
molecule* molI,
molecule* molJ
) const;
//- Disallow default bitwise copy construct
moleculeCloud(const moleculeCloud&);
inline bool evaluatePotentialLimit
(
molecule* molReal,
referredMolecule* molRef
) const;
//- Disallow default bitwise assignment
void operator=(const moleculeCloud&);
void calculateTetherForce();
void calculateExternalForce();
void removeHighEnergyOverlaps();
//- Disallow default bitwise copy construct
moleculeCloud(const moleculeCloud&);
//- Disallow default bitwise assignment
void operator=(const moleculeCloud&);
public:
@ -142,22 +168,31 @@ public:
// Member Functions
//- Determine which molecules are in which cells
void buildCellOccupancy();
//- Evolve the molecules (move, calculate forces, control state etc)
void evolve();
//- Integrate the equations of motion using algorithm selected at
// runtime from a dictionary. This will also call the function
// to calculate the intermolecular forces (calculatePairForce()).
void calculateForce();
void applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
);
// Access
inline const polyMesh& mesh() const;
inline const interactionLists& interactionLists() const;
inline const List<DynamicList<molecule*> >& cellOccupancy() const;
inline const List<molecule::constantProperties> constProps() const;
inline const molecule::constantProperties& constProps(label id) const;
// Move to potential class and store reference to it ->
inline scalar potentialEnergyLimit() const;
inline const labelList& removalOrder() const;
@ -167,10 +202,7 @@ public:
inline const tetherPotentialList& tetherPotentials() const;
inline const vector& gravity() const;
inline const List<DynamicList<molecule*> >& cellOccupancy() const;
// <- Move to potential class and store reference to it
// Member Operators

285
src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H Normal file
View File

@ -0,0 +1,285 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::moleculeCloud::evaluatePair
(
molecule* molI,
molecule* molJ
)
{
idI = molI->id();
idJ = molJ->id();
rIJ = molI->position() - molJ->position();
rIJMagSq = magSqr(rIJ);
if(pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
{
rIJMag = mag(rIJ);
fIJ = (rIJ/rIJMag)*pairPotentials_.force(idI, idJ, rIJMag);
// Acceleration increment for molI
molI->A() += fIJ/(molI->mass());
// Acceleration increment for molJ
molJ->A() += -fIJ/(molJ->mass());
scalar potentialEnergy
(
pairPotentials_.energy(idI, idJ, rIJMag)
);
molI->potentialEnergy() += 0.5*potentialEnergy;
molJ->potentialEnergy() += 0.5*potentialEnergy;
molI->rf() += rIJ * fIJ;
molJ->rf() += rIJ * fIJ;
}
}
inline void Foam::moleculeCloud::evaluatePair
(
molecule* molReal,
referredMolecule* molRef
)
{
idReal = molReal->id();
idRef = molRef->id();
rRealRef = molReal->position() - molRef->position();
rRealRefMagSq = magSqr(rRealRef);
if (pairPotentials_.rCutSqr(idReal, idRef, rRealRefMagSq))
{
rRealRefMag = mag(rRealRef);
fRealRef = (rRealRef/rRealRefMag)
*pairPotentials_.force(idReal, idRef, rRealRefMag);
// Acceleration increment for molReal
molReal->A() += fRealRef/(molReal->mass());
// Adding a contribution of 1/2 of the potential energy
// from this interaction
molReal->potentialEnergy() +=
0.5*pairPotentials_.energy(idReal, idRef, rRealRefMag);
molReal->rf() += rRealRef * fRealRef;
}
}
inline bool Foam::moleculeCloud::evaluatePotentialLimit
(
label idI,
label idJ,
scalar rIJMagSq
) const
{
Info << "Move Foam::moleculeCloud::evaluatePotentialLimit "
<< "(label idI, label idJ, scalar rIJMagSq) const "
<< "to pairPotential class."
bool remove = false;
if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
{
rIJMag = mag(rIJ);
// Guard against pairPotentials_.energy being evaluated
// if rIJMag < SMALL. A floating point exception will
// happen otherwise.
if (rIJMag < SMALL)
{
WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
<< "Real molecule pair "
<< " idI = " << idI
<< " at position " << molI->position()
<< " idJ = " << idJ
<< " at position " << molJ->position()
<< " are closer than " << SMALL
<< ": mag separation = " << rIJMag
<< ". These may have been placed on top of each"
<< " other by a rounding error in molConfig in parallel"
<< " or a block filled with molecules twice."
<< " Removing one of the molecules."
<< endl;
remove = true;
}
// Guard against pairPotentials_.energy being
// evaluated if rIJMag < rMin. A tabulation lookup
// error will occur otherwise.
if (rIJMag < pairPotentials_.rMin(idI, idJ))
{
remove = true;
}
if (!remove)
{
if
(
pairPotentials_.energy(idI, idJ, rIJMag)
> potentialEnergyLimit_
)
{
remove = true;
}
}
}
return remove;
}
inline bool Foam::moleculeCloud::evaluatePotentialLimit
(
molecule* molI,
molecule* molJ
) const
{
idI = molI->id();
idJ = molJ->id();
rIJ = molI->position() - molJ->position();
rIJMagSq = magSqr(rIJ);
return evaluatePotentialLimit(idI, idJ, rIJMagSq);
}
inline bool Foam::moleculeCloud::evaluatePotentialLimit
(
molecule* molReal,
referredMolecule* molRef
) const
{
idReal = molReal->id();
idRef = molRef->id();
rRealRef = molReal->position() - molRef->position();
rRealRefMagSq = magSqr(rRealRef);
return evaluatePotentialRefimit(idReal, idRef, rRealRefMagSq);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::polyMesh& moleculeCloud::mesh() const
{
return mesh_;
}
inline const Foam::interactionLists&
Foam::moleculeCloud::interactionLists() const
{
return il_;
}
inline const Foam::List<DynamicList<molecule*> >&
Foam::moleculeCloud::cellOccupancy() const
{
return cellOccupancy_;
}
inline const Foam::List<molecule::constantProperties> constProps() const
{
return constPropList_;
}
inline const molecule::constantProperties& constProps(label id) const
{
return constPropList_[id];
}
// Move to potential class and store reference to it ->
inline Foam::scalar Foam::moleculeCloud::potentialEnergyLimit() const
{
return potentialEnergyLimit_;
}
inline Foam::label Foam::moleculeCloud::nPairPotentials() const
{
return pairPotentials_.size();
}
inline const Foam::labelList& Foam::moleculeCloud::removalOrder() const
{
return removalOrder_;
}
inline const Foam::pairPotentialList&
Foam::moleculeCloud::pairPotentials() const
{
return pairPotentials_;
}
inline const Foam::tetherPotentialList&
Foam::moleculeCloud::tetherPotentials() const
{
return tetherPotentials_;
}
inline const Foam::vector& Foam::moleculeCloud::gravity() const
{
return gravity_;
}
// <- Move to potential class and store reference to it
// ************************************************************************* //

4
src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
View File

@ -57,7 +57,7 @@ class pairPotentialList
// Private data
List<word> idList_;
scalar rCutMax_;
scalar rCutMaxSqr_;
@ -101,7 +101,7 @@ public:
~pairPotentialList();
// Member Functions
void buildPotentials
(
const dictionary& idListDict,