moleculeCloud.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 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_,
58  IOobject::MUST_READ_IF_MODIFIED,
59  IOobject::NO_WRITE,
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  {
81  FatalErrorInFunction
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::nonBlocking);
131 
132  // Start sending referred data
133  label startOfRequests = Pstream::nRequests();
134  il_.sendReferredData(cellOccupancy(), pBufs);
135 
136  molecule* molI = NULL;
137  molecule* molJ = NULL;
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 = NULL;
272  molecule* molJ = NULL;
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::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  {
488  FatalErrorInFunction
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  {
528  FatalErrorInFunction
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  {
550  WarningInFunction
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  {
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  label cell = -1;
750  label tetFace = -1;
751  label tetPt = -1;
752 
753  mesh_.findCellFacePt
754  (
755  globalPosition,
756  cell,
757  tetFace,
758  tetPt
759  );
760 
761  if (findIndex(zone, cell) != -1)
762  {
763  createMolecule
764  (
765  globalPosition,
766  cell,
767  tetFace,
768  tetPt,
769  id,
770  tethered,
771  temperature,
772  bulkVelocity
773  );
774  }
775  }
776  }
777  else
778  {
779  // Place top and bottom caps.
780 
781  labelVector unitCellLatticePosition(0,0,0);
782 
783  for
784  (
785  unitCellLatticePosition.z() = -n;
786  unitCellLatticePosition.z() <= n;
787  unitCellLatticePosition.z() += 2*n
788  )
789  {
790  for
791  (
792  unitCellLatticePosition.y() = -n;
793  unitCellLatticePosition.y() <= n;
794  unitCellLatticePosition.y()++
795  )
796  {
797  for
798  (
799  unitCellLatticePosition.x() = -n;
800  unitCellLatticePosition.x() <= n;
801  unitCellLatticePosition.x()++
802  )
803  {
804  forAll(latticePositions, p)
805  {
806  label id = findIndex
807  (
808  pot_.idList(),
809  latticeIds[p]
810  );
811 
812  const vector& latticePosition =
813  vector
814  (
815  unitCellLatticePosition.x(),
816  unitCellLatticePosition.y(),
817  unitCellLatticePosition.z()
818  )
819  + latticePositions[p];
820 
821  point globalPosition =
822  anchor
823  + (
824  R
825  & (
826  latticeCellShape
827  & latticePosition
828  )
829  );
830 
831  partOfLayerInBounds =
832  mesh_.bounds().contains
833  (
834  globalPosition
835  );
836 
837  label cell = -1;
838  label tetFace = -1;
839  label tetPt = -1;
840 
841  mesh_.findCellFacePt
842  (
843  globalPosition,
844  cell,
845  tetFace,
846  tetPt
847  );
848 
849  if (findIndex(zone, cell) != -1)
850  {
851  createMolecule
852  (
853  globalPosition,
854  cell,
855  tetFace,
856  tetPt,
857  id,
858  tethered,
859  temperature,
860  bulkVelocity
861  );
862  }
863  }
864  }
865  }
866  }
867 
868  for
869  (
870  unitCellLatticePosition.z() = -(n-1);
871  unitCellLatticePosition.z() <= (n-1);
872  unitCellLatticePosition.z()++
873  )
874  {
875  for (label iR = 0; iR <= 2*n -1; iR++)
876  {
877  unitCellLatticePosition.x() = n;
878 
879  unitCellLatticePosition.y() = -n + (iR + 1);
880 
881  for (label iK = 0; iK < 4; iK++)
882  {
883  forAll(latticePositions, p)
884  {
885  label id = findIndex
886  (
887  pot_.idList(),
888  latticeIds[p]
889  );
890 
891  const vector& latticePosition =
892  vector
893  (
894  unitCellLatticePosition.x(),
895  unitCellLatticePosition.y(),
896  unitCellLatticePosition.z()
897  )
898  + latticePositions[p];
899 
900  point globalPosition =
901  anchor
902  + (
903  R
904  & (
905  latticeCellShape
906  & latticePosition
907  )
908  );
909 
910  partOfLayerInBounds =
911  mesh_.bounds().contains
912  (
913  globalPosition
914  );
915 
916  label cell = -1;
917  label tetFace = -1;
918  label tetPt = -1;
919 
920  mesh_.findCellFacePt
921  (
922  globalPosition,
923  cell,
924  tetFace,
925  tetPt
926  );
927 
928  if (findIndex(zone, cell) != -1)
929  {
930  createMolecule
931  (
932  globalPosition,
933  cell,
934  tetFace,
935  tetPt,
936  id,
937  tethered,
938  temperature,
939  bulkVelocity
940  );
941  }
942  }
943 
944  unitCellLatticePosition =
945  labelVector
946  (
947  - unitCellLatticePosition.y(),
948  unitCellLatticePosition.x(),
949  unitCellLatticePosition.z()
950  );
951  }
952  }
953  }
954  }
955 
956  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
957  // end of placement of molecules
958  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
959 
960  if
961  (
962  totalZoneMols == 0
963  && !partOfLayerInBounds
964  )
965  {
966  WarningInFunction
967  << "A whole layer of unit cells was placed "
968  << "outside the bounds of the mesh, but no "
969  << "molecules have been placed in zone '"
970  << zone.name()
971  << "'. This is likely to be because the zone "
972  << "has few cells ("
973  << zone.size()
974  << " in this case) and no lattice position "
975  << "fell inside them. "
976  << "Aborting filling this zone."
977  << endl;
978 
979  break;
980  }
981 
982  molsPlacedThisIteration =
983  this->size() - sizeBeforeIteration;
984 
985  totalZoneMols += molsPlacedThisIteration;
986 
987  n++;
988  }
989  }
990  }
991  }
992 }
993 
994 
995 void Foam::moleculeCloud::createMolecule
996 (
997  const point& position,
998  label cell,
999  label tetFace,
1000  label tetPt,
1001  label id,
1002  bool tethered,
1003  scalar temperature,
1004  const vector& bulkVelocity
1005 )
1006 {
1007  if (cell == -1)
1008  {
1009  mesh_.findCellFacePt(position, cell, tetFace, tetPt);
1010  }
1011 
1012  if (cell == -1)
1013  {
1014  FatalErrorInFunction
1015  << "Position specified does not correspond to a mesh cell." << nl
1016  << abort(FatalError);
1017  }
1018 
1019  point specialPosition(vector::zero);
1020 
1021  label special = 0;
1022 
1023  if (tethered)
1024  {
1025  specialPosition = position;
1026 
1027  special = molecule::SPECIAL_TETHERED;
1028  }
1029 
1030  const molecule::constantProperties& cP(constProps(id));
1031 
1032  vector v = equipartitionLinearVelocity(temperature, cP.mass());
1033 
1034  v += bulkVelocity;
1035 
1036  vector pi = vector::zero;
1037 
1038  tensor Q = I;
1039 
1040  if (!cP.pointMolecule())
1041  {
1042  pi = equipartitionAngularMomentum(temperature, cP);
1043 
1044  scalar phi(rndGen_.scalar01()*twoPi);
1045 
1046  scalar theta(rndGen_.scalar01()*twoPi);
1047 
1048  scalar psi(rndGen_.scalar01()*twoPi);
1049 
1050  Q = tensor
1051  (
1052  cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
1053  cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
1054  sin(psi)*sin(theta),
1055  - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
1056  - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
1057  cos(psi)*sin(theta),
1058  sin(theta)*sin(phi),
1059  - sin(theta)*cos(phi),
1060  cos(theta)
1061  );
1062  }
1063 
1064  addParticle
1065  (
1066  new molecule
1067  (
1068  mesh_,
1069  position,
1070  cell,
1071  tetFace,
1072  tetPt,
1073  Q,
1074  v,
1075  vector::zero,
1076  pi,
1077  vector::zero,
1078  specialPosition,
1079  constProps(id),
1080  special,
1081  id
1082  )
1083  );
1084 }
1085 
1086 
1087 Foam::label Foam::moleculeCloud::nSites() const
1088 {
1089  label n = 0;
1090 
1091  forAllConstIter(moleculeCloud, *this, mol)
1092  {
1093  n += constProps(mol().id()).nSites();
1094  }
1095 
1096  return n;
1097 }
1098 
1099 
1100 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1101 
1102 Foam::moleculeCloud::moleculeCloud
1103 (
1104  const polyMesh& mesh,
1105  const potential& pot,
1106  bool readFields
1107 )
1108 :
1109  Cloud<molecule>(mesh, "moleculeCloud", false),
1110  mesh_(mesh),
1111  pot_(pot),
1112  cellOccupancy_(mesh_.nCells()),
1113  il_(mesh_, pot_.pairPotentials().rCutMax(), false),
1114  constPropList_(),
1115  rndGen_(clock::getTime())
1116 {
1117  if (readFields)
1118  {
1119  molecule::readFields(*this);
1120  }
1121 
1122  buildConstProps();
1123 
1124  setSiteSizesAndPositions();
1125 
1126  removeHighEnergyOverlaps();
1127 
1128  calculateForce();
1129 }
1130 
1131 
1132 Foam::moleculeCloud::moleculeCloud
1133 (
1134  const polyMesh& mesh,
1135  const potential& pot,
1136  const IOdictionary& mdInitialiseDict,
1137  bool readFields
1138 )
1139 :
1140  Cloud<molecule>(mesh, "moleculeCloud", false),
1141  mesh_(mesh),
1142  pot_(pot),
1143  il_(mesh_, 0.0, false),
1144  constPropList_(),
1145  rndGen_(clock::getTime())
1146 {
1147  if (readFields)
1148  {
1149  molecule::readFields(*this);
1150  }
1151 
1152  clear();
1153 
1154  buildConstProps();
1155 
1156  initialiseMolecules(mdInitialiseDict);
1157 }
1158 
1159 
1160 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1161 
1162 void Foam::moleculeCloud::evolve()
1163 {
1164  molecule::trackingData td0(*this, 0);
1165  Cloud<molecule>::move(td0, mesh_.time().deltaTValue());
1166 
1167  molecule::trackingData td1(*this, 1);
1168  Cloud<molecule>::move(td1, mesh_.time().deltaTValue());
1169 
1170  molecule::trackingData td2(*this, 2);
1171  Cloud<molecule>::move(td2, mesh_.time().deltaTValue());
1172 
1173  calculateForce();
1174 
1175  molecule::trackingData td3(*this, 3);
1176  Cloud<molecule>::move(td3, mesh_.time().deltaTValue());
1177 }
1178 
1179 
1180 void Foam::moleculeCloud::calculateForce()
1181 {
1182  buildCellOccupancy();
1183 
1184  // Set accumulated quantities to zero
1185  forAllIter(moleculeCloud, *this, mol)
1186  {
1187  mol().siteForces() = vector::zero;
1188 
1189  mol().potentialEnergy() = 0.0;
1190 
1191  mol().rf() = tensor::zero;
1192  }
1193 
1194  calculatePairForce();
1195 
1196  calculateTetherForce();
1197 
1198  calculateExternalForce();
1199 }
1200 
1201 
1202 void Foam::moleculeCloud::applyConstraintsAndThermostats
1203 (
1204  const scalar targetTemperature,
1205  const scalar measuredTemperature
1206 )
1207 {
1208  scalar temperatureCorrectionFactor =
1209  sqrt(targetTemperature/max(VSMALL, measuredTemperature));
1210 
1211  Info<< "----------------------------------------" << nl
1212  << "Temperature equilibration" << nl
1213  << "Target temperature = "
1214  << targetTemperature << nl
1215  << "Measured temperature = "
1216  << measuredTemperature << nl
1217  << "Temperature correction factor = "
1218  << temperatureCorrectionFactor << nl
1219  << "----------------------------------------"
1220  << endl;
1221 
1222  forAllIter(moleculeCloud, *this, mol)
1223  {
1224  mol().v() *= temperatureCorrectionFactor;
1225 
1226  mol().pi() *= temperatureCorrectionFactor;
1227  }
1228 }
1229 
1230 
1231 void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
1232 {
1233  OFstream os(fName);
1234 
1235  os << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
1236 
1237  forAllConstIter(moleculeCloud, *this, mol)
1238  {
1239  const molecule::constantProperties& cP = constProps(mol().id());
1240 
1241  forAll(mol().sitePositions(), i)
1242  {
1243  const point& sP = mol().sitePositions()[i];
1244 
1245  os << pot_.siteIdList()[cP.siteIds()[i]]
1246  << ' ' << sP.x()*1e10
1247  << ' ' << sP.y()*1e10
1248  << ' ' << sP.z()*1e10
1249  << nl;
1250  }
1251  }
1252 }
1253 
1254 
1255 // ************************************************************************* //
Foam::moleculeCloud::buildCellOccupancy
void buildCellOccupancy()
Determine which molecules are in which cells.
Definition: moleculeCloud.C:109
Foam::IOdictionary
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
Foam::moleculeCloud::setSiteSizesAndPositions
void setSiteSizesAndPositions()
Definition: moleculeCloud.C:96
Foam::Tensor
Templated 3D tensor derived from VectorSpace adding construction from 9 components,...
Definition: complexI.H:224
Foam::tetherPotentialList::force
vector force(const label a, const vector rIT) const
Definition: tetherPotentialList.C:139
Foam::IOobject
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
mathematicalConstants.H
Foam::DiagTensor< scalar >
Foam::Vector< scalar >::zero
static const Vector zero
Definition: Vector.H:80
p
p
Definition: pEqn.H:62
Foam::word
A class for handling words, derived from string.
Definition: word.H:59
Foam::particle::origId
label origId() const
Return const access to the particle id on originating processor.
Definition: particleI.H:840
forAllIter
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:431
Foam::fileName
A class for handling file names.
Definition: fileName.H:69
Foam::moleculeCloud::nSites
label nSites() const
Definition: moleculeCloud.C:1087
clear
UEqn clear()
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
Foam::findIndex
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:56
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::tetherPotentialList::energy
scalar energy(const label a, const vector rIT) const
Definition: tetherPotentialList.C:149
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:255
Foam::molecule::constantProperties::siteIds
const List< label > & siteIds() const
Definition: moleculeI.H:352
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Foam::PstreamBuffers
Buffers for inter-processor communications streams (UOPstream, UIPstream).
Definition: PstreamBuffers.H:85
Foam::UPstream::parRun
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:377
Foam::dictionary::lookup
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:449
Foam::zone
Base class for zones.
Definition: zone.H:57
Foam::molecule::id
label id() const
Definition: moleculeI.H:598
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Foam::constant::mathematical::twoPi
const scalar twoPi(2 *pi)
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:179
Foam::readFields
This function object reads fields from the time directories and adds them to the mesh database for fu...
Definition: readFields.H:104
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Foam::IOobject::NO_WRITE
@ NO_WRITE
Definition: IOobject.H:118
Foam::cellZone
A subset of mesh cells.
Definition: cellZone.H:61
forAllConstIter
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
Foam::labelVector
Vector< label > labelVector
Vector of labels.
Definition: labelVector.H:49
Foam::Q
This function object calculates and outputs the second invariant of the velocity gradient tensor [1/s...
Definition: Q.H:119
Foam::molecule::constantProperties::pointMolecule
bool pointMolecule() const
Definition: moleculeI.H:428
Foam::IOobject::MUST_READ_IF_MODIFIED
@ MUST_READ_IF_MODIFIED
Definition: IOobject.H:109
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:43
moleculeCloud.H
R
#define R(A, B, C, D, E, F, K, M)
Foam::UPstream::nonBlocking
@ nonBlocking
Definition: UPstream.H:68
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::moleculeCloud::calculateTetherForce
void calculateTetherForce()
Definition: moleculeCloud.C:216
Foam::dictionary::found
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:304
Foam::molecule
Foam::molecule.
Definition: molecule.H:56
Foam::IDLList
Intrusive doubly-linked list.
Definition: IDLList.H:47
cellOccupancy
const List< DynamicList< molecule * > > & cellOccupancy
Definition: calculateMDFields.H:1
Foam::Vector< scalar >::one
static const Vector one
Definition: Vector.H:81
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:71
Foam::nl
static const char nl
Definition: Ostream.H:260
Foam::Info
messageStream Info
Foam::clock::getTime
static time_t getTime()
Get the current clock time in seconds.
Definition: clock.C:43
Foam::moleculeCloud::buildConstProps
void buildConstProps()
Definition: moleculeCloud.C:41
Foam::molecule::trackingData
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:157
Foam::ZoneMesh< cellZone, polyMesh >
Foam::moleculeCloud::evolve
void evolve()
Evolve the molecules (move, calculate forces, control state etc)
Definition: moleculeCloud.C:1162
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:73
Foam::moleculeCloud::calculatePairForce
void calculatePairForce()
Definition: moleculeCloud.C:128
Foam::I
static const sphericalTensor I(1)
Foam::moleculeCloud::writeXYZ
void writeXYZ(const fileName &fName) const
Write molecule sites in XYZ format.
Definition: moleculeCloud.C:1231
Foam::Vector::x
const Cmpt & x() const
Definition: VectorI.H:65
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:137
Foam::potential
Definition: potential.H:53
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
Foam::Tensor::zero
static const Tensor zero
Definition: Tensor.H:80
fvMesh.H
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Foam::Vector::z
const Cmpt & z() const
Definition: VectorI.H:77
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::molecule::constantProperties::nSites
label nSites() const
Definition: moleculeI.H:457
Foam::zone::name
const word & name() const
Return name.
Definition: zone.H:150
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Foam::DynamicList::append
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Foam::OFstream
Output to file stream.
Definition: OFstream.H:81
Foam::moleculeCloud::initialiseMolecules
void initialiseMolecules(const IOdictionary &mdInitialiseDict)
Definition: moleculeCloud.C:476
constProp
word constProp(initialConditions.lookup("constantProperty"))
Foam::max
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Foam::molecule::readFields
static void readFields(Cloud< molecule > &mC)
Definition: moleculeIO.C:94
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:318
Foam::constant::mathematical
mathematical constants.
psi
const volScalarField & psi
Definition: setRegionFluidFields.H:13
Foam::UPstream::nRequests
static label nRequests()
Get number of outstanding requests.
Definition: UPstream.C:96
Foam::sumOp
Definition: ops.H:162
Foam::tab
static const char tab
Definition: Ostream.H:259
Foam::moleculeCloud::calculateExternalForce
void calculateExternalForce()
Definition: moleculeCloud.C:243
Foam::readScalar
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
Foam::molecule::constantProperties::mass
scalar mass() const
Definition: moleculeI.H:451
Foam::tetherPotentialList
Definition: tetherPotentialList.H:50
Foam::moleculeCloud
Definition: moleculeCloud.H:56
Foam::Vector< scalar >
Foam::moleculeCloud::applyConstraintsAndThermostats
void applyConstraintsAndThermostats(const scalar targetTemperature, const scalar measuredTemperature)
Definition: moleculeCloud.C:1203
Foam::List< word >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:142
Foam::moleculeCloud::createMolecule
void createMolecule(const point &position, label cell, label tetFace, label tetPt, label id, bool tethered, scalar temperature, const vector &bulkVelocity)
Definition: moleculeCloud.C:996
Foam::Cloud< molecule >
Foam::molecule::constantProperties
Class to hold molecule constant properties.
Definition: molecule.H:81
Foam::constant::mathematical::pi
const scalar pi(M_PI)
Foam::Vector::y
const Cmpt & y() const
Definition: VectorI.H:71
Foam::moleculeCloud::calculateForce
void calculateForce()
Definition: moleculeCloud.C:1180
Foam::molecule::SPECIAL_TETHERED
@ SPECIAL_TETHERED
Definition: molecule.H:74
Foam::det
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:60
Foam::List::size
void size(const label)
Override size to be inconsistent with allocated storage.
Foam::defineTemplateTypeNameAndDebug
defineTemplateTypeNameAndDebug(IOPtrList< ensightPart >, 0)
Foam::dictionary::subDict
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:631
Foam::moleculeCloud::removeHighEnergyOverlaps
void removeHighEnergyOverlaps()
Definition: moleculeCloud.C:252
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:259
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
Foam::moleculeCloud::moleculeCloud
moleculeCloud(const moleculeCloud &)
Disallow default bitwise copy construct.
Foam::Cloud::move
void move(TrackData &td, const scalar trackTime)
Move the particles.
Definition: Cloud.C:187
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:256

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17