moleculeCloud.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2025 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "moleculeCloud.H"
27 #include "fvMesh.H"
28 #include "meshSearch.H"
29 #include "mathematicalConstants.H"
30 
31 using namespace Foam::constant::mathematical;
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(moleculeCloud, 0);
38 }
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::moleculeCloud::buildConstProps()
43 {
44  Info<< nl << "Reading moleculeProperties dictionary." << endl;
45 
46  const List<word>& idList(pot_.idList());
47 
48  constPropList_.setSize(idList.size());
49 
50  const List<word>& siteIdList(pot_.siteIdList());
51 
52  IOdictionary moleculePropertiesDict
53  (
54  IOobject
55  (
56  "moleculeProperties",
57  mesh_.time().constant(),
58  mesh_,
59  IOobject::MUST_READ_IF_MODIFIED,
60  IOobject::NO_WRITE,
61  false
62  )
63  );
64 
65  forAll(idList, i)
66  {
67  const word& id = idList[i];
68  const dictionary& molDict = moleculePropertiesDict.subDict(id);
69 
70  List<word> siteIdNames = molDict.lookup("siteIds");
71 
72  List<label> siteIds(siteIdNames.size());
73 
74  forAll(siteIdNames, sI)
75  {
76  const word& siteId = siteIdNames[sI];
77 
78  siteIds[sI] = findIndex(siteIdList, siteId);
79 
80  if (siteIds[sI] == -1)
81  {
82  FatalErrorInFunction
83  << siteId << " site not found."
84  << nl << abort(FatalError);
85  }
86  }
87 
88  molecule::constantProperties& constProp = constPropList_[i];
89 
90  constProp = molecule::constantProperties(molDict);
91 
92  constProp.siteIds() = siteIds;
93  }
94 }
95 
96 
97 void Foam::moleculeCloud::setSiteSizesAndPositions()
98 {
99  forAllIter(moleculeCloud, *this, mol)
100  {
101  const molecule::constantProperties& cP = constProps(mol().id());
102 
103  mol().setSiteSizes(cP.nSites());
104 
105  mol().setSitePositions(mesh(), cP);
106  }
107 }
108 
109 
110 void Foam::moleculeCloud::buildCellOccupancy()
111 {
112  forAll(cellOccupancy_, cO)
113  {
114  cellOccupancy_[cO].clear();
115  }
116 
117  forAllIter(moleculeCloud, *this, mol)
118  {
119  cellOccupancy_[mol().cell()].append(&mol());
120  }
121 
122  forAll(cellOccupancy_, cO)
123  {
124  cellOccupancy_[cO].shrink();
125  }
126 }
127 
128 
129 void Foam::moleculeCloud::calculatePairForce()
130 {
131  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
132 
133  // Start sending referred data
134  label startOfRequests = Pstream::nRequests();
135  il_.sendReferredData(cellOccupancy(), pBufs);
136 
137  molecule* molI = nullptr;
138  molecule* molJ = nullptr;
139 
140  {
141  // Real-Real interactions
142 
143  const labelListList& dil = il_.dil();
144 
145  forAll(dil, d)
146  {
147  forAll(cellOccupancy_[d],cellIMols)
148  {
149  molI = cellOccupancy_[d][cellIMols];
150 
151  forAll(dil[d], interactingCells)
152  {
153  List<molecule*> cellJ =
154  cellOccupancy_[dil[d][interactingCells]];
155 
156  forAll(cellJ, cellJMols)
157  {
158  molJ = cellJ[cellJMols];
159 
160  evaluatePair(*molI, *molJ);
161  }
162  }
163 
164  forAll(cellOccupancy_[d], cellIOtherMols)
165  {
166  molJ = cellOccupancy_[d][cellIOtherMols];
167 
168  if (molJ > molI)
169  {
170  evaluatePair(*molI, *molJ);
171  }
172  }
173  }
174  }
175  }
176 
177  // Receive referred data
178  il_.receiveReferredData(pBufs, startOfRequests);
179 
180  {
181  // Real-Referred interactions
182 
183  const labelListList& ril = il_.ril();
184 
185  List<IDLList<molecule>>& referredMols = il_.referredParticles();
186 
187  forAll(ril, r)
188  {
189  const List<label>& realCells = ril[r];
190 
191  IDLList<molecule>& refMols = referredMols[r];
192 
193  forAllIter
194  (
195  IDLList<molecule>,
196  refMols,
197  refMol
198  )
199  {
200  forAll(realCells, rC)
201  {
202  List<molecule*> celli = cellOccupancy_[realCells[rC]];
203 
204  forAll(celli, cellIMols)
205  {
206  molI = celli[cellIMols];
207 
208  evaluatePair(*molI, refMol());
209  }
210  }
211  }
212  }
213  }
214 }
215 
216 
217 void Foam::moleculeCloud::calculateTetherForce()
218 {
219  const tetherPotentialList& tetherPot(pot_.tetherPotentials());
220 
221  forAllIter(moleculeCloud, *this, mol)
222  {
223  if (mol().tethered())
224  {
225  vector rIT = mol().position(mesh()) - mol().specialPosition();
226 
227  label idI = mol().id();
228 
229  scalar massI = constProps(idI).mass();
230 
231  vector fIT = tetherPot.force(idI, rIT);
232 
233  mol().a() += fIT/massI;
234 
235  mol().potentialEnergy() += tetherPot.energy(idI, rIT);
236 
237  // What to do here?
238  // mol().rf() += rIT*fIT;
239  }
240  }
241 }
242 
243 
244 void Foam::moleculeCloud::calculateExternalForce()
245 {
246  forAllIter(moleculeCloud, *this, mol)
247  {
248  mol().a() += pot_.gravity();
249  }
250 }
251 
252 
253 void Foam::moleculeCloud::removeHighEnergyOverlaps()
254 {
255  Info<< nl << "Removing high energy overlaps, limit = "
256  << pot_.potentialEnergyLimit()
257  << nl << "Removal order:";
258 
259  forAll(pot_.removalOrder(), rO)
260  {
261  Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
262  }
263 
264  Info<< nl ;
265 
266  label initialSize = this->size();
267 
268  buildCellOccupancy();
269 
270  // Real-Real interaction
271 
272  molecule* molI = nullptr;
273  molecule* molJ = nullptr;
274 
275  {
276  DynamicList<molecule*> molsToDelete;
277 
278  const labelListList& dil(il_.dil());
279 
280  forAll(dil, d)
281  {
282  forAll(cellOccupancy_[d],cellIMols)
283  {
284  molI = cellOccupancy_[d][cellIMols];
285 
286  forAll(dil[d], interactingCells)
287  {
288  List<molecule*> cellJ =
289  cellOccupancy_[dil[d][interactingCells]];
290 
291  forAll(cellJ, cellJMols)
292  {
293  molJ = cellJ[cellJMols];
294 
295  if (evaluatePotentialLimit(*molI, *molJ))
296  {
297  label idI = molI->id();
298 
299  label idJ = molJ->id();
300 
301  if
302  (
303  idI == idJ
304  || findIndex(pot_.removalOrder(), idJ)
305  < findIndex(pot_.removalOrder(), idI)
306  )
307  {
308  if (findIndex(molsToDelete, molJ) == -1)
309  {
310  molsToDelete.append(molJ);
311  }
312  }
313  else if (findIndex(molsToDelete, molI) == -1)
314  {
315  molsToDelete.append(molI);
316  }
317  }
318  }
319  }
320  }
321 
322  forAll(cellOccupancy_[d], cellIOtherMols)
323  {
324  molJ = cellOccupancy_[d][cellIOtherMols];
325 
326  if (molJ > molI)
327  {
328  if (evaluatePotentialLimit(*molI, *molJ))
329  {
330  label idI = molI->id();
331 
332  label idJ = molJ->id();
333 
334  if
335  (
336  idI == idJ
337  || findIndex(pot_.removalOrder(), idJ)
338  < findIndex(pot_.removalOrder(), idI)
339  )
340  {
341  if (findIndex(molsToDelete, molJ) == -1)
342  {
343  molsToDelete.append(molJ);
344  }
345  }
346  else if (findIndex(molsToDelete, molI) == -1)
347  {
348  molsToDelete.append(molI);
349  }
350  }
351  }
352  }
353  }
354 
355  forAll(molsToDelete, mTD)
356  {
357  deleteParticle(*(molsToDelete[mTD]));
358  }
359  }
360 
361  buildCellOccupancy();
362 
363  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
364 
365  // Start sending referred data
366  label startOfRequests = Pstream::nRequests();
367 
368  il_.sendReferredData(cellOccupancy(), pBufs);
369 
370  // Receive referred data
371  il_.receiveReferredData(pBufs, startOfRequests);
372 
373  // Real-Referred interaction
374 
375  {
376  DynamicList<molecule*> molsToDelete;
377 
378  const labelListList& ril(il_.ril());
379 
380  List<IDLList<molecule>>& referredMols = il_.referredParticles();
381 
382  forAll(ril, r)
383  {
384  IDLList<molecule>& refMols = referredMols[r];
385 
386  forAllIter
387  (
388  IDLList<molecule>,
389  refMols,
390  refMol
391  )
392  {
393  molJ = &refMol();
394 
395  const List<label>& realCells = ril[r];
396 
397  forAll(realCells, rC)
398  {
399  label celli = realCells[rC];
400 
401  List<molecule*> cellIMols = cellOccupancy_[celli];
402 
403  forAll(cellIMols, cIM)
404  {
405  molI = cellIMols[cIM];
406 
407  if (evaluatePotentialLimit(*molI, *molJ))
408  {
409  label idI = molI->id();
410 
411  label idJ = molJ->id();
412 
413  if
414  (
415  findIndex(pot_.removalOrder(), idI)
416  < findIndex(pot_.removalOrder(), idJ)
417  )
418  {
419  if (findIndex(molsToDelete, molI) == -1)
420  {
421  molsToDelete.append(molI);
422  }
423  }
424  else if
425  (
426  findIndex(pot_.removalOrder(), idI)
427  == findIndex(pot_.removalOrder(), idJ)
428  )
429  {
430  // Remove one of the molecules
431  // arbitrarily, assuring that a
432  // consistent decision is made for
433  // both real-referred pairs.
434 
435  if (molI->origId() > molJ->origId())
436  {
437  if (findIndex(molsToDelete, molI) == -1)
438  {
439  molsToDelete.append(molI);
440  }
441  }
442  }
443  }
444  }
445  }
446  }
447  }
448 
449  forAll(molsToDelete, mTD)
450  {
451  deleteParticle(*(molsToDelete[mTD]));
452  }
453  }
454 
455  buildCellOccupancy();
456 
457  // Start sending referred data
458  startOfRequests = Pstream::nRequests();
459 
460  il_.sendReferredData(cellOccupancy(), pBufs);
461 
462  // Receive referred data
463  il_.receiveReferredData(pBufs, startOfRequests);
464 
465  label molsRemoved = initialSize - this->size();
466 
467  if (Pstream::parRun())
468  {
469  reduce(molsRemoved, sumOp<label>());
470  }
471 
472  Info<< tab << molsRemoved << " molecules removed" << endl;
473 }
474 
475 
476 void Foam::moleculeCloud::initialiseMolecules
477 (
478  const IOdictionary& mdInitialiseDict
479 )
480 {
481  Info<< nl
482  << "Initialising molecules in each zone specified in mdInitialiseDict."
483  << endl;
484 
485  const cellZoneList& cellZones = mesh_.cellZones();
486 
487  const meshSearch& searchEngine = meshSearch::New(mesh());
488 
489  if (!cellZones.size())
490  {
491  FatalErrorInFunction
492  << "No cellZones found in the mesh."
493  << abort(FatalError);
494  }
495 
496  label nLocateBoundaryHits = 0;
497 
498  forAll(cellZones, z)
499  {
500  const cellZone& zone(cellZones[z]);
501 
502  if (zone.size())
503  {
504  if (!mdInitialiseDict.found(zone.name()))
505  {
506  Info<< "No specification subDictionary for zone "
507  << zone.name() << " found, skipping." << endl;
508  }
509  else
510  {
511  const dictionary& zoneDict =
512  mdInitialiseDict.subDict(zone.name());
513 
514  const scalar temperature
515  (
516  zoneDict.template lookup<scalar>("temperature")
517  );
518 
519  const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
520 
521  List<word> latticeIds
522  (
523  zoneDict.lookup("latticeIds")
524  );
525 
526  List<vector> latticePositions
527  (
528  zoneDict.lookup("latticePositions")
529  );
530 
531  if (latticeIds.size() != latticePositions.size())
532  {
533  FatalErrorInFunction
534  << "latticeIds and latticePositions must be the same "
535  << " size." << nl
536  << abort(FatalError);
537  }
538 
539  diagTensor latticeCellShape
540  (
541  zoneDict.lookup("latticeCellShape")
542  );
543 
544  scalar latticeCellScale = 0.0;
545 
546  if (zoneDict.found("numberDensity"))
547  {
548  scalar numberDensity =
549  zoneDict.lookup<scalar>("numberDensity");
550 
551  if (numberDensity < vSmall)
552  {
553  WarningInFunction
554  << "numberDensity too small, not filling zone "
555  << zone.name() << endl;
556 
557  continue;
558  }
559 
560  latticeCellScale = pow
561  (
562  latticeIds.size()/(det(latticeCellShape)*numberDensity),
563  (1.0/3.0)
564  );
565  }
566  else if (zoneDict.found("massDensity"))
567  {
568  scalar unitCellMass = 0.0;
569 
570  forAll(latticeIds, i)
571  {
572  label id = findIndex(pot_.idList(), latticeIds[i]);
573 
574  const molecule::constantProperties& cP(constProps(id));
575 
576  unitCellMass += cP.mass();
577  }
578 
579  scalar massDensity = zoneDict.lookup<scalar>("massDensity");
580 
581  if (massDensity < vSmall)
582  {
583  WarningInFunction
584  << "massDensity too small, not filling zone "
585  << zone.name() << endl;
586 
587  continue;
588  }
589 
590 
591  latticeCellScale = pow
592  (
593  unitCellMass/(det(latticeCellShape)*massDensity),
594  (1.0/3.0)
595  );
596  }
597  else
598  {
599  FatalErrorInFunction
600  << "massDensity or numberDensity not specified " << nl
601  << abort(FatalError);
602  }
603 
604  latticeCellShape *= latticeCellScale;
605 
606  vector anchor(zoneDict.lookup("anchor"));
607 
608  bool tethered = false;
609 
610  if (zoneDict.found("tetherSiteIds"))
611  {
612  tethered = bool
613  (
614  List<word>(zoneDict.lookup("tetherSiteIds")).size()
615  );
616  }
617 
618  const vector orientationAngles
619  (
620  zoneDict.lookup("orientationAngles")
621  );
622 
623  scalar phi(orientationAngles.x()*pi/180.0);
624 
625  scalar theta(orientationAngles.y()*pi/180.0);
626 
627  scalar psi(orientationAngles.z()*pi/180.0);
628 
629  const tensor R
630  (
631  cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
632  cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
633  sin(psi)*sin(theta),
634  - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
635  - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
636  cos(psi)*sin(theta),
637  sin(theta)*sin(phi),
638  - sin(theta)*cos(phi),
639  cos(theta)
640  );
641 
642  // Find the optimal anchor position. Finding the approximate
643  // mid-point of the zone of cells and snapping to the nearest
644  // lattice location.
645 
646  vector zoneMin = vGreat*vector::one;
647 
648  vector zoneMax = -vGreat*vector::one;
649 
650  forAll(zone, cell)
651  {
652  const point cellCentre = mesh_.cellCentres()[zone[cell]];
653 
654  if (cellCentre.x() > zoneMax.x())
655  {
656  zoneMax.x() = cellCentre.x();
657  }
658  if (cellCentre.x() < zoneMin.x())
659  {
660  zoneMin.x() = cellCentre.x();
661  }
662  if (cellCentre.y() > zoneMax.y())
663  {
664  zoneMax.y() = cellCentre.y();
665  }
666  if (cellCentre.y() < zoneMin.y())
667  {
668  zoneMin.y() = cellCentre.y();
669  }
670  if (cellCentre.z() > zoneMax.z())
671  {
672  zoneMax.z() = cellCentre.z();
673  }
674  if (cellCentre.z() < zoneMin.z())
675  {
676  zoneMin.z() = cellCentre.z();
677  }
678  }
679 
680  point zoneMid = 0.5*(zoneMin + zoneMax);
681 
682  point latticeMid = inv(latticeCellShape) & (R.T()
683  & (zoneMid - anchor));
684 
685  point latticeAnchor
686  (
687  label(latticeMid.x() + 0.5*sign(latticeMid.x())),
688  label(latticeMid.y() + 0.5*sign(latticeMid.y())),
689  label(latticeMid.z() + 0.5*sign(latticeMid.z()))
690  );
691 
692  anchor += (R & (latticeCellShape & latticeAnchor));
693 
694  // Continue trying to place molecules as long as at
695  // least one molecule is placed in each iteration.
696  // The "|| totalZoneMols == 0" condition means that the
697  // algorithm will continue if the origin is outside the
698  // zone.
699 
700  label n = 0;
701 
702  label totalZoneMols = 0;
703 
704  label molsPlacedThisIteration = 0;
705 
706  while
707  (
708  molsPlacedThisIteration != 0
709  || totalZoneMols == 0
710  )
711  {
712  label sizeBeforeIteration = this->size();
713 
714  bool partOfLayerInBounds = false;
715 
716  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
717  // start of placement of molecules
718  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
719 
720  if (n == 0)
721  {
722  // Special treatment is required for the first position,
723  // i.e. iteration zero.
724 
725  labelVector unitCellLatticePosition(0,0,0);
726 
727  forAll(latticePositions, p)
728  {
729  label id = findIndex(pot_.idList(), latticeIds[p]);
730 
731  const vector& latticePosition =
732  vector
733  (
734  unitCellLatticePosition.x(),
735  unitCellLatticePosition.y(),
736  unitCellLatticePosition.z()
737  )
738  + latticePositions[p];
739 
740  point globalPosition =
741  anchor
742  + (R & (latticeCellShape & latticePosition));
743 
744  partOfLayerInBounds = mesh_.bounds().contains
745  (
746  globalPosition
747  );
748 
749  const label cell =
750  searchEngine.findCell(globalPosition);
751 
752  if (findIndex(zone, cell) != -1)
753  {
754  createMolecule
755  (
756  searchEngine,
757  globalPosition,
758  cell,
759  nLocateBoundaryHits,
760  id,
761  tethered,
762  temperature,
763  bulkVelocity
764  );
765  }
766  }
767  }
768  else
769  {
770  // Place top and bottom caps.
771 
772  labelVector unitCellLatticePosition(0,0,0);
773 
774  for
775  (
776  unitCellLatticePosition.z() = -n;
777  unitCellLatticePosition.z() <= n;
778  unitCellLatticePosition.z() += 2*n
779  )
780  {
781  for
782  (
783  unitCellLatticePosition.y() = -n;
784  unitCellLatticePosition.y() <= n;
785  unitCellLatticePosition.y()++
786  )
787  {
788  for
789  (
790  unitCellLatticePosition.x() = -n;
791  unitCellLatticePosition.x() <= n;
792  unitCellLatticePosition.x()++
793  )
794  {
795  forAll(latticePositions, p)
796  {
797  label id = findIndex
798  (
799  pot_.idList(),
800  latticeIds[p]
801  );
802 
803  const vector& latticePosition =
804  vector
805  (
806  unitCellLatticePosition.x(),
807  unitCellLatticePosition.y(),
808  unitCellLatticePosition.z()
809  )
810  + latticePositions[p];
811 
812  point globalPosition =
813  anchor
814  + (
815  R
816  & (
817  latticeCellShape
818  & latticePosition
819  )
820  );
821 
822  partOfLayerInBounds =
823  mesh_.bounds().contains
824  (
825  globalPosition
826  );
827 
828  const label cell =
829  searchEngine.findCell
830  (
831  globalPosition
832  );
833 
834  if (findIndex(zone, cell) != -1)
835  {
836  createMolecule
837  (
838  searchEngine,
839  globalPosition,
840  cell,
841  nLocateBoundaryHits,
842  id,
843  tethered,
844  temperature,
845  bulkVelocity
846  );
847  }
848  }
849  }
850  }
851  }
852 
853  for
854  (
855  unitCellLatticePosition.z() = -(n-1);
856  unitCellLatticePosition.z() <= (n-1);
857  unitCellLatticePosition.z()++
858  )
859  {
860  for (label iR = 0; iR <= 2*n -1; iR++)
861  {
862  unitCellLatticePosition.x() = n;
863 
864  unitCellLatticePosition.y() = -n + (iR + 1);
865 
866  for (label iK = 0; iK < 4; iK++)
867  {
868  forAll(latticePositions, p)
869  {
870  label id = findIndex
871  (
872  pot_.idList(),
873  latticeIds[p]
874  );
875 
876  const vector& latticePosition =
877  vector
878  (
879  unitCellLatticePosition.x(),
880  unitCellLatticePosition.y(),
881  unitCellLatticePosition.z()
882  )
883  + latticePositions[p];
884 
885  point globalPosition =
886  anchor
887  + (
888  R
889  & (
890  latticeCellShape
891  & latticePosition
892  )
893  );
894 
895  partOfLayerInBounds =
896  mesh_.bounds().contains
897  (
898  globalPosition
899  );
900 
901  const label cell =
902  searchEngine.findCell
903  (
904  globalPosition
905  );
906 
907  if (findIndex(zone, cell) != -1)
908  {
909  createMolecule
910  (
911  searchEngine,
912  globalPosition,
913  cell,
914  nLocateBoundaryHits,
915  id,
916  tethered,
917  temperature,
918  bulkVelocity
919  );
920  }
921  }
922 
923  unitCellLatticePosition =
924  labelVector
925  (
926  - unitCellLatticePosition.y(),
927  unitCellLatticePosition.x(),
928  unitCellLatticePosition.z()
929  );
930  }
931  }
932  }
933  }
934 
935  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
936  // end of placement of molecules
937  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
938 
939  if
940  (
941  totalZoneMols == 0
942  && !partOfLayerInBounds
943  )
944  {
945  WarningInFunction
946  << "A whole layer of unit cells was placed "
947  << "outside the bounds of the mesh, but no "
948  << "molecules have been placed in zone '"
949  << zone.name()
950  << "'. This is likely to be because the zone "
951  << "has few cells ("
952  << zone.size()
953  << " in this case) and no lattice position "
954  << "fell inside them. "
955  << "Aborting filling this zone."
956  << endl;
957 
958  break;
959  }
960 
961  molsPlacedThisIteration =
962  this->size() - sizeBeforeIteration;
963 
964  totalZoneMols += molsPlacedThisIteration;
965 
966  n++;
967  }
968  }
969  }
970  }
971 
972  reduce(nLocateBoundaryHits, sumOp<label>());
973  if (nLocateBoundaryHits != 0)
974  {
975  WarningInFunction
976  << "Initialisation of cloud " << this->name()
977  << " did not accurately locate " << nLocateBoundaryHits
978  << " particles" << endl;
979  }
980 }
981 
982 
983 void Foam::moleculeCloud::createMolecule
984 (
985  const meshSearch& searchEngine,
986  const point& position,
987  label cell,
988  label& nLocateBoundaryHits,
989  label id,
990  bool tethered,
991  scalar temperature,
992  const vector& bulkVelocity
993 )
994 {
995  point specialPosition(Zero);
996 
997  label special = 0;
998 
999  if (tethered)
1000  {
1001  specialPosition = position;
1002 
1003  special = molecule::SPECIAL_TETHERED;
1004  }
1005 
1006  const molecule::constantProperties& cP(constProps(id));
1007 
1008  vector v = equipartitionLinearVelocity(temperature, cP.mass());
1009 
1010  v += bulkVelocity;
1011 
1012  vector pi = Zero;
1013 
1014  tensor Q = I;
1015 
1016  if (!cP.pointMolecule())
1017  {
1018  pi = equipartitionAngularMomentum(temperature, cP);
1019 
1020  scalar phi(rndGen_.scalar01()*twoPi);
1021 
1022  scalar theta(rndGen_.scalar01()*twoPi);
1023 
1024  scalar psi(rndGen_.scalar01()*twoPi);
1025 
1026  Q = tensor
1027  (
1028  cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
1029  cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
1030  sin(psi)*sin(theta),
1031  - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
1032  - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
1033  cos(psi)*sin(theta),
1034  sin(theta)*sin(phi),
1035  - sin(theta)*cos(phi),
1036  cos(theta)
1037  );
1038  }
1039 
1040  addParticle
1041  (
1042  new molecule
1043  (
1044  searchEngine,
1045  position,
1046  cell,
1047  nLocateBoundaryHits,
1048  Q,
1049  v,
1050  Zero,
1051  pi,
1052  Zero,
1053  specialPosition,
1054  constProps(id),
1055  special,
1056  id
1057  )
1058  );
1059 }
1060 
1061 
1062 Foam::label Foam::moleculeCloud::nSites() const
1063 {
1064  label n = 0;
1065 
1066  forAllConstIter(moleculeCloud, *this, mol)
1067  {
1068  n += constProps(mol().id()).nSites();
1069  }
1070 
1071  return n;
1072 }
1073 
1074 
1075 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1076 
1077 Foam::moleculeCloud::moleculeCloud
1078 (
1079  const polyMesh& mesh,
1080  const potential& pot,
1081  bool readFields
1082 )
1083 :
1084  lagrangian::Cloud<molecule>(mesh, "moleculeCloud", false),
1085  mesh_(mesh),
1086  pot_(pot),
1087  cellOccupancy_(mesh_.nCells()),
1088  il_(mesh_, pot_.pairPotentials().rCutMax(), false),
1089  constPropList_(),
1090  rndGen_(clock::getTime()),
1091  stdNormal_(rndGen_.generator())
1092 {
1093  if (readFields)
1094  {
1095  molecule::readFields(*this);
1096  }
1097 
1098  buildConstProps();
1099 
1100  setSiteSizesAndPositions();
1101 
1102  removeHighEnergyOverlaps();
1103 
1104  calculateForce();
1105 }
1106 
1107 
1108 Foam::moleculeCloud::moleculeCloud
1109 (
1110  const polyMesh& mesh,
1111  const potential& pot,
1112  const IOdictionary& mdInitialiseDict,
1113  bool readFields
1114 )
1115 :
1116  lagrangian::Cloud<molecule>(mesh, "moleculeCloud", false),
1117  mesh_(mesh),
1118  pot_(pot),
1119  il_(mesh_, 0.0, false),
1120  constPropList_(),
1121  rndGen_(clock::getTime()),
1122  stdNormal_(rndGen_.generator())
1123 {
1124  if (readFields)
1125  {
1126  molecule::readFields(*this);
1127  }
1128 
1129  clear();
1130 
1131  buildConstProps();
1132 
1133  initialiseMolecules(mdInitialiseDict);
1134 }
1135 
1136 
1137 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1138 
1139 void Foam::moleculeCloud::evolve()
1140 {
1141  molecule::trackingData td(*this);
1142 
1143  td.part() = molecule::trackingData::tpVelocityHalfStep0;
1144  lagrangian::Cloud<molecule>::move(*this, td);
1145 
1146  td.part() = molecule::trackingData::tpLinearTrack;
1147  lagrangian::Cloud<molecule>::move(*this, td);
1148 
1149  td.part() = molecule::trackingData::tpRotationalTrack;
1150  lagrangian::Cloud<molecule>::move(*this, td);
1151 
1152  calculateForce();
1153 
1154  td.part() = molecule::trackingData::tpVelocityHalfStep1;
1155  lagrangian::Cloud<molecule>::move(*this, td);
1156 }
1157 
1158 
1159 void Foam::moleculeCloud::calculateForce()
1160 {
1161  buildCellOccupancy();
1162 
1163  // Set accumulated quantities to zero
1164  forAllIter(moleculeCloud, *this, mol)
1165  {
1166  mol().siteForces() = Zero;
1167 
1168  mol().potentialEnergy() = 0.0;
1169 
1170  mol().rf() = Zero;
1171  }
1172 
1173  calculatePairForce();
1174 
1175  calculateTetherForce();
1176 
1177  calculateExternalForce();
1178 }
1179 
1180 
1181 void Foam::moleculeCloud::applyConstraintsAndThermostats
1182 (
1183  const scalar targetTemperature,
1184  const scalar measuredTemperature
1185 )
1186 {
1187  scalar temperatureCorrectionFactor =
1188  sqrt(targetTemperature/max(vSmall, measuredTemperature));
1189 
1190  Info<< "----------------------------------------" << nl
1191  << "Temperature equilibration" << nl
1192  << "Target temperature = "
1193  << targetTemperature << nl
1194  << "Measured temperature = "
1195  << measuredTemperature << nl
1196  << "Temperature correction factor = "
1197  << temperatureCorrectionFactor << nl
1198  << "----------------------------------------"
1199  << endl;
1200 
1201  forAllIter(moleculeCloud, *this, mol)
1202  {
1203  mol().v() *= temperatureCorrectionFactor;
1204 
1205  mol().pi() *= temperatureCorrectionFactor;
1206  }
1207 }
1208 
1209 
1210 void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
1211 {
1212  OFstream os(fName);
1213 
1214  os << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
1215 
1216  forAllConstIter(moleculeCloud, *this, mol)
1217  {
1218  const molecule::constantProperties& cP = constProps(mol().id());
1219 
1220  forAll(mol().sitePositions(), i)
1221  {
1222  const point& sP = mol().sitePositions()[i];
1223 
1224  os << pot_.siteIdList()[cP.siteIds()[i]]
1225  << ' ' << sP.x()*1e10
1226  << ' ' << sP.y()*1e10
1227  << ' ' << sP.z()*1e10
1228  << nl;
1229  }
1230  }
1231 }
1232 
1233 
1234 // ************************************************************************* //
n
label n
Definition: TABSMDCalcMethod2.H:31
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:449
forAllIter
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:474
forAllConstIter
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:492
cellOccupancy
const List< DynamicList< molecule * > > & cellOccupancy
Definition: calculateMDFields.H:1
Cloud
Base cloud calls templated on particle type.
Foam::IOdictionary
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Foam::IOobject::MUST_READ_IF_MODIFIED
@ MUST_READ_IF_MODIFIED
Definition: IOobject.H:119
Foam::OFstream
Output to file stream.
Definition: OFstream.H:87
Foam::UPstream::nRequests
static label nRequests()
Get number of outstanding requests.
Definition: UPstream.C:137
Foam::UPstream::parRun
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
Foam::VectorSpace::one
static const Form one
Definition: VectorSpace.H:119
Foam::Vector::z
const Cmpt & z() const
Definition: VectorI.H:87
Foam::Vector::y
const Cmpt & y() const
Definition: VectorI.H:81
Foam::Vector::x
const Cmpt & x() const
Definition: VectorI.H:75
Foam::clock
Read access to the system clock with formatting.
Definition: clock.H:51
Foam::fileName
A class for handling file names.
Definition: fileName.H:82
Foam::lagrangian::Cloud::move
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td)
Move the particles.
Definition: Cloud.C:292
Foam::meshSearch::New
static const meshSearch & New(const polyMesh &mesh, const pointInCellShapes=pointInCellShapes::tets)
Lookup or construct from mesh and cell decomposition option.
Definition: meshSearch.C:61
Foam::moleculeCloud::writeXYZ
void writeXYZ(const fileName &fName) const
Write molecule sites in XYZ format.
Definition: moleculeCloud.C:1210
Foam::moleculeCloud::evolve
void evolve()
Evolve the molecules (move, calculate forces, control state etc)
Definition: moleculeCloud.C:1139
Foam::moleculeCloud::applyConstraintsAndThermostats
void applyConstraintsAndThermostats(const scalar targetTemperature, const scalar measuredTemperature)
Definition: moleculeCloud.C:1182
Foam::moleculeCloud::moleculeCloud
moleculeCloud(const polyMesh &mesh, const potential &pot, bool readFields=true)
Construct given mesh and potential references.
Definition: moleculeCloud.C:1078
Foam::molecule::constantProperties
Class to hold molecule constant properties.
Definition: molecule.H:91
Foam::molecule::constantProperties::siteIds
const List< label > & siteIds() const
Definition: moleculeI.H:353
Foam::molecule::trackingData
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:169
Foam::molecule::trackingData::part
trackPart part() const
Return the part of the tracking operation taking place.
Definition: molecule.H:203
Foam::molecule
Foam::molecule.
Definition: molecule.H:68
Foam::molecule::SPECIAL_TETHERED
@ SPECIAL_TETHERED
Definition: molecule.H:83
Foam::molecule::readFields
static void readFields(lagrangian::Cloud< molecule > &mC)
Definition: moleculeIO.C:89
List
Motion of the mesh specified as a list of pointMeshMovers.
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:78
mesh
Foam::fvMesh mesh(Foam::IOobject(regionName, runTime.name(), runTime, Foam::IOobject::MUST_READ), false)
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:250
Foam::constant::mathematical
mathematical constants.
Foam::constant::mathematical::twoPi
const scalar twoPi(2 *pi)
Foam::dimensions::temperature
const dimensionSet temperature
Foam::tracking::position
point position(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label facei, const label faceTrii, const scalar stepFraction)
Return the position given the coordinates and tet topology.
Definition: trackingI.H:224
Foam
Namespace for OpenFOAM.
Definition: atmosphericBoundaryLayer.C:32
Foam::Zero
static const zero Zero
Definition: zero.H:97
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:179
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:288
Foam::readFields
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const HashSet< word > &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the specified type.
Definition: ReadFields.C:244
Foam::labelVector
Vector< label > labelVector
Vector of labels.
Definition: labelVector.H:49
Foam::tab
static const char tab
Definition: Ostream.H:296
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:277
Foam::Info
messageStream Info
Foam::I
static const Identity< scalar > I
Definition: Identity.H:93
Foam::point
vector point
Point is a vector.
Definition: point.H:41
Foam::diagTensor
DiagTensor< scalar > diagTensor
A scalar version of the templated DiagTensor.
Definition: diagTensor.H:48
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:49
Foam::inv
void inv(pointPatchField< tensor > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:297
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Foam::R
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
Foam::name
word name(const LagrangianState state)
Return a string representation of a Lagrangian state enumeration.
Definition: LagrangianState.C:30
Foam::det
void det(pointPatchField< scalar > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:290
Foam::pow
tmp< DimensionedField< typename powProduct< Type, r >::type, GeoMesh, Field > > pow(const DimensionedField< Type, GeoMesh, PrimitiveField > &df, typename powProduct< Type, r >::type)
Definition: DimensionedFieldFunctions.C:61
Foam::findIndex
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Foam::FatalError
error FatalError
Foam::sqrt
void sqrt(LagrangianPatchField< scalar > &f, const LagrangianPatchField< scalar > &f1)
Definition: LagrangianPatchFieldFunctions.H:88
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(atmosphericBoundaryLayer, 0)
Foam::nl
static const char nl
Definition: Ostream.H:297
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:278
Foam::max
dimensioned< Type > max(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
Definition: DimensionedFieldFunctions.C:373
p
volScalarField & p
Definition: createFieldRefs.H:4