46 void Foam::moleculeCloud::buildConstProps()
48 Info<<
nl <<
"Reading moleculeProperties dictionary." <<
endl;
50 const List<word>& idList(pot_.
idList());
52 constPropList_.
setSize(idList.size());
54 const List<word>& siteIdList(pot_.
siteIdList());
56 IOdictionary moleculePropertiesDict
71 const word& id(idList[i]);
73 const dictionary& molDict(moleculePropertiesDict.subDict(
id));
75 List<word> siteIdNames = molDict.lookup(
"siteIds");
77 List<label> siteIds(siteIdNames.size());
81 const word& siteId = siteIdNames[sI];
83 siteIds[sI] =
findIndex(siteIdList, siteId);
85 if (siteIds[sI] == -1)
88 << siteId <<
" site not found."
93 molecule::constantProperties& constProp = constPropList_[i];
95 constProp = molecule::constantProperties(molDict);
97 constProp.siteIds() = siteIds;
102 void Foam::moleculeCloud::setSiteSizesAndPositions()
104 iterator mol(this->begin());
106 for (mol = this->begin(); mol != this->end(); ++mol)
108 const molecule::constantProperties& cP = constProps(mol().
id());
110 mol().setSiteSizes(cP.nSites());
112 mol().setSitePositions(cP);
117 void Foam::moleculeCloud::buildCellOccupancy()
119 forAll(cellOccupancy_, cO)
121 cellOccupancy_[cO].clear();
124 iterator mol(this->begin());
133 cellOccupancy_[mol().cell()].append(&mol());
136 forAll(cellOccupancy_, cO)
138 cellOccupancy_[cO].shrink();
145 void Foam::moleculeCloud::calculatePairForce()
147 iterator mol(this->begin());
152 molecule* molI = &mol();
154 molecule* molJ = &mol();
156 const directInteractionList& dil(il_.dil());
160 forAll(cellOccupancy_[
d],cellIMols)
162 molI = cellOccupancy_[
d][cellIMols];
164 forAll(dil[d], interactingCells)
166 List< molecule* > cellJ =
167 cellOccupancy_[dil[
d][interactingCells]];
171 molJ = cellJ[cellJMols];
173 evaluatePair(molI, molJ);
177 forAll(cellOccupancy_[d],cellIOtherMols)
179 molJ = cellOccupancy_[
d][cellIOtherMols];
183 evaluatePair(molI, molJ);
193 molecule* molK = &mol();
195 referredCellList& ril(il_.ril());
199 const List<label>& realCells = ril[r].realCellsForInteraction();
203 referredMolecule* molL(&(ril[r][refMols]));
207 List<molecule*> cellK = cellOccupancy_[realCells[rC]];
211 molK = cellK[cellKMols];
213 evaluatePair(molK, molL);
222 void Foam::moleculeCloud::calculateTetherForce()
224 const tetherPotentialList& tetherPot(pot_.tetherPotentials());
226 iterator mol(this->begin());
228 for (mol = this->begin(); mol != this->end(); ++mol)
230 if (mol().tethered())
232 vector rIT = mol().position() - mol().specialPosition();
234 label idI = mol().id();
236 scalar massI = constProps(idI).mass();
238 vector fIT = tetherPot.force(idI, rIT);
240 mol().a() += fIT/massI;
242 mol().potentialEnergy() += tetherPot.energy(idI, rIT);
251 void Foam::moleculeCloud::calculateExternalForce()
253 iterator mol(this->begin());
255 for (mol = this->begin(); mol != this->end(); ++mol)
257 mol().a() += pot_.gravity();
262 void Foam::moleculeCloud::removeHighEnergyOverlaps()
264 Info<<
nl <<
"Removing high energy overlaps, limit = "
265 << pot_.potentialEnergyLimit()
266 <<
nl <<
"Removal order:";
268 forAll(pot_.removalOrder(), rO)
270 Info<<
' ' << pot_.idList()[pot_.removalOrder()[rO]];
275 label initialSize = this->size();
277 buildCellOccupancy();
280 iterator mol(this->begin());
285 molecule* molI = &mol();
287 molecule* molJ = &mol();
289 DynamicList<molecule*> molsToDelete;
291 const directInteractionList& dil(il_.dil());
295 forAll(cellOccupancy_[
d],cellIMols)
297 molI = cellOccupancy_[
d][cellIMols];
299 forAll(dil[d], interactingCells)
301 List< molecule* > cellJ =
302 cellOccupancy_[dil[
d][interactingCells]];
306 molJ = cellJ[cellJMols];
308 if (evaluatePotentialLimit(molI, molJ))
310 label idI = molI->id();
312 label idJ = molJ->id();
323 molsToDelete.append(molJ);
326 else if (
findIndex(molsToDelete, molI) == -1)
328 molsToDelete.append(molI);
335 forAll(cellOccupancy_[d],cellIOtherMols)
337 molJ = cellOccupancy_[
d][cellIOtherMols];
341 if (evaluatePotentialLimit(molI, molJ))
343 label idI = molI->id();
345 label idJ = molJ->id();
356 molsToDelete.append(molJ);
359 else if (
findIndex(molsToDelete, molI) == -1)
361 molsToDelete.append(molI);
368 forAll (molsToDelete, mTD)
370 deleteParticle(*(molsToDelete[mTD]));
374 buildCellOccupancy();
379 molecule* molK = &mol();
381 DynamicList<molecule*> molsToDelete;
383 referredCellList& ril(il_.ril());
387 referredCell& refCell = ril[r];
391 referredMolecule* molL = &(refCell[refMols]);
393 List <label> realCells = refCell.realCellsForInteraction();
397 label cellK = realCells[rC];
399 List<molecule*> cellKMols = cellOccupancy_[cellK];
403 molK = cellKMols[cKM];
405 if (evaluatePotentialLimit(molK, molL))
407 label idK = molK->id();
409 label idL = molL->id();
419 molsToDelete.append(molK);
432 && cellK <= refCell.sourceCell()
437 molsToDelete.append(molK);
448 molsToDelete.append(molK);
458 forAll (molsToDelete, mTD)
460 deleteParticle(*(molsToDelete[mTD]));
464 buildCellOccupancy();
466 label molsRemoved = initialSize - this->size();
470 reduce(molsRemoved, sumOp<label>());
473 Info<<
tab << molsRemoved <<
" molecules removed" <<
endl;
477 void Foam::moleculeCloud::initialiseMolecules
479 const IOdictionary& mdInitialiseDict
483 <<
"Initialising molecules in each zone specified in mdInitialiseDict."
488 if (!cellZones.size())
490 FatalErrorIn(
"void Foam::moleculeCloud::initialiseMolecules")
491 <<
"No cellZones found in the mesh."
497 const cellZone& zone(cellZones[z]);
501 if (!mdInitialiseDict.found(zone.name()))
503 Info<<
"No specification subDictionary for zone "
504 << zone.name() <<
" found, skipping." <<
endl;
508 const dictionary& zoneDict =
509 mdInitialiseDict.subDict(zone.name());
511 const scalar temperature
516 const vector bulkVelocity(zoneDict.lookup(
"bulkVelocity"));
518 List<word> latticeIds
520 zoneDict.lookup(
"latticeIds")
523 List<vector> latticePositions
525 zoneDict.lookup(
"latticePositions")
528 if(latticeIds.size() != latticePositions.size())
530 FatalErrorIn(
"Foam::moleculeCloud::initialiseMolecules")
531 <<
"latticeIds and latticePositions must be the same "
538 zoneDict.lookup(
"latticeCellShape")
541 scalar latticeCellScale = 0.0;
543 if (zoneDict.found(
"numberDensity"))
547 zoneDict.lookup(
"numberDensity")
550 if (numberDensity < VSMALL)
552 WarningIn(
"moleculeCloud::initialiseMolecules")
553 <<
"numberDensity too small, not filling zone "
554 << zone.name() <<
endl;
559 latticeCellScale =
pow
561 latticeIds.size()/(
det(latticeCellShape)*numberDensity),
565 else if (zoneDict.found(
"massDensity"))
567 scalar unitCellMass = 0.0;
571 label
id =
findIndex(pot_.idList(), latticeIds[i]);
573 const molecule::constantProperties& cP(constProps(
id));
575 unitCellMass += cP.mass();
580 zoneDict.lookup(
"massDensity")
583 if (massDensity < VSMALL)
585 WarningIn(
"moleculeCloud::initialiseMolecules")
586 <<
"massDensity too small, not filling zone "
587 << zone.name() <<
endl;
593 latticeCellScale =
pow
595 unitCellMass/(
det(latticeCellShape)*massDensity),
601 FatalErrorIn(
"Foam::moleculeCloud::initialiseMolecules")
602 <<
"massDensity or numberDensity not specified " << nl
606 latticeCellShape *= latticeCellScale;
608 vector anchor(zoneDict.lookup(
"anchor"));
610 bool tethered =
false;
612 if (zoneDict.found(
"tetherSiteIds"))
616 List<word>(zoneDict.lookup(
"tetherSiteIds")).size()
620 const vector orientationAngles
622 zoneDict.lookup(
"orientationAngles")
663 const point cellCentre = mesh_.cellCentres()[zone[cell]];
665 if (cellCentre.x() > zoneMax.x())
667 zoneMax.
x() = cellCentre.x();
669 if (cellCentre.x() < zoneMin.x())
671 zoneMin.x() = cellCentre.x();
673 if (cellCentre.y() > zoneMax.y())
675 zoneMax.y() = cellCentre.y();
677 if (cellCentre.y() < zoneMin.y())
679 zoneMin.y() = cellCentre.y();
681 if (cellCentre.z() > zoneMax.z())
683 zoneMax.z() = cellCentre.z();
685 if (cellCentre.z() < zoneMin.z())
687 zoneMin.z() = cellCentre.z();
691 point zoneMid = 0.5*(zoneMin + zoneMax);
693 point latticeMid =
inv(latticeCellShape) & (
R.T()
694 & (zoneMid - anchor));
698 label(latticeMid.x() + 0.5*
sign(latticeMid.x())),
699 label(latticeMid.y() + 0.5*
sign(latticeMid.y())),
700 label(latticeMid.z() + 0.5*
sign(latticeMid.z()))
703 anchor += (
R & (latticeCellShape & latticeAnchor));
713 label totalZoneMols = 0;
715 label molsPlacedThisIteration = 0;
719 molsPlacedThisIteration != 0
720 || totalZoneMols == 0
723 label sizeBeforeIteration = this->size();
725 bool partOfLayerInBounds =
false;
740 label
id =
findIndex(pot_.idList(), latticeIds[
p]);
742 const vector& latticePosition =
745 unitCellLatticePosition.x(),
746 unitCellLatticePosition.y(),
747 unitCellLatticePosition.z()
749 + latticePositions[
p];
751 point globalPosition =
752 anchor + (
R & (latticeCellShape & latticePosition));
754 partOfLayerInBounds =
755 mesh_.bounds().contains(globalPosition);
757 label cell = mesh_.findCell(globalPosition);
781 unitCellLatticePosition.z() = -n;
782 unitCellLatticePosition.z() <= n;
783 unitCellLatticePosition.z() += 2*n
788 unitCellLatticePosition.y() = -n;
789 unitCellLatticePosition.y() <= n;
790 unitCellLatticePosition.y()++
795 unitCellLatticePosition.x() = -n;
796 unitCellLatticePosition.x() <= n;
797 unitCellLatticePosition.x()++
808 const vector& latticePosition =
811 unitCellLatticePosition.x(),
812 unitCellLatticePosition.y(),
813 unitCellLatticePosition.z()
815 + latticePositions[
p];
817 point globalPosition = anchor
819 &(latticeCellShape & latticePosition));
821 partOfLayerInBounds =
822 mesh_.bounds().contains(globalPosition);
824 label cell = mesh_.findCell
848 unitCellLatticePosition.z() = -(n-1);
849 unitCellLatticePosition.z() <= (n-1);
850 unitCellLatticePosition.z()++
853 for (label iR = 0; iR <= 2*n -1; iR++)
855 unitCellLatticePosition.x() = n;
857 unitCellLatticePosition.y() = -n + (iR + 1);
859 for (label iK = 0; iK < 4; iK++)
869 const vector& latticePosition =
872 unitCellLatticePosition.x(),
873 unitCellLatticePosition.y(),
874 unitCellLatticePosition.z()
876 + latticePositions[
p];
878 point globalPosition = anchor
880 &(latticeCellShape & latticePosition));
882 partOfLayerInBounds =
883 mesh_.bounds().contains(globalPosition);
885 label cell = mesh_.findCell
904 unitCellLatticePosition =
907 - unitCellLatticePosition.y(),
908 unitCellLatticePosition.x(),
909 unitCellLatticePosition.z()
923 && !partOfLayerInBounds
926 WarningIn(
"Foam::moleculeCloud::initialiseMolecules()")
927 <<
"A whole layer of unit cells was placed "
928 <<
"outside the bounds of the mesh, but no "
929 <<
"molecules have been placed in zone '"
931 <<
"'. This is likely to be because the zone "
934 <<
" in this case) and no lattice position "
935 <<
"fell inside them. "
936 <<
"Aborting filling this zone."
942 molsPlacedThisIteration =
943 this->size() - sizeBeforeIteration;
945 totalZoneMols += molsPlacedThisIteration;
955 void Foam::moleculeCloud::createMolecule
957 const point& position,
962 const vector& bulkVelocity
965 const Cloud<molecule>& cloud = *
this;
969 cell = mesh_.findCell(position);
975 <<
"Position specified does not correspond to a mesh cell." << nl
985 specialPosition = position;
990 const molecule::constantProperties& cP(constProps(
id));
992 vector v = equipartitionLinearVelocity(temperature, cP.mass());
1000 if (!cP.pointMolecule())
1002 pi = equipartitionAngularMomentum(temperature, cP);
1045 Foam::label Foam::moleculeCloud::nSites()
const
1049 const_iterator mol(this->begin());
1051 for (mol = this->begin(); mol != this->end(); ++mol)
1053 n += constProps(mol().
id()).nSites();
1062 Foam::moleculeCloud::moleculeCloud
1072 cellOccupancy_(mesh_.nCells()),
1073 il_(mesh_, pot_.pairPotentials().rCutMaxSqr(),
false),
1084 setSiteSizesAndPositions();
1086 removeHighEnergyOverlaps();
1092 Foam::moleculeCloud::moleculeCloud
1116 initialiseMolecules(mdInitialiseDict);
1142 buildCellOccupancy();
1147 for (mol = this->begin(); mol != this->end(); ++mol)
1151 mol().potentialEnergy() = 0.0;
1156 calculatePairForce();
1158 calculateTetherForce();
1160 calculateExternalForce();
1166 const scalar targetTemperature,
1167 const scalar measuredTemperature
1170 scalar temperatureCorrectionFactor =
1171 sqrt(targetTemperature/measuredTemperature);
1173 Info<<
"----------------------------------------" << nl
1174 <<
"Temperature equilibration" << nl
1175 <<
"Target temperature = "
1176 << targetTemperature << nl
1177 <<
"Measured temperature = "
1178 << measuredTemperature << nl
1179 <<
"Temperature correction factor = "
1180 << temperatureCorrectionFactor << nl
1181 <<
"----------------------------------------"
1186 for (mol = this->begin(); mol != this->end(); ++mol)
1188 mol().v() *= temperatureCorrectionFactor;
1190 mol().pi() *= temperatureCorrectionFactor;
1199 str << nSites() << nl <<
"moleculeCloud site positions in angstroms" <<
nl;
1203 for (mol = this->begin(); mol != this->end(); ++mol)
1207 forAll(mol().sitePositions(), i)
1209 const point& sP = mol().sitePositions()[i];
1211 str << pot_.siteIdList()[cP.
siteIds()[i]]
1212 <<
' ' << sP.
x()*1e10
1213 <<
' ' << sP.
y()*1e10
1214 <<
' ' << sP.
z()*1e10