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