MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PeriodicBoundaryHandler.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
27 #include <iostream>
28 #include <iomanip>
29 #include <vector>
30 #include "Math/Helpers.h"
31 
34 #include "DPMBase.h"
35 #include "ParticleHandler.h"
36 #include <set>
38 #include "Logger.h"
40 
43 {
45  logger(DEBUG, "PeriodicBoundaryHandler::PeriodicBoundaryHandler() finished");
46 }
47 
54 {
55  objects_.clear();
56  setDPMBase(nullptr);
58  logger(DEBUG, "PeriodicBoundaryHandler::PeriodicBoundaryHandler(const PeriodicBoundaryHandler &BH) finished");
59 }
60 
67 {
68  if (this != &rhs)
69  {
70  objects_.clear();
72  }
73  logger(DEBUG, "PeriodicBoundaryHandler PeriodicBoundaryHandler::operator =(const BoundaryHandler& rhs)");
74  return *this;
75 }
76 
79 {
83  objects_.clear();
84 
85  logger(DEBUG, "PeriodicBoundaryHandler::~PeriodicBoundaryHandler() finished");
86 }
87 
91 {
92 #ifdef MERCURY_USE_MPI
93  if (NUMBER_OF_PROCESSORS == 1)
94  {
95  return;
96  }
97  //Puts the particle in the Particle list
98  objects_.push_back(P);
99  //set the particleHandler pointer
100  P->setPeriodicHandler(this);
101  //Remove all ghost particles
103  //Initialise ghost particles
104  addNewParticles();
105 #endif
106 }
107 
110 {
111 }
112 
118 {
119  return "PeriodicBoundaryHandler";
120 }
121 
129 {
130  interactionDistance_ = interactionDistance;
131 }
132 
138 {
139  return interactionDistance_;
140 }
141 
151 void PeriodicBoundaryHandler::updateStatus(std::set<BaseParticle*>& particlesToBeDeleted)
152 {
153  //Step 1: collect position and velocity data of the real particle and send the data
155 
156  //Step 2: finalise the data transision
158 
159  //Step 3: Update the status of the ghost and real particles
160  updateParticleStatus(particlesToBeDeleted);
161 }
162 
171 {
172  shiftParticle(particle, particle->getPeriodicComplexity());
173 }
174 
183 void PeriodicBoundaryHandler::shiftParticle(BaseParticle* particle, const std::vector<int>& complexity)
184 {
185  int boundaryIndex = 0;
186  for (BasePeriodicBoundary* boundary : *this)
187  {
188  if (complexity[boundaryIndex] == -1)
189  {
190  boundary->shiftPosition(particle);
191  }
192  boundaryIndex++;
193  }
194 }
195 
210 void
211 PeriodicBoundaryHandler::computePeriodicComplexity(std::vector<int>& periodicComplexity, int& totalPeriodicComplexity,
212  Vec3D position)
213 {
214  //Compute the new position
215  periodicComplexity.resize(this->getSize());
216  totalPeriodicComplexity = 0;
217  int index = 0;
218  for (BasePeriodicBoundary* boundary : *this)
219  {
220  Mdouble distance = boundary->getDistance(position);
221  if (std::abs(distance) <= interactionDistance_)
222  {
223  if (distance > 0)
224  {
225  //real particle
226  periodicComplexity[index] = 1;
227  totalPeriodicComplexity++;
228  }
229  else if (distance < 0)
230  {
231  //ghost particle
232  periodicComplexity[index] = -1;
233  }
234  else //distance == 0 In extremely rare cases (test cases mostly)
235  {
236  //Check if on left or right side of boundary
237  if (boundary->isClosestToLeftBoundary(position))
238  {
239  //ghost particle on the left side
240  periodicComplexity[index] = -1;
241  }
242  else
243  {
244  //real particle on the right side
245  periodicComplexity[index] = 1;
246  totalPeriodicComplexity++;
247  }
248  }
249  }
250  else
251  {
252  if (distance < 0)
253  {
254  periodicComplexity[index] = -2;
255  }
256  else
257  {
258  periodicComplexity[index] = 2;
259  }
260  }
261  index++;
262  }
263 /*
264  //Modify periodic complexity in case of maser
265  for (int b = 0; b < getSize(); b++)
266  {
267  objects_[b]->modifyPeriodicComplexity(periodicComplexity, position, b);
268  }
269 */
270 }
271 
283 {
284  std::vector<int> periodicComplexity;
285  int totalPeriodicComplexity;
286 
287  computePeriodicComplexity(periodicComplexity, totalPeriodicComplexity, position);
288 
289  return periodicComplexity;
290 }
291 
301 {
302 #ifdef MERCURY_USE_MPI
303  if (NUMBER_OF_PROCESSORS == 1)
304  {
305  return;
306  }
307  //Step 0: Perform actions before finding new particles. I.e. opening maser boundary
309 
310  //Step 1: Find new particles, such that the target domains can be found
312 
313  //Step 3: Communicate number of particles/interactions to other dommains
315 
316  //Step 4: Perform data transmission
318 
319  //Step 5: Finalise data transmission
321 #endif
322 }
323 
334 {
335 #ifdef MERCURY_USE_MPI
336  if (NUMBER_OF_PROCESSORS == 1)
337  {
338  return;
339  }
340  if (getDPMBase()->getCurrentDomain()->containsParticle(particle))
341  {
342  //Step 1: Find if the new particle needs to be added
343  findNewParticle(particle);
344  }
345  //Follow the standard communication routine
349 #endif
350 }
351 
358 {
359  unsigned int sum = 0;
360  for (auto& index : periodicGhostList_)
361  {
362  sum += index.size();
363  }
364  return sum;
365 }
366 
374 {
375  int sum = 0;
376  for (auto& index : periodicGhostList_)
377  {
378  int numberOfMPIParticles = 0;
379  for (int pIndex = 0; pIndex < index.size(); pIndex++)
380  {
381  if (index[pIndex]->particle->isMPIParticle() == true)
382  {
383  numberOfMPIParticles++;
384  }
385  }
386  sum += (index.size() - numberOfMPIParticles);
387  }
388  return sum;
389 }
390 
401 void PeriodicBoundaryHandler::getMPIFlags(BaseParticle* particle, bool& isInMPIDomain, bool& isMPIParticle)
402 {
403  //Check if the particle is actually in this domain. If not the case it is an "MG"-particle
405  if (domain->containsParticle(particle))
406  {
407  //If the particle is in the inner domain, it is not in the mpi doamin
408  if (domain->isInInnerDomain(particle))
409  {
410  isInMPIDomain = false;
411  isMPIParticle = false;
412  }
413  else
414  {
415  isInMPIDomain = true;
416  isMPIParticle = false;
417  }
418  }
419  else
420  {
421  //If the particle is outside the communication zone it is not in the mpi domain
422  if (domain->isInGreaterDomain(particle))
423  {
424  isInMPIDomain = true;
425  isMPIParticle = true;
426  }
427  else
428  {
429  isInMPIDomain = false;
430  isMPIParticle = true;
431  }
432  }
433 }
434 
444 {
445  bool isInMPIDomain;
446  bool isMPIParticle;
447  getMPIFlags(particle, isInMPIDomain, isMPIParticle);
448  particle->setInMPIDomain(isInMPIDomain);
449  particle->setMPIParticle(isMPIParticle);
450 }
451 
464 void
465 PeriodicBoundaryHandler::generateGhosts(std::vector<std::vector<int> >& list, const std::vector<int> periodicComplexity,
466  std::vector<int>& complexity, int level)
467 {
468  //Add the complexity to the list
469  if (level == 0)
470  {
471  //logger(VERBOSE,"complexity: [%, %, %]",complexity[0], complexity[1], complexity[2]);
472  list.push_back(complexity);
473  }
474  else
475  {
476  //Vary complexity values of level
477  if (periodicComplexity[level - 1] == 1)
478  {
479  complexity[level - 1] = 1;
480  generateGhosts(list, periodicComplexity, complexity, level - 1);
481 
482  complexity[level - 1] = -1;
483  generateGhosts(list, periodicComplexity, complexity, level - 1);
484  }
485  else
486  {
487  generateGhosts(list, periodicComplexity, complexity, level - 1);
488  }
489  }
490 }
491 
498 {
499  //For all new target lists
500  unsigned long numberOfTargets = sendTargetList_.size();
501  periodicGhostParticleSend_.resize(numberOfTargets);
502  periodicGhostComplexitySend_.resize(numberOfTargets);
503  for (int i = 0; i < numberOfTargets; i++)
504  {
506  {
507  //For all new ghost particles copy particle data and periodic complexity
508  for (int j = 0; j < numberOfNewPeriodicGhostParticlesSend_[i]; j++)
509  {
511  BaseParticle* particle = ppid->particle;
513  for (int k = 0; k < getSize(); k++)
514  {
516  }
517  }
518  }
519  }
520 }
521 
529 {
530  //For all send target lists
532  unsigned long numberOfTargets = sendTargetList_.size();
533  interactionDataSend_.resize(numberOfTargets);
534  for (int i = 0; i < numberOfTargets; i++)
535  {
536  //Allocate enough space
538 
539  //Fill in the vector
540  unsigned int indexInteraction = 0;
541  for (BaseInteraction* interaction : newInteractionList_[i])
542  {
543  interaction->getMPIInteraction(interactionDataSend_[i], indexInteraction);
544  indexInteraction++;
545  }
546  }
547 }
548 
559 void
560 PeriodicBoundaryHandler::processReceivedGhostParticleData(int targetIndex, std::vector<BaseParticle*>& newParticles)
561 {
562  for (int j = 0; j < numberOfNewPeriodicGhostParticlesReceive_[targetIndex]; j++)
563  {
564  //Create the ghost particle and copy basic information
567  particle, &(getDPMBase()->particleHandler));
568 
569  //Obtain real periodic complexity
570  std::vector<int> realPeriodicComplexity = computePeriodicComplexity(particle->getPosition());
571 
572  //Obtain and set the ghost periodic complexity
573  std::vector<int> ghostPeriodicComplexity(getSize());
574  for (int k = 0; k < getSize(); k++)
575  {
576  ghostPeriodicComplexity[k] = periodicGhostComplexityReceive_[targetIndex][getSize() * j + k];
577  }
578  particle->setPeriodicComplexity(ghostPeriodicComplexity);
579 
580  //Shift the ghost particle to it's correct positions and velocities
581  shiftParticle(particle, ghostPeriodicComplexity);
582 
583  //Add particle to simulation
584  logger(VERBOSE, "Adding a ghost at position %", particle->getPosition());
586 
587  //Set the correct flags
589  pGhost->setPeriodicGhostParticle(true);
590  pGhost->setInPeriodicDomain(true);
591  //Give the correct mpi flags
592  setMPIFlags(pGhost);
593 
594  //Create the periodic ID
596  gpid->realPeriodicComplexity = realPeriodicComplexity;
597  gpid->particle = pGhost;
598  periodicGhostList_[receiveTargetList_[targetIndex]].push_back(gpid);
599 
600  //Add to the newParticle list used for possible interactions
601  newParticles.push_back(pGhost);
602  }
603 }
604 
615 void PeriodicBoundaryHandler::processReceivedInteractionData(int targetIndex, std::vector<BaseParticle*>& newParticles)
616 {
618  for (int l = 0; l < numberOfNewInteractionsReceive_[targetIndex]; l++)
619  {
620  unsigned int identificationP;
621  unsigned int identificationI;
622  bool isWallInteraction;
623  unsigned timeStamp;
624 
625  //Get the general information required to setup a new interaction
626  iH.getInteractionDetails(interactionDataReceive_[targetIndex], l, identificationP, identificationI,
627  isWallInteraction, timeStamp);
628 
629  //Obtain the particle pointer of the ghost
630  BaseParticle* pGhost;
631  int idOther;
632  //Check if the particle is P
633  for (BaseParticle* particle : newParticles)
634  {
635  if (particle->getId() == identificationP)
636  {
637  pGhost = particle;
638  idOther = identificationI;
639  break;
640  }
641 
642  if (particle->getId() == identificationI)
643  {
644  pGhost = particle;
645  idOther = identificationP;
646  break;
647  }
648  }
649 
650  //If it is a wall interaction, do stuff
651  if (isWallInteraction)
652  {
653  BaseInteractable* I = getDPMBase()->wallHandler.getObjectById(identificationI);
654  //Create interactions
655  BaseInteraction* i = I->getInteractionWith(pGhost, timeStamp, &iH);
656  if (i!= nullptr) i->setMPIInteraction(interactionDataReceive_[targetIndex], l, false);
657  }
658  else
659  {
660  //Obtain potential interaction particles
661  std::vector<BaseParticle*> interactingParticleList;
662  getDPMBase()->hGridGetInteractingParticleList(pGhost, interactingParticleList);
663 
664  //Find the other interacting particle
665  BaseParticle* otherParticle;
666  for (BaseParticle* p2 : interactingParticleList)
667  {
668  if (p2->getId() == idOther)
669  {
670  otherParticle = p2;
671  break;
672  }
673  }
674 
675  //Add the interaction
676  BaseInteraction* i = pGhost->getInteractionWith(otherParticle, timeStamp, &iH);
677  if (i!= nullptr) i->setMPIInteraction(interactionDataReceive_[targetIndex], l, false);
678  }
679  }
680 }
681 
692 void PeriodicBoundaryHandler::processLocalInteractionData(std::vector<BaseParticle*>& newParticles)
693 {
695  for (int i = 0; i < sendTargetList_.size(); i++)
696  {
698  {
699  for (int l = 0; l < numberOfNewInteractionsSend_[i]; l++)
700  {
701  unsigned int identificationP;
702  unsigned int identificationI;
703  bool isWallInteraction;
704  unsigned timeStamp;
705 
706  //Get the general information required to setup a new interaction
707  iH.getInteractionDetails(interactionDataSend_[i], l, identificationP, identificationI,
708  isWallInteraction, timeStamp);
709 
710  //Obtain the particle pointer of the ghost
711  BaseParticle* pGhost;
712  int idOther;
713  //Check if the particle is P
714  for (BaseParticle* particle : newParticles)
715  {
716  if (particle->getId() == identificationP)
717  {
718  pGhost = particle;
719  idOther = identificationI;
720  break;
721  }
722 
723  if (particle->getId() == identificationI)
724  {
725  pGhost = particle;
726  idOther = identificationP;
727  break;
728  }
729  }
730 
731  //If it is a wall interaction, do stuff
732  if (isWallInteraction)
733  {
734  BaseWall* I = getDPMBase()->wallHandler.getObjectById(identificationI);
735  //Create interactions
736  BaseInteraction* j = I->getInteractionWith(pGhost, timeStamp, &iH);
737  if (j!= nullptr) j->setMPIInteraction(interactionDataSend_[i], l, false);
738  logger(VERBOSE, "Wall interaction added!");
739  }
740  else
741  {
742  //Obtain potential interaction particles
743  std::vector<BaseParticle*> interactingParticleList;
744  getDPMBase()->hGridGetInteractingParticleList(pGhost, interactingParticleList);
745 
746  if (interactingParticleList.empty())
747  {
748  logger(VERBOSE, "Failed in creating an interaction :(");
749  }
750 
751  //Find the other interacting particle
752  BaseParticle* otherParticle = nullptr;
753  for (BaseParticle* p2 : interactingParticleList)
754  {
755  if (p2->getId() == idOther)
756  {
757  otherParticle = p2;
758  break;
759  }
760  }
761  if (otherParticle == nullptr)
762  {
763  //The interacting object can't be found
764  continue;
765  }
766 
767  //Add the interaction
768  BaseInteraction* j = pGhost->getInteractionWith(otherParticle, timeStamp, &iH);
769  if (j!= nullptr) j->setMPIInteraction(interactionDataSend_[i], l, false);
770  }
771  }
772  }
773  }
774 }
775 
776 
783 {
784  for (int i : sendTargetList_)
785  {
786  for (int j = 0; j < newPeriodicParticleList_[i].size(); j++)
787  {
788  //Update the particle status
790  ppid->particle->setInPeriodicDomain(true);
791 
792  //make new entry in the list
793  periodicParticleList_[i].push_back(ppid);
794  }
795  }
796 }
797 
805 void PeriodicBoundaryHandler::processLocalGhostParticles(std::vector<BaseParticle*>& newParticles)
806 {
807  for (int i = 0; i < sendTargetList_.size(); i++)
808  {
809  //Check if the processor is sending ghosts to itself
811  {
812  for (int j = 0; j < numberOfNewPeriodicGhostParticlesSend_[i]; j++)
813  {
815  BaseParticle* particle = ppid->particle;
816 
817  //Create ghost particle
818  BaseParticle* pGhost = particle->copy();
819 
820  //Obtain and set the ghost periodic complexity
822 
823  //Shift the particle
824  shiftParticle(pGhost);
825 
826  //Set the correct flags
827  pGhost->setPeriodicGhostParticle(true);
828  pGhost->setInPeriodicDomain(true);
829  //Note: mpi flags dont have to be set for local particles
830  //these flags are copied from the original particle
831  logger(VERBOSE, "Adding a ghost with id % at position %", pGhost->getId(), pGhost->getPosition());
832 
833  //Finally add it to the particle handler
835 
836  //Do some bookkeeping
838  gpid->particle = pGhost;
839  gpid->otherParticle = particle;
840  gpid->realPeriodicComplexity = particle->getPeriodicComplexity();
841  periodicGhostList_[PROCESSOR_ID].push_back(gpid);
842 
843  //Flag real particle and do some bookkeeping
844  particle->setInPeriodicDomain(true);
845  ppid->otherParticle = pGhost;
846 
847  //Add to the new particle list for an interaction update
848  newParticles.push_back(pGhost);
849  }
850  }
851  }
852 }
853 
861 {
862  MPIContainer& communicator = MPIContainer::Instance();
863 
864  //For all lists that contain ghost particles
865  //The variable dataIndex indicates which index in update<...>Receive_ the data is located
866  int dataIndex = -1;
867  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
868  {
869  unsigned long numberOfParticles = periodicGhostList_[i].size();
870  bool global;
871  if (numberOfParticles > 0)
872  {
873  //Check if the update is global or from the current processor
874  if (i != PROCESSOR_ID)
875  {
876  global = true;
877  dataIndex++;
878  }
879  else
880  {
881  global = false;
882  }
883 
884  //Update all particles
885  for (int p = 0; p < numberOfParticles; p++)
886  {
888  BaseParticle* pGhost = pgid->particle;
889 
890  //Depending on where the particle is located the data is stored in a different place
891  if (global)
892  {
893  //logger(VERBOSE,"i: %, numberOfParticles: %, dataIndex: %",i,numberOfParticles, dataIndex);
894  //Check if this particle really belongs to the data that is send
895  logger.assert(pGhost->getId() == updatePositionDataReceive_[dataIndex][p].id,
896  "Periodic particle lists are not in syc");
897 
898  //Add the real position and velocity of the particles
899  //Note: The before updating the position is the position of the ghost
900  //Note: The received position and velocity values are of the real particle
901  //Note: It will be shifted to the correct values after the complexity is computed
902  pGhost->setPreviousPosition(pGhost->getPosition());
903  pGhost->setPosition(updatePositionDataReceive_[dataIndex][p].position);
904  pGhost->setOrientation(updatePositionDataReceive_[dataIndex][p].orientation);
905  pGhost->setVelocity(updateVelocityDataReceive_[dataIndex][p].velocity);
906  pGhost->setAngularVelocity(updateVelocityDataReceive_[dataIndex][p].angularVelocity);
907 
908  }
909  else
910  {
911  //MpiPeriodicGhostParticleID * pgip = periodicGhostList_[i][p];
912  pGhost->setPreviousPosition(pGhost->getPosition());
913  pGhost->setPosition(pgid->otherParticle->getPosition());
914  pGhost->setOrientation(pgid->otherParticle->getOrientation());
915  pGhost->setVelocity(pgid->otherParticle->getVelocity());
917  }
918 
919  //Move current periodic complexity to previous periodic complexity
921 
922  //Compute the new realPeriodicComplexity
923  std::vector<int> previousRealPeriodicComplexity = periodicGhostList_[i][p]->realPeriodicComplexity;
924  //The ghost particle has the real position at the moment, so compute the real periodic complexity here
925  std::vector<int> realPeriodicComplexity = computePeriodicComplexity(pGhost->getPosition());
926  std::vector<int> periodicComplexity(getSize());
927  for (int b = 0; b < getSize(); b++)
928  {
929  int sign = mathsFunc::sign(previousRealPeriodicComplexity[b]
930  * realPeriodicComplexity[b]
931  * pGhost->getPreviousPeriodicComplexity()[b]);
932  periodicComplexity[b] = sign * abs(realPeriodicComplexity[b]);
933  //The maser boundary needs this correction
934  if (periodicComplexity[b] == -3)
935  {
936  periodicComplexity[b] = 1;
937  }
938 
939  }
940  pGhost->setPeriodicComplexity(periodicComplexity);
941 
942  for (int b = 0; b < getSize(); b++)
943  {
944  objects_[b]->modifyGhostAfterCreation(pGhost, b);
945  }
946  //Shift the particle to correct position
947  //Note: If the real particle changed complexity this position will be calculated incorrectly
948  //Hence the previous complexity is used.
949  shiftParticle(pGhost, pGhost->getPreviousPeriodicComplexity());
950 
951  //Update hGrid
952  Vec3D displacement = pGhost->getPreviousPosition() - pGhost->getPosition();
953  getDPMBase()->hGridUpdateMove(pGhost, displacement.getLengthSquared());
954 
955  //Do some book keeping
956  periodicGhostList_[i][p]->realPeriodicComplexity = realPeriodicComplexity;
957 
958  }
959  }
960  }
961 
962  //Update periodic complexity of periodic particles
963  for (BaseParticle* particle : getDPMBase()->particleHandler)
964  {
965  if ((particle->isInPeriodicDomain()) && !(particle->isPeriodicGhostParticle()))
966  {
967  //Update the periodicComplexity of the real particles
968  particle->setPreviousPeriodicComplexity(particle->getPeriodicComplexity());
969  std::vector<int> periodicComplexity;
970  int totalPeriodicComplexity;
971  computePeriodicComplexity(periodicComplexity, totalPeriodicComplexity, particle->getPosition());
972  //Modify periodic complexity tailored to specific boundary requirements
973  for (int b = 0; b < getSize(); b++)
974  {
975  objects_[b]->modifyPeriodicComplexity(periodicComplexity, totalPeriodicComplexity, particle, b);
976  }
977  particle->setPeriodicComplexity(periodicComplexity);
978  }
979  }
980 }
981 
982 
990 bool PeriodicBoundaryHandler::checkIsReal(const std::vector<int> complexity)
991 {
992  bool isReal = true;
993  for (int i = 0; i < getSize(); i++)
994  {
995  if (mathsFunc::sign(complexity[i]) == -1)
996  {
997  isReal = false;
998  }
999  }
1000 
1001  return isReal;
1002 }
1003 
1014 bool PeriodicBoundaryHandler::checkChanged(const std::vector<int> previousComplexity,
1015  const std::vector<int> currentComplexity)
1016 {
1017  bool changed = false;
1018 
1019  for (int i = 0; i < currentComplexity.size(); i++)
1020  {
1021  if (previousComplexity[i] != currentComplexity[i])
1022  {
1023  changed = true;
1024  }
1025  }
1026 
1027  return changed;
1028 }
1029 
1038 void PeriodicBoundaryHandler::updateParticleStatus(std::set<BaseParticle*>& particlesToBeDeleted)
1039 {
1040  MPIContainer& communicator = MPIContainer::Instance();
1041  int processorID = communicator.getProcessorID();
1042  int numberOfProcessors = communicator.getNumberOfProcessors();
1043  std::set<MpiPeriodicParticleID*> deletePeriodicIDList;
1044  std::set<MpiPeriodicGhostParticleID*> deletePeriodicGhostIDList;
1045  std::set<BaseParticle*> specialUpdateParticleList;
1046 
1047  //For all domains
1048  for (int i = 0; i < numberOfProcessors; i++)
1049  {
1050  int numberOfPeriodicParticles = periodicParticleList_[i].size();
1051  int numberOfPeriodicGhostParticles = periodicGhostList_[i].size();
1052 
1053  //Loop over all periodic particles to see if their complexity changed
1054  for (int p = 0; p < numberOfPeriodicParticles; p++)
1055  {
1057  BaseParticle* particle = ppid->particle;
1058 
1059  //Check particle status
1060  bool isReal = checkIsReal(particle->getPeriodicComplexity());
1061  bool changed = checkChanged(particle->getPreviousPeriodicComplexity(), particle->getPeriodicComplexity());
1062 
1063  //Only if the particle changed we need to undertake action
1064  if (changed)
1065  {
1066  if (isReal)
1067  {
1068  //Flag this particle as normal, it will be re-introduced when finding new periodic particles
1069  logger(VERBOSE, "Real particle % changed complexity at: %", particle->getId(),
1070  particle->getPosition());
1071  particle->setInPeriodicDomain(false);
1072  //Incase of a special flag 3, perform update action
1073  updateMaserParticle(particle);
1074  }
1075  else
1076  {
1077  //Oh noes, the particle became a ghost. Kill it with balefire!!... if it is necessary
1078  logger(VERBOSE, "Real particle % changed to ghost at: %", particle->getId(),
1079  particle->getPosition());
1080  particlesToBeDeleted.insert(particle);
1081  }
1082 
1083  //Delete the ID
1084  deletePeriodicIDList.insert(ppid);
1085  periodicParticleList_[i][p] = nullptr;
1086  }
1087 
1088 
1089  //If a PM particle changes from to PMG it will need to be deleted.
1090  //The deletion will be done by the M boundary, but it still needs to be flushed from the periodic lists
1091  if (particle->isInMPIDomain())
1092  {
1093  //Store old values and compute new to see if the status has changed
1094  bool isMPIParticleOld = particle->isMPIParticle();
1095  bool isMPIParticle;
1096  bool isInMPIDomain;
1097  getMPIFlags(particle, isInMPIDomain, isMPIParticle);
1098 
1099  //Particle needs to be removed from lists if it becomes an MG particle
1100  if (isMPIParticleOld != isMPIParticle)
1101  {
1102  logger(VERBOSE, "PM to PMG: Flush from boundary");
1103  deletePeriodicIDList.insert(ppid);
1104  periodicParticleList_[i][p] = nullptr;
1105  }
1106 
1107  //MG particle needs to be removed from lists if it moves to another domain
1108  if (!isInMPIDomain)
1109  {
1110  if (particle->isMPIParticle())
1111  {
1112  logger(VERBOSE, "PMG leaves domain: Flush from boundary");
1113  deletePeriodicIDList.insert(ppid);
1114  periodicParticleList_[i][p] = nullptr;
1115  }
1116  }
1117  }
1118  }
1119 
1120  //Loop over all periodic ghost particles to see if their complexity changed
1121  for (int p = 0; p < numberOfPeriodicGhostParticles; p++)
1122  {
1124  BaseParticle* pGhost = pgid->particle;
1125 
1126  //Check particle status
1127  bool isReal = checkIsReal(pGhost->getPeriodicComplexity());
1128  bool changed = checkChanged(pGhost->getPreviousPeriodicComplexity(), pGhost->getPeriodicComplexity());
1129 
1130  //Update mixed particles, particles that also are in the mpi domain also need an update
1131  //Note that these particles are not listed in the domain lists, they are taken care here.
1132  //Store old values and update status of pGhost
1133  bool isInMPIDomainOld = pGhost->isInMPIDomain();
1134  bool isMPIParticleOld = pGhost->isMPIParticle();
1135  setMPIFlags(pGhost);
1136  if (isInMPIDomainOld)
1137  {
1138 
1139  //Case 1: pGhost changed from real to mpi particle
1140  if (isMPIParticleOld != pGhost->isMPIParticle())
1141  {
1142  //Case 1: turned from M ghost to M
1143  //The correct flags have been set above already
1144 
1145  //Case 2: Turned from M to M ghost
1146  if (pGhost->isMPIParticle())
1147  {
1148  logger(VERBOSE, "PGM to PGMG: Deleting particle.");
1149  particlesToBeDeleted.insert(pGhost);
1150  deletePeriodicGhostIDList.insert(pgid);
1151  periodicGhostList_[i][p] = nullptr;
1152  }
1153 
1154  }
1155 
1156  //Case 2: pGhost left the mpi domain
1157  if (pGhost->isInMPIDomain() != isInMPIDomainOld)
1158  {
1159  //Case 1: Moved inside the current domain
1160  //The correct flags have been set above already
1161 
1162  //Case 2: Moved to a neighbour domain
1163  if (pGhost->isMPIParticle())
1164  {
1165  logger(VERBOSE, "PGMG moved out of domain: Deleting particle.");
1166  particlesToBeDeleted.insert(pGhost);
1167  deletePeriodicGhostIDList.insert(pgid);
1168  periodicGhostList_[i][p] = nullptr;
1169  }
1170  }
1171  }
1172 
1173  //Check if the particles need to be deleted based on their periodic complexity
1174  if (changed)
1175  {
1176  if (isReal)
1177  {
1178  //Check if the complexity of the particle is truely real based on it's current position
1179  std::vector<int> pc = computePeriodicComplexity(pGhost->getPosition());
1180  int tpc = 0;
1181  for (int b = 0; b < getSize(); b++)
1182  {
1183  objects_[b]->modifyPeriodicComplexity(pc, tpc, pGhost, b);
1184  }
1185  if (!checkIsReal(pc))
1186  {
1187  logger(ERROR, "Round-off error detected.");
1188  //logger(WARN,"Round-off error corrected, phew!");
1189  }
1190 
1191  //There are two cases
1192  //Case 1: PGMG particles turning PMG; These are still not real
1193  //Case 2: PGM particles turning PG; These are truely real
1194  if (pGhost->isMPIParticle())
1195  {
1196  logger(VERBOSE, "PGMG turned PMG: delete it");
1197  particlesToBeDeleted.insert(pGhost);
1198  }
1199  else
1200  {
1201  //Turn the particle real
1202  logger(VERBOSE, "Ghost particle changed to real at position: %", pGhost->getPosition());
1203  pGhost->setInPeriodicDomain(false);
1204  pGhost->setPeriodicGhostParticle(false);
1205  //Make sure this particle can be detected now by the parallel boundaries
1206  pGhost->setInMPIDomain(false);
1207  }
1208  }
1209  else
1210  {
1211  //Delete the particle
1212  logger(VERBOSE, "Ghost particle changed complexity at position: %", pGhost->getPosition());
1213  particlesToBeDeleted.insert(pGhost);
1214  }
1215 
1216  //Delete the ID
1217  deletePeriodicGhostIDList.insert(pgid);
1218  periodicGhostList_[i][p] = nullptr;
1219  }
1220  }
1221 
1222  }
1223 
1224  //Delete IDs
1225  for (auto ppid_it : deletePeriodicIDList)
1226  {
1227  delete ppid_it;
1228  }
1229 
1230  for (auto pgid_it : deletePeriodicGhostIDList)
1231  {
1232  delete pgid_it;
1233  }
1234  unsigned int nextId = getDPMBase()->particleHandler.getNextId();
1235  communicator.broadcast(nextId);
1237 }
1238 
1253 int PeriodicBoundaryHandler::findTargetProcessor(const std::vector<int>& complexity)
1254 {
1255  //Find the middle of this domain (or any valid point anyway)
1256  Vec3D middlePosition = getDPMBase()->domainHandler.getCurrentDomain()->getMiddle();
1257 
1258  //Create the particle with a target position
1259  SphericalParticle particle;
1260  particle.setPosition(middlePosition);
1261  shiftParticle(&particle, complexity);
1262 
1263  //Obtain target domain
1264  int targetGlobalIndex = getDPMBase()->domainHandler.getParticleDomainGlobalIndex(&particle);
1265  return getDPMBase()->domainHandler.getParticleProcessor(targetGlobalIndex);
1266 }
1267 
1275 {
1276  //Skip the ghost particles
1277  if (checkIfAddNewParticle(particle))
1278  {
1279  //If the particle was not yet flagged to be in the periodic domain, check if it is in the domain
1280  //if(!particle->isInPeriodicDomain())
1281  {
1282  int totalPeriodicComplexity;
1283  std::vector<int> periodicComplexity;
1284  computePeriodicComplexity(periodicComplexity, totalPeriodicComplexity, particle->getPosition());
1285 
1286  //Modify periodic complexity tailored to specific boundary requirements
1287  for (int b = 0; b < getSize(); b++)
1288  {
1289  objects_[b]->modifyPeriodicComplexity(periodicComplexity, totalPeriodicComplexity, particle, b);
1290  }
1291 
1292  //Set periodicComplexity
1293  particle->setPeriodicComplexity(periodicComplexity);
1294 
1295  //Create ghost particle ID's
1296  std::vector<std::vector<int> > list(0);
1297  if (totalPeriodicComplexity > 0)
1298  {
1299  periodicComplexity = particle->getPeriodicComplexity();
1300 
1301  //Generating all possible complexities.
1302  generateGhosts(list, periodicComplexity, periodicComplexity, getSize());
1303  //logger(VERBOSE,"New ghost particles: %",list.size() - 1);
1304 
1305  //Note: the first one in the list is the real particle such that we skip it
1306  for (int i = 1; i < list.size(); i++)
1307  {
1308  //Compute target domain
1309  //logger(VERBOSE,"Particle position: %",particle->getPosition());
1310  int targetProcessor = findTargetProcessor(list[i]);
1311 
1312  //Create periodic particle ID
1314  ppid->particle = particle;
1315  ppid->periodicComplexity = periodicComplexity;
1316  ppid->targetPeriodicComplexity = list[i];
1317  ppid->targetProcessor = targetProcessor;
1318  newPeriodicParticleList_[targetProcessor].push_back(ppid);
1319  logger(VERBOSE, "Adding a periodic particle with id % and real particle position: %",
1320  particle->getId(), particle->getPosition());
1321  }
1322  }
1323  }
1324  }
1325 }
1326 
1332 {
1333  //Check for all particles if there are unassigned periodic particles
1334  for (auto particle_it = getDPMBase()->particleHandler.begin();
1335  particle_it != getDPMBase()->particleHandler.end(); particle_it++)
1336  {
1337  findNewParticle(*particle_it);
1338  }
1339 
1340 }
1341 
1347 {
1348  //Check all new particle lists
1349  newInteractionList_.resize(sendTargetList_.size());
1350  for (int i = 0; i < sendTargetList_.size(); i++)
1351  {
1352  for (int p = 0; p < newPeriodicParticleList_[sendTargetList_[i]].size(); p++)
1353  {
1354  BaseParticle* particle = newPeriodicParticleList_[sendTargetList_[i]][p]->particle;
1355 
1356  //Loop over all its interactions
1357  std::vector<BaseInteraction*> interactions = particle->getInteractions();
1358  for (BaseInteraction* interaction : interactions)
1359  {
1360  //Find out what the new particle is interacting with
1361  BaseParticle* particleP = dynamic_cast<BaseParticle*>(interaction->getP());
1362  BaseParticle* objectI = dynamic_cast<BaseParticle*>(interaction->getI());
1363 
1364  //If the P in the interaction structure is the new particle, find I
1365  if (particle == particleP)
1366  {
1367  //Check if the new particle is interacting with a wall
1368  if (!objectI)
1369  {
1370  newInteractionList_[i].push_back(interaction);
1371  }
1372  else //is I a particle
1373  {
1374  newInteractionList_[i].push_back(interaction);
1375  }
1376  }
1377  else //newBoundaryParticle is I in the interaction, P can only be a particle
1378  {
1379  newInteractionList_[i].push_back(interaction);
1380  }
1381  }
1382  }
1383 
1384  //Store the number of interactions
1386  }
1387 }
1388 
1398 {
1399 #ifdef MERCURY_USE_MPI
1400  //Step 1: Check to how many domains this domain has to send particles
1401  int numberOfTargetDomainsLocal = 0;
1402  for (int index = 0; index < NUMBER_OF_PROCESSORS; index++)
1403  {
1404  if (newPeriodicParticleList_[index].size() > 0)
1405  {
1406  numberOfTargetDomainsLocal++;
1407  sendTargetList_.push_back(index);
1409  }
1410  }
1411 
1412  //Communicate with other domains to how many domains this domain is sending
1413  std::vector<int> numberOfTargetDomains(NUMBER_OF_PROCESSORS);
1414  MPIContainer::Instance().allGather(numberOfTargetDomainsLocal,1,numberOfTargetDomains,1);
1415 
1416  //Step 2: Communicate to other domains what the target domains are
1417  std::vector<int> receiveTargetList;
1418  for (int index = 0; index < NUMBER_OF_PROCESSORS; index++)
1419  {
1420  if (numberOfTargetDomains[index] > 0)
1421  {
1422  if (index == PROCESSOR_ID)
1423  {
1424  MPIContainer::Instance().broadcast(sendTargetList_.data(), numberOfTargetDomainsLocal, index);
1425  }
1426  else
1427  {
1428  receiveTargetList.resize(numberOfTargetDomains[index], 0);
1429  MPIContainer::Instance().broadcast(receiveTargetList.data(),numberOfTargetDomains[index],index);
1430  }
1431 
1432  //If this domain is a target, flag it.
1433  for (int i = 0; i < receiveTargetList.size(); i++)
1434  {
1435  if (receiveTargetList[i] == PROCESSOR_ID)
1436  {
1437  receiveTargetList_.push_back(index);
1438  }
1439  }
1440  }
1441  }
1442 #endif
1443 }
1444 
1452 {
1453  MPIContainer& communicator = MPIContainer::Instance();
1454  //Communicate to how many particles the target domains receive
1455  int tagSend;
1456  for (int i = 0; i < sendTargetList_.size(); i++)
1457  {
1458  if (sendTargetList_[i] != PROCESSOR_ID)
1459  {
1461  communicator.send(numberOfNewPeriodicGhostParticlesSend_[i], sendTargetList_[i], tagSend);
1462 
1464  communicator.send(numberOfNewInteractionsSend_[i], sendTargetList_[i], tagSend);
1465  }
1466  }
1467 
1468  //Perform the receive routines
1471  int tagReceive;
1472  for (int i = 0; i < receiveTargetList_.size(); i++)
1473  {
1475  {
1478 
1480  communicator.receive(numberOfNewInteractionsReceive_[i], receiveTargetList_[i], tagReceive);
1481  }
1482  }
1483 }
1484 
1494 {
1495  //Communicate which domains send and which domains receive particles
1497 
1498  //Find new interactions
1500 
1501  //Communicate the number of particles and interactions send to the other domains
1503 
1504  //Synchronise the communications
1506 }
1507 
1509 {
1510  MPIContainer& communicator = MPIContainer::Instance();
1511  int numberOfTargetsReceive = receiveTargetList_.size();
1512  int numberOfTargetsSend = sendTargetList_.size();
1513 
1514  //Make sure that the receiving vectors have the correct length
1515  periodicGhostParticleReceive_.resize(numberOfTargetsReceive);
1516  periodicGhostComplexityReceive_.resize(numberOfTargetsReceive);
1517  interactionDataReceive_.resize(numberOfTargetsReceive);
1519  for (int i = 0; i < numberOfTargetsReceive; i++)
1520  {
1524  }
1525 
1526  //Collect data for sending
1529 
1530  //Send particle and interaction data
1531  for (int i = 0; i < numberOfTargetsSend; i++)
1532  {
1533  if (sendTargetList_[i] != PROCESSOR_ID)
1534  {
1535  //Send particle data
1536  int sendCount = numberOfNewPeriodicGhostParticlesSend_[i];
1537  int processor = sendTargetList_[i];
1539  communicator.send(periodicGhostParticleSend_[i].data(), MercuryMPIType::PARTICLE, sendCount, processor,
1540  tagSend);
1541 
1542  //Send complexity
1545  communicator.send(periodicGhostComplexitySend_[i].data(), sendCount, processor, tagSend);
1546 
1547  //Send interactions
1549  sendCount = numberOfNewInteractionsSend_[i];
1550  if (sendCount > 0)
1551  {
1552  communicator.send(interactionDataSend_[i], MercuryMPIType::INTERACTION, sendCount, processor, tagSend);
1553  }
1554  }
1555  }
1556 
1557  //Receive data
1558  for (int i = 0; i < numberOfTargetsReceive; i++)
1559  {
1561  {
1562  //Receive particle data
1563  int receiveCount = numberOfNewPeriodicGhostParticlesReceive_[i];
1564  int processor = receiveTargetList_[i];
1567  receiveCount, processor, tagReceive);
1568 
1569  //Receive complexity
1572  communicator.receive(periodicGhostComplexityReceive_[i].data(), receiveCount, processor, tagReceive);
1573 
1574  //Send interactions
1576  receiveCount = numberOfNewInteractionsReceive_[i];
1577  if (receiveCount > 0)
1578  {
1579  communicator.receive(interactionDataReceive_[i], MercuryMPIType::INTERACTION, receiveCount, processor,
1580  tagReceive);
1581  }
1582  }
1583  }
1584 
1585  //Synchronise the communications
1586  communicator.sync();
1587 }
1588 
1590 {
1591  //Process the global received data
1592  for (int i = 0; i < receiveTargetList_.size(); i++)
1593  {
1594  std::vector<BaseParticle*> newGhostParticles;
1595  processReceivedGhostParticleData(i, newGhostParticles);
1596  processReceivedInteractionData(i, newGhostParticles);
1597  }
1598 
1599  //Process the local data
1600  std::vector<BaseParticle*> newGhostParticles;
1601  processLocalGhostParticles(newGhostParticles);
1602  processLocalInteractionData(newGhostParticles);
1603 
1604  //Process the periodic particles;
1606 
1607  //Clear lists
1608  receiveTargetList_.clear();
1609  sendTargetList_.clear();
1618  interactionDataSend_.clear();
1619  interactionDataReceive_.clear();
1620  newInteractionList_.clear();
1621  for (auto& i : newPeriodicParticleList_)
1622  {
1623  i.clear();
1624  }
1625  for (auto& i : newPeriodicParticleList_)
1626  {
1627  i.clear();
1628  }
1629 }
1630 
1632 {
1633  MPIContainer& communicator = MPIContainer::Instance();
1634  unsigned int numberOfProcessors = communicator.getNumberOfProcessors();
1635  unsigned int processorID = communicator.getProcessorID();
1636 
1637  //For all lists that contain periodic particles
1638  for (unsigned int i = 0; i < numberOfProcessors; i++)
1639  {
1640  unsigned int numberOfParticles = periodicParticleList_[i].size();
1641  if (numberOfParticles > 0 && i != processorID)
1642  {
1643  //Increase the vector size;
1644  updatePositionDataSend_.emplace_back(0);
1645  updateVelocityDataSend_.emplace_back(0);
1646 
1647  //Collect the data
1648  for (unsigned int p = 0; p < numberOfParticles; p++)
1649  {
1650  BaseParticle* particle = periodicParticleList_[i][p]->particle;
1651  updatePositionDataSend_.back().push_back(copyPositionFrom(particle));
1652  updateVelocityDataSend_.back().push_back(copyVelocityFrom(particle));
1653  }
1654 
1655  //Send position data
1656  unsigned int count = numberOfParticles;
1657  unsigned int processor = i;
1658  unsigned int tag = processorID * MAX_PROC + processor * 10 + MercuryMPITag::POSITION_DATA;
1659  communicator.send(updatePositionDataSend_.back().data(), MercuryMPIType::POSITION, count, processor, tag);
1660 
1661  //Send velocity data
1662  tag = processorID * MAX_PROC + processor * 10 + MercuryMPITag::VELOCITY_DATA;
1663  communicator.send(updateVelocityDataSend_.back().data(), MercuryMPIType::VELOCITY, count, processor, tag);
1664  }
1665  }
1666 
1667  //For all lists that contain ghost particles
1668  for (int i = 0; i < numberOfProcessors; i++)
1669  {
1670  unsigned int numberOfParticles = periodicGhostList_[i].size();
1671  if (numberOfParticles > 0 && i != processorID)
1672  {
1673  //Increase the vector size
1674  updatePositionDataReceive_.emplace_back(numberOfParticles);
1675  updateVelocityDataReceive_.emplace_back(numberOfParticles);
1676 
1677  //Receive position data
1678  int count = numberOfParticles;
1679  int processor = i;
1680  int tag = processor * MAX_PROC + processorID * 10 + MercuryMPITag::POSITION_DATA;
1681  communicator.receive(updatePositionDataReceive_.back().data(), MercuryMPIType::POSITION, count, processor,
1682  tag);
1683 
1684  //Receive velocity data
1685  tag = processor * MAX_PROC + processorID * 10 + MercuryMPITag::VELOCITY_DATA;
1686  communicator.receive(updateVelocityDataReceive_.back().data(), MercuryMPIType::POSITION, count, processor,
1687  tag);
1688  }
1689  }
1690 
1691  //Synchronise all the requests
1692  communicator.sync();
1693 }
1694 
1696 {
1697  //Update the positions, velocities and periodic complexities of particles
1698  updateParticles();
1699 
1700  //Delete vectors
1701  updatePositionDataSend_.clear();
1703  updateVelocityDataSend_.clear();
1705 }
1706 
1707 //TODO this is for the deletion boundary?
1708 void PeriodicBoundaryHandler::flushParticles(std::set<BaseParticle*>& particlesToBeFlushed)
1709 {
1714  std::set<MpiPeriodicParticleIDBase*> toBeDeleted;
1715  for (auto p_it : particlesToBeFlushed)
1716  {
1717  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
1718  {
1719  for (auto& p : periodicParticleList_[i])
1720  {
1721  if (p != nullptr)
1722  {
1723  if (p_it == p->particle)
1724  {
1725  toBeDeleted.insert(p);
1726  p = nullptr;
1727  }
1728  }
1729  }
1730  }
1731 
1732  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
1733  {
1734  for (auto& p : periodicGhostList_[i])
1735  {
1736  if (p != nullptr)
1737  {
1738  if (p_it == p->particle)
1739  {
1740  toBeDeleted.insert(p);
1741  p = nullptr;
1742  }
1743  }
1744  }
1745  }
1746  }
1747 
1748  //Delete ID's
1749  for (auto id_it : toBeDeleted)
1750  {
1751  delete id_it;
1752  }
1753 
1754 }
1755 
1762 void PeriodicBoundaryHandler::cleanCommunicationList(std::vector<MpiPeriodicParticleIDBase*>& list)
1763 {
1764  for (int i = 0; i < list.size(); i++)
1765  {
1766  if (list[i] == nullptr)
1767  {
1768  list[i] = list.back();
1769  list.pop_back();
1770  i--;
1771  }
1772  }
1773 }
1774 
1776 {
1777  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
1778  {
1781  }
1782 }
1783 
1789 //TODO optimise the target finding strategy such that only "NEW target/receives" are added
1791 {
1792  newPeriodicParticleList_ = std::vector<PeriodicList>(NUMBER_OF_PROCESSORS, PeriodicList(0));
1793  periodicParticleList_ = std::vector<PeriodicList>(NUMBER_OF_PROCESSORS, PeriodicList(0));
1794  periodicGhostList_ = std::vector<PeriodicGhostList>(NUMBER_OF_PROCESSORS, PeriodicGhostList(0));
1795 }
1796 
1798 {
1799  for (BasePeriodicBoundary* boundary : *this)
1800  {
1802  }
1803 }
1804 
1811 {
1812  if (NUMBER_OF_PROCESSORS > 1)
1813  {
1814  //Clear ID lists
1815  for (auto& i : periodicParticleList_)
1816  {
1817  //logger(INFO,"Size: %",periodicParticleList_[i].size());
1818  for (int j = 0; j < i.size(); j++)
1819  {
1820  delete i[j];
1821  }
1822  i.clear();
1823  }
1824 
1825  for (auto& i : periodicGhostList_)
1826  {
1827  for (int j = 0; j < i.size(); j++)
1828  {
1829  delete i[j];
1830  }
1831  i.clear();
1832  }
1833 
1834  //Collect all ghost particles and unflag
1836  std::set<BaseParticle*> toBeDeleted; //Set because I am too lazy to implement a vector based flushParticle function
1837  for (BaseParticle* particle : pH)
1838  {
1839  if (particle->isPeriodicGhostParticle())
1840  {
1841  toBeDeleted.insert(particle);
1842  }
1843 
1844  //Unflag its periodicity
1845  particle->setInPeriodicDomain(false);
1846  }
1847 
1848  //Flush from mpi boundary and clean the lists
1849  getDPMBase()->getCurrentDomain()->flushParticles(toBeDeleted);
1851 
1852  //Delete particles
1853  int index = 0;
1854  for (BaseParticle* particle : toBeDeleted)
1855  {
1856  pH.removeGhostObject(particle->getIndex());
1857  index++;
1858  }
1859  }
1860 }
1861 
1870 {
1871  if (particle->isMaserParticle())
1872  {
1873  for (int b = 0; b < getSize(); b++)
1874  {
1875  if (particle->getPeriodicComplexity(b) == 3)
1876  {
1877  particle->setMaserParticle(false);
1878 
1879  logger(VERBOSE, "particle % with position % goes into outflow domain, new ID = %", particle->getId(),
1880  particle->getPosition(), getDPMBase()->particleHandler.getNextId());
1881  const unsigned int newID = getDPMBase()->particleHandler.getNextId();
1882  logger(VERBOSE, "new id % position X %", newID, particle->getPosition().X);
1883  particle->setId(newID);
1885  }
1886  }
1887  }
1888 }
1889 
1891 {
1892  if (particle->isPeriodicGhostParticle())
1893  {
1894  return false;
1895  }
1896  if (particle->isInPeriodicDomain())
1897  {
1898  return false;
1899  }
1900  return true;
1901 }
1902 
1903 
1904 
1905 
1906 
1907 
1908 
1909 
1910 
1911 
1912 
std::vector< int > numberOfNewInteractionsSend_
Stores the number of new interactions to be send to target processor corresponding to sendTargetList_...
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
std::vector< int > realPeriodicComplexity
Definition: MpiDataClass.h:161
void performNewParticleTransmission()
Collects and sends the ghost particle data.
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
void addNewParticle(BaseParticle *particle)
Adds a new particle to the periodic list.
unsigned int getNextId()
Definition: BaseHandler.h:256
bool isInInnerDomain(BaseParticle *particle)
Check if the particle is in the current domain but not in the communication zone. ...
Definition: Domain.cc:430
std::vector< MpiPeriodicGhostParticleID * > PeriodicGhostList
std::vector< std::vector< int > > periodicGhostComplexityReceive_
Data container for periodic complexity that is being received from other processors.
bool isInPeriodicDomain() const
Indicates if the particle is in the periodic boundary communication zone.
void finalisePositionAndVelocityUpdate()
Communicates position and velocity data from periodic boundaries and updates ghost particles...
bool checkChanged(std::vector< int > previousComplexity, std::vector< int > complexity)
checks of two periodic complexities differ
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
Mdouble X
the vector components
Definition: Vector.h:65
Domain * getCurrentDomain()
Function that returns a pointer to the domain corresponding to the processor.
Definition: DPMBase.cc:5034
std::vector< MpiPeriodicParticleID * > PeriodicList
A basic particle.
void setPreviousPeriodicComplexity(std::vector< int > complexity)
Set the previous periodic communication complexity of the paritcle.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
void increaseId()
Definition: BaseHandler.h:251
virtual void hGridGetInteractingParticleList(BaseParticle *obj, std::vector< BaseParticle * > &list)
Creates a list of neighbour particles obtained from the hgrid.
Definition: DPMBase.h:936
bool isMaserParticle() const
Indicates if this particle belongs to the maser boundary.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
void setInMPIDomain(bool flag)
Flags the status of the particle if wether it is in the communication zone or not.
void updateParticleStatus(std::set< BaseParticle * > &particlesToBeDeleted)
Updates the status of periodic particles and ghost particles.
virtual Mdouble getDistance(const BaseParticle &particle) const =0
Returns the distance between a particle and the closest boundary, required for any periodic boundary...
void communicateTargetDomains()
Creats a list of send and receive targets for periodic/ghost particles.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void setDPMBase(DPMBase *DPMBase)
Sets the problem that is solved using this handler.
virtual void shiftPosition(BaseParticle *particle) const =0
Shifts the position (and velocity) of to the ghost particle.
void processLocalGhostParticles(std::vector< BaseParticle * > &newParticles)
Creates ghost particles of periodic particles that are located on the same processor.
void finaliseNewParticleTransmission()
creates the ghost particles and performs some bookkeeping to keep track of them
void * createMPIInteractionDataArray(unsigned int numberOfInteractions) const
creates an empty MPIInteractionDataArray
const std::vector< int > & getPreviousPeriodicComplexity() const
Sets the previous periodic communication complexity of the particle.
BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler) override
Checks if particle is in interaction with given particle P, and if so, returns vector of pointer to t...
T * getObjectById(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:565
std::vector< std::vector< int > > periodicGhostComplexitySend_
Data container for periodic complexity that is being send to other processors.
std::vector< std::vector< MPIParticleVelocity > > updateVelocityDataSend_
Data container for velocity data that is being send to other processors.
void updateMaserParticle(BaseParticle *particle)
Updates the maser flag of particles leaving the maser.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void flushParticles(std::set< BaseParticle * > &toBeDeletedList)
Particles that are going to be deleted from the simulation are flushed out of the communication bound...
Definition: Domain.cc:1669
std::enable_if< std::is_scalar< T >::value, void >::type receive(T &t, int from, int tag)
asynchronously receive a scalar from some other processor.
Definition: MpiContainer.h:213
virtual BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler)=0
Returns the interaction between this object and a given BaseParticle.
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:704
void addGhostObject(int fromPrcessor, int toProcessor, BaseParticle *p)
Adds a ghost particle located at fromProcessor to toProcessor.
bool containsParticle(BaseParticle *particle, Mdouble offset=0.0)
Check to see if a given particle is within the current domain.
Definition: Domain.cc:400
MpiPeriodicParticleIDBase MpiPeriodicParticleID
Definition: MpiDataClass.h:166
std::size_t getNumberOfProcessors() const
Get the total number of processors participating in this simulation.
Mdouble interactionDistance_
The interaction distance between a position and the boundary for which particles start to participate...
DomainHandler domainHandler
An object of the class DomainHandler which deals with parallel code.
Definition: DPMBase.h:1354
std::vector< PeriodicList > newPeriodicParticleList_
Mdouble getNumberOfTruePeriodicGhostParticles()
Returns the number of particles that are flagged as periodicGhostParticles, but not as MPIParticles...
void flushParticles(std::set< BaseParticle * > &particlesToBeFlushed)
Removes particles from the periodiocParticleList_ and periociGhostList_.
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
bool isInGreaterDomain(BaseParticle *particle)
Check to see if a given particle is in the current domain or in neighbouring communication zones...
Definition: Domain.cc:420
void findNewParticle(BaseParticle *particle)
Checks if a particle is in the periodic domain, but are not flagged as being in the periodic domain...
bool checkIsReal(std::vector< int > complexity)
checks if a periodic complexity is real
int sign(T val)
This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0...
Definition: ExtendedMath.h:95
std::vector< int > numberOfNewInteractionsReceive_
Stores the number of new interactions to be received from target processor corresponding to receiveTa...
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorfism.
int findTargetProcessor(const std::vector< int > &complexity)
For a given complexity this function returns the target processor.
void findNewParticles()
Loops over all particles in the simulation to check if they need to be added to the periodic lists...
void setPeriodicHandler(PeriodicBoundaryHandler *periodicHandler)
Sets the periodicBoundaryHandler, required for parallel periodic boundaries.
BaseInteraction * getInteractionWith(BaseParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler) override
Returns the interaction between this wall and a given particle, nullptr if there is no interaction...
Definition: BaseWall.cc:369
std::vector< std::vector< MPIParticle > > periodicGhostParticleSend_
Data container for particles that are being send to other processors.
std::vector< std::vector< MPIParticlePosition > > updatePositionDataReceive_
Data container for position data that is being received from other processors.
Stores information about interactions between two interactable objects; often particles but could be ...
void collectGhostParticleData()
Collects ghost particle data that needs to be be sent to other processors.
void performActionsBeforeAddingParticles()
Actions that boundaries perform before adding new periodic/ghost particles.
void initialise()
Initialises the communication list vectors as they can not be determined on compile time...
std::vector< void * > interactionDataReceive_
Stores the interaction data that is going to be received.
void generateGhosts(std::vector< std::vector< int > > &list, std::vector< int > periodicComplexity, std::vector< int > &complexity, int level)
generates a list of periodic complexities corresponding to a give real particle.
void setMPIFlags(BaseParticle *particle)
Sets the MPIParticle and isMPIParticle flags of a given particle.
std::vector< std::vector< BaseInteraction * > > newInteractionList_
Container to store pointers to all BasePeriodicBoundary objects.
std::vector< int > numberOfNewPeriodicGhostParticlesReceive_
A vector that stores how many new ghost particles will be received from other processors.
void copyDataFromMPIParticleToParticle(MPIParticle *bP, BaseParticle *p, ParticleHandler *particleHandler)
Copies data from an MPIParticle class to a BaseParticle and sets the particleHandler and species...
~PeriodicBoundaryHandler() override
Destructor, it destructs the PeriodicBoundaryHandler and all BasePeriodicBoundary it contains...
std::size_t getProcessorID()
Reduces a scalar on all processors to one scalar on a target processor.
void getInteractionDetails(void *interactionData, unsigned int index, unsigned int &identificationP, unsigned int &identificationI, bool &isWallInteraction, unsigned &timeStamp)
reads the basic interaction details from an MPIInteractionDataArray
void processLocalInteractionData(std::vector< BaseParticle * > &newParticles)
Process the interaction data for local ghosts.
void communicateNumberOfNewParticlesAndInteractions()
Communicate the number of new particles and interactions to target processors.
void findNewInteractions()
Finds interactions that accompany future ghost particles.
#define MAX_PROC
Definition: GeneralDefine.h:51
void processPeriodicParticles()
Creates a periodioc particle ID for book keeping and moves the ID to the correct list.
static BaseParticle * newParticle()
MPIParticle copyDataFromParticleToMPIParticle(BaseParticle *p)
Copies data from a SuperQuadricParticle to an MPIParticle class and returns this. ...
void cleanCommunicationLists()
Removes nullptrs from boundaryParticleList_ and boundaryParticleListNeighbour_.
Definition: Domain.cc:1713
bool isInMPIDomain()
Indicates if the particle is in the communication zone of the mpi domain.
PeriodicBoundaryHandler operator=(const PeriodicBoundaryHandler &rhs)
Assignment operator, copies only the vector of BasePeriodicBoundary and sets the other variables to 0...
void updateStatus(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
Updates the positions/velocity of ghost particles and accordingly the status of these particles...
void setPreviousPosition(const Vec3D &pos)
Sets the particle's position in the previous time step.
void addNewParticles()
Adds new particles to the periodic particle lists.
unsigned int getNumberOfPeriodicGhostParticles()
Returns the number of particles that are flagged is periodicGhostParticle.
void shiftParticle(BaseParticle *particle)
Shifts the position of the particle based on its current periodic complexity.
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
std::string getName() const override
Returns the name of the handler, namely the string "PeriodicBoundaryHandler".
Container to store Interaction objects.
void collectInteractionData()
Collects interaction data into an MPI data structure.
std::vector< int > targetPeriodicComplexity
Definition: MpiDataClass.h:160
std::vector< PeriodicGhostList > periodicGhostList_
A vector the size of the number of processors, each entry containing a vector of ghost periodioc part...
#define PROCESSOR_ID
Definition: GeneralDefine.h:63
std::vector< PeriodicList > periodicParticleList_
A vector the size of the number of processors, each entry containing a vector of periodic particle ID...
std::vector< std::vector< MPIParticlePosition > > updatePositionDataSend_
Data container for position data that is being send to other processors.
std::vector< int > periodicComplexity
Definition: MpiDataClass.h:159
std::vector< BasePeriodicBoundary * > objects_
The actual list of Object pointers.
Definition: BaseHandler.h:302
void getMPIFlags(BaseParticle *particle, bool &isInMPIDomain, bool &isMPIParticle)
Determines if a given particle is in the MPI domain and if it is an MPI Particle. ...
Basic class for walls.
Definition: BaseWall.h:47
void processReceivedInteractionData(int targetIndex, std::vector< BaseParticle * > &newParticles)
Process the received interaction data.
bool checkIfAddNewParticle(BaseParticle *particle)
std::vector< int > numberOfNewPeriodicGhostParticlesSend_
A vector that stores how many particles are going to be send to other processors. ...
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1359
void setNextId(unsigned int id)
Definition: BaseHandler.h:261
bool isPeriodicGhostParticle() const
Indicates if this particle is a ghost in the periodic boundary.
Container to store all BaseParticle.
#define NUMBER_OF_PROCESSORS
For the MPI communication routines this quantity is often required. defining this macro makes the cod...
Definition: GeneralDefine.h:62
std::vector< void * > interactionDataSend_
Stores the interaction data that is going to be send.
MPIParticlePosition copyPositionFrom(BaseParticle *particle)
Copies the position from a particle to an MPIParticlePosition class.
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1339
std::enable_if< std::is_scalar< T >::value, void >::type broadcast(T &t, int fromProcessor=0)
Broadcasts a scalar from the root to all other processors.
Definition: MpiContainer.h:429
int getParticleDomainGlobalIndex(BaseParticle *particle)
void setMaserParticle(bool flag)
Flags the status of the particle if it belongs to the maser boundary or not.
const std::vector< int > & getPeriodicComplexity()
Obtains the periodic communication complexity of the particle.
The simulation can be subdivided into Domain's used in parallel code.
Definition: Domain.h:63
PeriodicBoundaryHandler()
Default constructor, it creates an empty PeriodicBoundaryHandler.
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void sync()
Process all pending asynchronous communication requests before continuing.
Definition: MpiContainer.h:148
void prepareNewParticleTransmission()
Initial preparation work for sending ghost particles.
std::enable_if< std::is_scalar< T >::value, void >::type send(T &t, int to, int tag)
Asynchronously send a scalar to some other processor.
Definition: MpiContainer.h:167
Vec3D getMiddle() const
Gives the middle of the domain.
Definition: Domain.cc:1704
void updateParticles()
Updates position/velocity and periodic complexity of ghost particles.
Domain * getCurrentDomain()
Gets the domain assigned to the processor.
int getParticleProcessor(int globalIndex)
void setPeriodicComplexity(std::vector< int > complexity)
Set the periodic communication complexity of the particle.
Mdouble getInteractionDistance()
Returns the interaction distance.
const Vec3D & getPreviousPosition() const
Returns the particle's position in the previous time step.
Definition: BaseParticle.h:400
bool isMPIParticle() const
Indicates if this particle is a ghost in the MPI domain.
void setMPIParticle(bool flag)
Flags the mpi particle status.
const std::vector< BaseInteraction * > & getInteractions() const
Returns a list of interactions which belong to this interactable.
void setPeriodicGhostParticle(bool flag)
Flags the status of the particle to be a ghost in periodic boundary or not.
void readAndAddObject(std::istream &is) override
Pure virtual function needs implementation, but it does nothing for the periodicBoudnaryHandler.
std::vector< int > computePeriodicComplexity(Vec3D position)
Computes the periodic complexity based on a given position.
std::vector< std::vector< MPIParticleVelocity > > updateVelocityDataReceive_
Data container for velocity data that is being received from other processors.
Defines the basic properties that a interactable object can have.
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
virtual void setMPIInteraction(void *interactionDataArray, unsigned int index, bool resetPointers)
virtual void hGridUpdateMove(BaseParticle *, Mdouble)
Definition: DPMBase.cc:1893
T * getLastObject()
Gets a pointer to the last Object in this BaseHandler.
Definition: BaseHandler.h:634
virtual void performActionsBeforeAddingParticles()
Actions that need to be performed before adding new ghost particles.
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Definition: Vector.h:49
void addObject(BasePeriodicBoundary *P) override
Adds a BasePeriodicBoundary to the PeriodicBoundaryHandler.
BaseParticle * otherParticle
Definition: MpiDataClass.h:151
MPIParticleVelocity copyVelocityFrom(BaseParticle *particle)
Copies the velocity from a particle to an MPIParticleVelocity class.
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
void preparePositionAndVelocityUpdate()
Collects the position and velocity data from periodic boundaries.
void cleanCommunicationList(std::vector< MpiPeriodicParticleIDBase * > &list)
Removes the nullptr's from a communication list.
void setInPeriodicDomain(bool flag)
Flags the status of the particle wether it is in the periodic communication zone or not...
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void processReceivedGhostParticleData(int targetIndex, std::vector< BaseParticle * > &newParticles)
Processes the received ghost data, creates a ghost particle and does some book keeping.
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
std::vector< std::vector< MPIParticle > > periodicGhostParticleReceive_
Data container for particles that are being received from other processors.
MpiPeriodicParticleIDBase MpiPeriodicGhostParticleID
Definition: MpiDataClass.h:167
void clearCommunicationLists()
Removes all ghost particles and bookkeeping for a fresh start.
virtual bool isClosestToLeftBoundary(const Vec3D &position) const =0
Returns true if it is closer to the left boundary than the right boundary.
void setId(unsigned long id)
Assigns a unique identifier to each object in the handler (container) which remains constant even aft...
Definition: BaseObject.cc:72
std::vector< int > receiveTargetList_
A list that keeps track which target processors the current processor is receiving new particles from...
std::vector< int > sendTargetList_
A list that keeps track to which targets this processor is sending new particles to.
void setInteractionDistance(Mdouble interactionDistance)
Sets the interaction distance.