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-2018 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  defineTemplateTypeNameAndDebug(Cloud<molecule>, 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 cellZoneMesh& 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  readScalar(zoneDict.lookup("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 = readScalar
544  (
545  zoneDict.lookup("numberDensity")
546  );
547 
548  if (numberDensity < vSmall)
549  {
551  << "numberDensity too small, not filling zone "
552  << zone.name() << endl;
553 
554  continue;
555  }
556 
557  latticeCellScale = pow
558  (
559  latticeIds.size()/(det(latticeCellShape)*numberDensity),
560  (1.0/3.0)
561  );
562  }
563  else if (zoneDict.found("massDensity"))
564  {
565  scalar unitCellMass = 0.0;
566 
567  forAll(latticeIds, i)
568  {
569  label id = findIndex(pot_.idList(), latticeIds[i]);
570 
571  const molecule::constantProperties& cP(constProps(id));
572 
573  unitCellMass += cP.mass();
574  }
575 
576  scalar massDensity = readScalar
577  (
578  zoneDict.lookup("massDensity")
579  );
580 
581  if (massDensity < vSmall)
582  {
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  {
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  mesh_.cellTree().findInside(globalPosition);
751 
752  if (findIndex(zone, cell) != -1)
753  {
754  createMolecule
755  (
756  globalPosition,
757  cell,
758  id,
759  tethered,
760  temperature,
761  bulkVelocity
762  );
763  }
764  }
765  }
766  else
767  {
768  // Place top and bottom caps.
769 
770  labelVector unitCellLatticePosition(0,0,0);
771 
772  for
773  (
774  unitCellLatticePosition.z() = -n;
775  unitCellLatticePosition.z() <= n;
776  unitCellLatticePosition.z() += 2*n
777  )
778  {
779  for
780  (
781  unitCellLatticePosition.y() = -n;
782  unitCellLatticePosition.y() <= n;
783  unitCellLatticePosition.y()++
784  )
785  {
786  for
787  (
788  unitCellLatticePosition.x() = -n;
789  unitCellLatticePosition.x() <= n;
790  unitCellLatticePosition.x()++
791  )
792  {
793  forAll(latticePositions, p)
794  {
795  label id = findIndex
796  (
797  pot_.idList(),
798  latticeIds[p]
799  );
800 
801  const vector& latticePosition =
802  vector
803  (
804  unitCellLatticePosition.x(),
805  unitCellLatticePosition.y(),
806  unitCellLatticePosition.z()
807  )
808  + latticePositions[p];
809 
810  point globalPosition =
811  anchor
812  + (
813  R
814  & (
815  latticeCellShape
816  & latticePosition
817  )
818  );
819 
820  partOfLayerInBounds =
821  mesh_.bounds().contains
822  (
823  globalPosition
824  );
825 
826  const label cell =
827  mesh_.cellTree().findInside
828  (
829  globalPosition
830  );
831 
832  if (findIndex(zone, cell) != -1)
833  {
834  createMolecule
835  (
836  globalPosition,
837  cell,
838  id,
839  tethered,
840  temperature,
841  bulkVelocity
842  );
843  }
844  }
845  }
846  }
847  }
848 
849  for
850  (
851  unitCellLatticePosition.z() = -(n-1);
852  unitCellLatticePosition.z() <= (n-1);
853  unitCellLatticePosition.z()++
854  )
855  {
856  for (label iR = 0; iR <= 2*n -1; iR++)
857  {
858  unitCellLatticePosition.x() = n;
859 
860  unitCellLatticePosition.y() = -n + (iR + 1);
861 
862  for (label iK = 0; iK < 4; iK++)
863  {
864  forAll(latticePositions, p)
865  {
866  label id = findIndex
867  (
868  pot_.idList(),
869  latticeIds[p]
870  );
871 
872  const vector& latticePosition =
873  vector
874  (
875  unitCellLatticePosition.x(),
876  unitCellLatticePosition.y(),
877  unitCellLatticePosition.z()
878  )
879  + latticePositions[p];
880 
881  point globalPosition =
882  anchor
883  + (
884  R
885  & (
886  latticeCellShape
887  & latticePosition
888  )
889  );
890 
891  partOfLayerInBounds =
892  mesh_.bounds().contains
893  (
894  globalPosition
895  );
896 
897  const label cell =
898  mesh_.cellTree().findInside
899  (
900  globalPosition
901  );
902 
903  if (findIndex(zone, cell) != -1)
904  {
905  createMolecule
906  (
907  globalPosition,
908  cell,
909  id,
910  tethered,
911  temperature,
912  bulkVelocity
913  );
914  }
915  }
916 
917  unitCellLatticePosition =
919  (
920  - unitCellLatticePosition.y(),
921  unitCellLatticePosition.x(),
922  unitCellLatticePosition.z()
923  );
924  }
925  }
926  }
927  }
928 
929  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
930  // end of placement of molecules
931  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
932 
933  if
934  (
935  totalZoneMols == 0
936  && !partOfLayerInBounds
937  )
938  {
940  << "A whole layer of unit cells was placed "
941  << "outside the bounds of the mesh, but no "
942  << "molecules have been placed in zone '"
943  << zone.name()
944  << "'. This is likely to be because the zone "
945  << "has few cells ("
946  << zone.size()
947  << " in this case) and no lattice position "
948  << "fell inside them. "
949  << "Aborting filling this zone."
950  << endl;
951 
952  break;
953  }
954 
955  molsPlacedThisIteration =
956  this->size() - sizeBeforeIteration;
957 
958  totalZoneMols += molsPlacedThisIteration;
959 
960  n++;
961  }
962  }
963  }
964  }
965 }
966 
967 
968 void Foam::moleculeCloud::createMolecule
969 (
970  const point& position,
971  label cell,
972  label id,
973  bool tethered,
974  scalar temperature,
975  const vector& bulkVelocity
976 )
977 {
978  point specialPosition(Zero);
979 
980  label special = 0;
981 
982  if (tethered)
983  {
984  specialPosition = position;
985 
986  special = molecule::SPECIAL_TETHERED;
987  }
988 
989  const molecule::constantProperties& cP(constProps(id));
990 
991  vector v = equipartitionLinearVelocity(temperature, cP.mass());
992 
993  v += bulkVelocity;
994 
995  vector pi = Zero;
996 
997  tensor Q = I;
998 
999  if (!cP.pointMolecule())
1000  {
1001  pi = equipartitionAngularMomentum(temperature, cP);
1002 
1003  scalar phi(rndGen_.scalar01()*twoPi);
1004 
1005  scalar theta(rndGen_.scalar01()*twoPi);
1006 
1007  scalar psi(rndGen_.scalar01()*twoPi);
1008 
1009  Q = tensor
1010  (
1011  cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
1012  cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
1013  sin(psi)*sin(theta),
1014  - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
1015  - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
1016  cos(psi)*sin(theta),
1017  sin(theta)*sin(phi),
1018  - sin(theta)*cos(phi),
1019  cos(theta)
1020  );
1021  }
1022 
1023  addParticle
1024  (
1025  new molecule
1026  (
1027  mesh_,
1028  position,
1029  cell,
1030  Q,
1031  v,
1032  Zero,
1033  pi,
1034  Zero,
1035  specialPosition,
1036  constProps(id),
1037  special,
1038  id
1039  )
1040  );
1041 }
1042 
1043 
1044 Foam::label Foam::moleculeCloud::nSites() const
1045 {
1046  label n = 0;
1047 
1048  forAllConstIter(moleculeCloud, *this, mol)
1049  {
1050  n += constProps(mol().id()).nSites();
1051  }
1052 
1053  return n;
1054 }
1055 
1056 
1057 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1058 
1061  const polyMesh& mesh,
1062  const potential& pot,
1063  bool readFields
1064 )
1065 :
1066  Cloud<molecule>(mesh, "moleculeCloud", false),
1067  mesh_(mesh),
1068  pot_(pot),
1069  cellOccupancy_(mesh_.nCells()),
1070  il_(mesh_, pot_.pairPotentials().rCutMax(), false),
1071  constPropList_(),
1072  rndGen_(clock::getTime())
1073 {
1074  if (readFields)
1075  {
1076  molecule::readFields(*this);
1077  }
1078 
1079  buildConstProps();
1080 
1081  setSiteSizesAndPositions();
1082 
1083  removeHighEnergyOverlaps();
1084 
1085  calculateForce();
1086 }
1087 
1088 
1091  const polyMesh& mesh,
1092  const potential& pot,
1093  const IOdictionary& mdInitialiseDict,
1094  bool readFields
1095 )
1096 :
1097  Cloud<molecule>(mesh, "moleculeCloud", false),
1098  mesh_(mesh),
1099  pot_(pot),
1100  il_(mesh_, 0.0, false),
1101  constPropList_(),
1102  rndGen_(clock::getTime())
1103 {
1104  if (readFields)
1105  {
1106  molecule::readFields(*this);
1107  }
1108 
1109  clear();
1110 
1111  buildConstProps();
1112 
1113  initialiseMolecules(mdInitialiseDict);
1114 }
1115 
1116 
1117 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1118 
1120 {
1121  molecule::trackingData td0(*this, 0);
1122  Cloud<molecule>::move(*this, td0, mesh_.time().deltaTValue());
1123 
1124  molecule::trackingData td1(*this, 1);
1125  Cloud<molecule>::move(*this, td1, mesh_.time().deltaTValue());
1126 
1127  molecule::trackingData td2(*this, 2);
1128  Cloud<molecule>::move(*this, td2, mesh_.time().deltaTValue());
1129 
1130  calculateForce();
1131 
1132  molecule::trackingData td3(*this, 3);
1133  Cloud<molecule>::move(*this, td3, mesh_.time().deltaTValue());
1134 }
1135 
1136 
1138 {
1139  buildCellOccupancy();
1140 
1141  // Set accumulated quantities to zero
1142  forAllIter(moleculeCloud, *this, mol)
1143  {
1144  mol().siteForces() = Zero;
1145 
1146  mol().potentialEnergy() = 0.0;
1147 
1148  mol().rf() = Zero;
1149  }
1150 
1151  calculatePairForce();
1152 
1153  calculateTetherForce();
1154 
1155  calculateExternalForce();
1156 }
1157 
1158 
1161  const scalar targetTemperature,
1162  const scalar measuredTemperature
1163 )
1164 {
1165  scalar temperatureCorrectionFactor =
1166  sqrt(targetTemperature/max(vSmall, measuredTemperature));
1167 
1168  Info<< "----------------------------------------" << nl
1169  << "Temperature equilibration" << nl
1170  << "Target temperature = "
1171  << targetTemperature << nl
1172  << "Measured temperature = "
1173  << measuredTemperature << nl
1174  << "Temperature correction factor = "
1175  << temperatureCorrectionFactor << nl
1176  << "----------------------------------------"
1177  << endl;
1178 
1179  forAllIter(moleculeCloud, *this, mol)
1180  {
1181  mol().v() *= temperatureCorrectionFactor;
1182 
1183  mol().pi() *= temperatureCorrectionFactor;
1184  }
1185 }
1186 
1187 
1188 void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
1189 {
1190  OFstream os(fName);
1191 
1192  os << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
1193 
1194  forAllConstIter(moleculeCloud, *this, mol)
1195  {
1196  const molecule::constantProperties& cP = constProps(mol().id());
1197 
1198  forAll(mol().sitePositions(), i)
1199  {
1200  const point& sP = mol().sitePositions()[i];
1201 
1202  os << pot_.siteIdList()[cP.siteIds()[i]]
1203  << ' ' << sP.x()*1e10
1204  << ' ' << sP.y()*1e10
1205  << ' ' << sP.z()*1e10
1206  << nl;
1207  }
1208  }
1209 }
1210 
1211 
1212 // ************************************************************************* //
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
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
surfaceScalarField & phi
Class to hold molecule constant properties.
Definition: molecule.H:89
static const char tab
Definition: Ostream.H:264
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
tUEqn clear()
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:256
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:52
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)
static const zero Zero
Definition: zero.H:97
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
defineTemplateTypeNameAndDebug(IOPtrList< ensightPart >, 0)
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:265
const word constProp(initialConditions.lookup("constantProperty"))
const Time & time() const
Return time.
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:165
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
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
Namespace for OpenFOAM.