Initialize at openlb-1-3

1 files changed, 2257 insertions, 0 deletions

diff --git a/src/particles/superParticleSystem3D.hh b/src/particles/superParticleSystem3D.hh
new file mode 100644
index 0000000..f6450f7
--- /dev/null
+++ b/src/particles/superParticleSystem3D.hh
@@ -0,0 +1,2257 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2016 Thomas Henn
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+ */
+
+#ifndef SUPERPARTICLESYSTEM_3D_HH
+#define SUPERPARTICLESYSTEM_3D_HH
+
+#define shadows
+
+#include <utility>
+#include <string>
+#include <array>
+#include <iomanip>
+#include <ios>
+#include <iostream>
+#include <random>
+#include "io/fileName.h"
+#include "superParticleSystem3D.h"
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+namespace olb {
+
+template<typename T, template<typename U> class PARTICLETYPE>
+SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
+  CuboidGeometry3D<T>& cuboidGeometry, LoadBalancer<T>& loadBalancer,
+  SuperGeometry3D<T>& superGeometry)
+  : SuperStructure3D<T>(cuboidGeometry, loadBalancer,
+                        superGeometry.getOverlap()),
+    clout(std::cout, "SuperParticleSystem3d"),
+    _superGeometry(superGeometry),
+    _overlap(0)
+{
+  init();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
+  SuperGeometry3D<T>& superGeometry)
+  : SuperStructure3D<T>(superGeometry.getCuboidGeometry(),
+                        superGeometry.getLoadBalancer(),
+                        superGeometry.getOverlap()),
+    clout(std::cout, "SuperParticleSystem3d"),
+    _superGeometry(superGeometry),
+    _overlap(0)
+{
+  init();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::init()
+{
+  int rank = 0;
+  int size = 0;
+
+  _stopSorting = 0;
+
+#ifdef PARALLEL_MODE_MPI
+  rank = singleton::mpi().getRank();
+  size = singleton::mpi().getSize();
+#endif
+  for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) {
+    //clout << i << std::endl;
+    if (this->_loadBalancer.rank(i) == rank) {
+      auto dummy = new ParticleSystem3D<T, PARTICLETYPE>(i, _superGeometry);
+      this->_cuboidGeometry.get(i).getOrigin().toStdVector()[0];
+      std::vector<T> physPos = this->_cuboidGeometry.get(i).getOrigin().toStdVector();
+      std::vector<T> physExtend(3, 0);
+      T physR = this->_cuboidGeometry.get(i).getDeltaR();
+      for (int j = 0; j < 3; j++) {
+        physPos[j] -= .5 * physR;
+        physExtend[j] = (this->_cuboidGeometry.get(i).getExtend()[j] + 1)
+                        * physR;
+      }
+      dummy->setPosExt(physPos, physExtend);
+      _pSystems.push_back(dummy);
+    }
+  }
+  //singleton::exit(0);
+
+#ifdef PARALLEL_MODE_MPI
+  for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) {
+    if (this->_loadBalancer.rank(i) == rank) {
+      std::vector<int> dummy;
+      this->getCuboidGeometry().getNeighbourhood(i, dummy, 3);
+      _rankNeighbours.insert(_rankNeighbours.end(), dummy.begin(), dummy.end());
+      _cuboidNeighbours.push_back(dummy);
+    }
+  }
+  for (auto& N : _rankNeighbours) {
+    N = this->_loadBalancer.rank(N);
+  }
+#endif
+
+  /* Ein jeder ist sein eigener Nachbar*/
+  if (rank == 0) {
+    _rankNeighbours.push_back(size - 1);
+  }
+  if (rank == size - 1) {
+    _rankNeighbours.push_back(0);
+  }
+  _rankNeighbours.push_back(rank);
+  _rankNeighbours.sort();
+  _rankNeighbours.unique();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
+  SuperParticleSystem3D<T, PARTICLETYPE>& spSys)
+  : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer,
+                        int(spSys._overlap)),
+    clout(std::cout, "SuperParticleSystem3d"),
+    _pSystems(spSys._pSystems),
+    _superGeometry(spSys._superGeometry),
+    _rankNeighbours(spSys._rankNeighbours),
+    _overlap(spSys._overlap)
+{
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
+  SuperParticleSystem3D<T, PARTICLETYPE> const& spSys)
+  : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer,
+                        int(spSys._overlap)),
+    clout(std::cout, "SuperParticleSystem3d"),
+    _pSystems(spSys._pSystems),
+    _superGeometry(spSys._superGeometry),
+    _rankNeighbours(spSys._rankNeighbours),
+    _overlap(spSys._overlap)
+{
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
+  SuperParticleSystem3D<T, PARTICLETYPE> && spSys)
+  : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer,
+                        int(spSys._overlap)),
+    clout(std::cout, "SuperParticleSystem3d"),
+    _superGeometry(spSys._superGeometry),
+    _rankNeighbours(spSys._rankNeighbours),
+    _overlap(spSys._overlap)
+{
+  _pSystems = std::move(spSys._pSystems);
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::print()
+{
+  int no = globalNumOfParticles();
+  int active = globalNumOfActiveParticles();
+  clout << "activeParticles= " << active << " (" << no << ") " << std::endl;
+  //  cout << "[SuperParticleSystem3D] " << _pSystems.size()
+  //      << " pSystems on rank " << singleton::mpi().getRank() << "\n";
+  //  for (auto pS : _pSystems) {
+  //    cout << pS->_particles.size() << " ";
+  //  }
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::printDeep(std::string message)
+{
+  clout << "=========================================================================================================================" << std::endl;
+  clout << "printDeep diagnostic tool" << message << std::endl;
+  clout << "=========================================================================================================================" << std::endl;
+  for (auto pS : _pSystems) {
+    pS->printDeep ( std::to_string(_pSystems.size()) + " pSystems on rank " + std::to_string(singleton::mpi().getRank()) + ":  " );
+  }
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::print(std::list<int> mat)
+{
+  std::list<int>::iterator _matIter;
+  int no = globalNumOfParticles();
+  clout << "globalNumOfParticles=" << no;
+  int active = globalNumOfActiveParticles();
+  clout << "; activeParticles=" << active;
+  for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
+    clout << "; material" << *_matIter << "=" << countMaterial(*_matIter);
+  }
+  clout << std::endl;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::captureEscapeRate(
+  std::list<int> mat)
+{
+  std::list<int>::iterator _matIter;
+  T sum = T();
+  // capture rate
+  for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
+    sum += (T) countMaterial(*_matIter);
+  }
+  clout << "captureRate=" << 1. - sum / globalNumOfParticles()
+        << "; escapeRate=" << sum / globalNumOfParticles() << std::endl;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate(
+  std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate,
+  int iT, int iTConsole, T genPartPerTimeStep)
+{
+  std::list<int>::iterator _matIter;
+  T pSumOutletNew = T();
+  T maxDiffEscapeRateNew = maxDiffEscapeRate;
+  T diffERtmp = T(); // temporal differencial escape rate
+
+  // count particles at outlet
+  for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
+    pSumOutletNew += (T) countMaterial(*_matIter);
+  }
+
+  // calculate diff. escape rate
+  if (globalNumOfParticles() > globalPSum) {
+    T diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum);
+    diffEscapeRate += diffERtmp;
+  } else {
+    if (genPartPerTimeStep != 0.) {
+      diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep);
+      diffEscapeRate += diffERtmp;
+    }
+  }
+  // calculate max. diff. escape rate
+  if (diffERtmp > maxDiffEscapeRateNew) {
+    maxDiffEscapeRateNew = diffERtmp;
+  }
+  // console output
+  if (iT % iTConsole == 0) {
+    diffEscapeRate /= iTConsole;
+    if (globalNumOfParticles() > globalPSum) {
+      clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
+    } else {
+      if (genPartPerTimeStep != 0.) {
+        clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical "
+              << genPartPerTimeStep  << " generated particles per phys. time step"
+              << " diffEscapeRate = " << diffEscapeRate << std::endl;
+      } else {
+        clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl;
+      }
+    }
+    if (maxDiffEscapeRateNew > maxDiffEscapeRate) {
+      clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl;
+    }
+    diffEscapeRate = 0. ;
+  }
+  maxDiffEscapeRate = maxDiffEscapeRateNew;
+  pSumOutlet = pSumOutletNew;
+  globalPSum = globalNumOfParticles();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate(
+  std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate,
+  int iT, int iTConsole, T genPartPerTimeStep,
+  T& avDiffEscapeRate, T latticeTimeStart, T latticeTimeEnd)
+{
+  std::list<int>::iterator _matIter;
+  T pSumOutletNew = T();
+  T maxDiffEscapeRateNew = maxDiffEscapeRate;
+  T avDiffEscapeRateNew = T();
+
+  // count particle at outlet
+  for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
+    pSumOutletNew += (T) countMaterial(*_matIter);
+  }
+  // calculate diff. escape rate
+  if (globalNumOfParticles() > globalPSum) {
+    avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum);
+    diffEscapeRate += avDiffEscapeRateNew;
+  } else {
+    if (genPartPerTimeStep != 0.) {
+      avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep);
+      diffEscapeRate += avDiffEscapeRateNew;
+    }
+  }
+  // calculate max. diff. escape rate
+  if (avDiffEscapeRateNew > maxDiffEscapeRateNew) {
+    maxDiffEscapeRateNew = avDiffEscapeRateNew;
+  }
+  // calculate average diff. escape rate between tStart and tEnd
+  if ((iT >= latticeTimeStart) && (iT <= latticeTimeEnd)) {
+    avDiffEscapeRate += avDiffEscapeRateNew;
+  }
+  if (iT == latticeTimeEnd) {
+    avDiffEscapeRate /= (latticeTimeEnd - latticeTimeStart);
+    clout << "average diffEscapeRate between t = " << latticeTimeStart << "s and t = " << latticeTimeEnd << "s : "
+          << avDiffEscapeRate << std::endl;
+  }
+  // console output
+  if (iT % iTConsole == 0) {
+    diffEscapeRate /= iTConsole;
+    if (globalNumOfParticles() > globalPSum) {
+      clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
+    } else {
+      if (genPartPerTimeStep != 0.) {
+        clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical "
+              << genPartPerTimeStep  << " generated particles per phys. time step" << std::endl;
+        clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
+      } else {
+        clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl;
+      }
+    }
+    if (maxDiffEscapeRateNew > maxDiffEscapeRate) {
+      clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl;
+    }
+    diffEscapeRate = 0. ;
+  }
+  maxDiffEscapeRate = maxDiffEscapeRateNew;
+  pSumOutlet = pSumOutletNew;
+  globalPSum = globalNumOfParticles();
+}
+
+// Output_txt.file for tracer particles (particles that have an id != 0
+/* _filename is filename of output txt-file
+ * _file is contact to output file
+ *  */
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::getOutput(std::string _filename,
+    int it, T conversionFactorTime, unsigned short _properties )
+{
+
+  enum particleProperties
+    : unsigned short {position = 1,
+                    velocity = 2,
+                    radius = 4,
+                    mass = 8,
+                    force = 16,
+                    storeForce = 32,
+                   };
+
+  /// initialize Parameters for the desired ParticleOutput:
+  /// ID, Radius, Position, Velocity, Force, storeForces
+  std::vector < T > id, rad;
+  std::vector<std::array<T, 3>> pos, vel, forces, storeForces;
+  /// initialize iteration Parameters for the loops
+  int k, j, numOfTracerParticles = globalNumOfTracerParticles();
+  int m = 0;
+  /// setPrecision for the decimal places in the txtOutputFile
+  std::setprecision(9);
+  /// set size of list in correlation to the number of TracerParticles
+  for (k = 0; k < numOfTracerParticles; k++) {
+    id.push_back(T());
+    rad.push_back(T());
+    pos.push_back( { T(), T(), T() });
+    vel.push_back( { T(), T(), T() });
+    forces.push_back( { T(), T(), T() });
+    storeForces.push_back( { T(), T(), T() });
+  }
+
+  k = 0;
+  /// get Information about each Particle in each particleSystem
+  std::vector<ParticleSystem3D<T, PARTICLETYPE>*> _pSystems = getPSystems();
+  for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
+    std::deque<PARTICLETYPE<T>*> particles =
+      _pSystems[pS]->getParticlesPointer();
+    for (auto p : particles) {
+      if (p->getID() != 0) {
+        /// check whether it is a tracerParticle(!=0) or not (==0)
+        id[k] = p->getID();
+        rad[k] = p->getRad();
+        for (j = 0; j < 3; j++) { /// loop for each coordinate [x,y,z]
+          pos[k][j] = p->getPos()[j];
+          vel[k][j] = p->getVel()[j];
+          forces[k][j] = p->getForce()[j];
+          storeForces[k][j] = p->getStoreForce()[j];
+        }
+       // p->resetStoreForce();
+        k++;
+        p++;
+      }
+    }
+  }
+
+#ifdef PARALLEL_MODE_MPI
+  /// prepare data vectors for parallel computing (mpi)
+  for (m = 0; m < numOfTracerParticles; m++) {
+    singleton::mpi().reduceAndBcast(id[m], MPI_SUM);
+    singleton::mpi().reduceAndBcast(rad[m], MPI_SUM);
+    for (j = 0; j < 3; j++) { /// loop for each coordinate [x,y,z]
+      singleton::mpi().reduceAndBcast(pos[m][j], MPI_SUM);
+      singleton::mpi().reduceAndBcast(vel[m][j], MPI_SUM);
+      singleton::mpi().reduceAndBcast(forces[m][j], MPI_SUM);
+      singleton::mpi().reduceAndBcast(storeForces[m][j], MPI_SUM);
+    }
+  }
+#endif
+
+  /// Particle Output to txtFile just on main process ID (==0)
+  if (!singleton::mpi().getRank()) {
+    std::ofstream _file;
+
+    int i = 0;
+    if (it == 0) {
+      /// write headers of each column at timeStep zero
+      /// after header name, there is the number of the column
+      _file.open(_filename, std::ios::out | std::ios::trunc);
+      _file << "Timestep" << i + 1 << " " << "physTime" << i + 2;
+      i = i + 2;
+
+      for (m = 0; m < numOfTracerParticles; m++) {
+        _file << " id" << i + 1;
+        i = i + 1;
+        _file << " rad" << i + 1;
+        i = i + 1;
+        if (_properties & particleProperties::position) {
+          _file << " pos0_" << i + 1 << " pos1_" << " pos2_" << i + 3;
+          i = i + 3;
+        }
+        if (_properties & particleProperties::force) {
+          _file << " forc0_" << i + 1 << " forc1_" << i + 2 << " forc2_"
+                << i + 3;
+          i = i + 3;
+        }
+        if (_properties & particleProperties::velocity) {
+          _file << " vel0_" << i + 1 << " vel1_" << i + 2 << " vel2_" << i + 3;
+          i = i + 3;
+        }
+        if (_properties & particleProperties::storeForce) {
+          _file << " hforc0_" << i + 1 << " hforc1_" << i + 2 << " hforc2_"
+                << i + 3;
+          i = i + 3;
+        }
+      }
+      _file << "\n";
+      _file.close();
+    }
+
+    if (it >= 0) {
+      /// write the results in regart to the above defined headers
+      _file.open(_filename, std::ios::out | std::ios::app);
+      _file << it << " " << conversionFactorTime*it << " ";
+
+      for (m = 0; m < numOfTracerParticles; m++) {
+        _file << id[m] << " ";
+        _file << rad[m] << " ";
+        if (_properties & particleProperties::position) {
+          _file << pos[m][0] << " " << pos[m][1] << " " << pos[m][2] << " ";
+        }
+        if (_properties & particleProperties::force) {
+          _file << forces[m][0] << " " << forces[m][1] << " " << forces[m][2]
+                << " ";
+        }
+        if (_properties & particleProperties::velocity) {
+          _file << vel[m][0] << " " << vel[m][1] << " " << vel[m][2] << " ";
+        }
+        if (_properties & particleProperties::storeForce) {
+          _file << storeForces[m][0] << " " << storeForces[m][1] << " "
+                << storeForces[m][2] << " ";
+        }
+      }
+    }
+    _file << "\n";
+    _file.close();
+  }
+}
+
+
+template<typename T, template<typename U> class PARTICLETYPE>
+std::vector<ParticleSystem3D<T, PARTICLETYPE>*>& SuperParticleSystem3D<T,
+    PARTICLETYPE>::getPSystems()    //+*
+{
+  return _pSystems;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::simulate(T dT, bool scale)
+{
+
+  //  for (auto pS : _pSystems) {
+  //      pS->velocityVerlet1(dT);
+  //      if (_overlap > 0) {
+  //        pS->makeTree();
+  //        cout << "makeTree" << endl;
+  //      }
+  //  }
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //      pS->velocityVerlet2(dT);
+  //  }
+
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //    pS->computeForce();
+  //    pS->predictorCorrector1(dT);
+  //  }
+  //
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //    pS->computeForce();
+  //    pS->predictorCorrector2(dT);
+  //  }
+
+  //  for (auto pS : _pSystems) {
+  //    pS->rungeKutta4_1(dT);
+  //  }
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //    pS->rungeKutta4_2(dT);
+  //  }
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //    pS->rungeKutta4_3(dT);
+  //  }
+  //  updateParticleDistribution();
+  //  for (auto pS : _pSystems) {
+  //    pS->rungeKutta4_4(dT);
+  //  }
+  //  updateParticleDistribution();
+  for (auto pS : _pSystems) {
+    time_t delta = clock();
+    pS->_contactDetection->sort();
+    _stopSorting += clock() - delta;
+    pS->simulate(dT, scale);
+    pS->computeBoundary();
+
+  }
+  updateParticleDistribution();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Mathias ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps, bool resetExternalField, bool scale )
+{
+  // reset external field
+  if (resetExternalField)
+    backCoupling.resetExternalField(material);
+
+  for (auto pS : _pSystems) {
+    time_t delta = clock();
+    pS->_contactDetection->sort();
+    _stopSorting += clock() - delta;
+    pS->simulateWithTwoWayCoupling_Mathias(dT, forwardCoupling, backCoupling, material, subSteps, scale);
+    pS->computeBoundary();
+
+  }
+  updateParticleDistribution();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Davide ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps, bool resetExternalField, bool scale )
+{
+  // reset external field
+  if (resetExternalField)
+    backCoupling.resetExternalField(material);
+
+  for (auto pS : _pSystems) {
+    time_t delta = clock();
+    pS->_contactDetection->sort();
+    _stopSorting += clock() - delta;
+    pS->simulateWithTwoWayCoupling_Davide(dT, forwardCoupling, backCoupling, material, subSteps, scale);
+    pS->computeBoundary();
+  }
+  updateParticleDistribution();
+}
+
+// multiple collision models
+template<>
+void SuperParticleSystem3D<double, MagneticParticle3D>::simulate(double dT, std::set<int> sActivityOfParticle, bool scale)
+{
+  for (auto pS : _pSystems) {
+    time_t delta = clock();
+    if (pS->getIGeometry() == singleton::mpi().getRank()) {
+      pS->_contactDetection->sort();
+    }
+    _stopSorting += clock() - delta;
+    if (pS->getIGeometry() == singleton::mpi().getRank()) {
+      pS->simulate(dT, sActivityOfParticle, scale);
+      pS->computeBoundary();
+    }
+  }
+  updateParticleDistribution();
+}
+
+template<>
+bool SuperParticleSystem3D<double, MagneticParticle3D>::particleSActivityTest(int sActivity)
+{
+  for (auto pS : _pSystems) {
+    for (auto p : pS->_particles) {
+      if (p.getSActivity() == sActivity) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::setOverlap(T overlap)
+{
+  if ( int(overlap) + 1 > _superGeometry.getOverlap() ) {
+    clout << "Overlap of SuperParticleSystem3D should be < overlap "
+          "of SuperStructure3D" << std::endl;
+    exit(-1);
+  }
+  _overlap = overlap;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+T SuperParticleSystem3D<T, PARTICLETYPE>::getOverlap()
+{
+  return _overlap;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::numOfPSystems()
+{
+  return _pSystems.size();
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfParticles()
+{
+#ifdef PARALLEL_MODE_MPI
+  // cout << "return1" << endl;
+  int buffer = rankNumOfParticles();
+  // cout << "return2" << endl;
+  singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
+  // cout << "return3" << endl;
+  return buffer;
+#else
+  // cout << "return4" << endl;
+  return rankNumOfParticles();
+#endif
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfShadowParticles()
+{
+#ifdef PARALLEL_MODE_MPI
+  int buffer = rankNumOfShadowParticles();
+  singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
+  return buffer;
+#else
+  return rankNumOfShadowParticles();
+#endif
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfParticles()
+{
+  int num = 0;
+  for (auto pS : _pSystems) {
+    num += pS->size();
+  }
+  return num;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfShadowParticles()
+{
+  int num = 0;
+  for (auto pS : _pSystems) {
+    num += pS->_shadowParticles.size();
+  }
+  return num;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfActiveParticles()
+{
+  int num = 0;
+  for (auto pS : _pSystems) {
+    num += pS->numOfActiveParticles();
+  }
+  return num;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::countLocMaterial(int mat)
+{
+  int num = 0;
+  for (auto pS : _pSystems) {
+    num += pS->countMaterial(mat);
+  }
+  return num;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::countMaterial(int mat)
+{
+#ifdef PARALLEL_MODE_MPI
+  int buffer = countLocMaterial(mat);
+  singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
+  return buffer;
+#else
+  return countLocMaterial(mat);
+#endif
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfActiveParticles()
+{
+#ifdef PARALLEL_MODE_MPI
+  int buffer = rankNumOfActiveParticles();
+  singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
+  return buffer;
+#else
+  return rankNumOfActiveParticles();
+#endif
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfTracerParticles()
+{
+#if PARALLEL_MODE_MPI
+  int buffer = rankNumOfTracerParticles();
+  singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
+  return buffer;
+#else
+  return rankNumOfTracerParticles();
+#endif
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfTracerParticles()
+{
+  int num = 0;
+  for (auto pS : _pSystems) {
+    std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer();
+    int pSNum = 0;
+    for (auto p : particles) {
+      if (p->getID() != 0) {
+        pSNum++;
+      }
+    }
+    num += pSNum;
+  }
+  return num;
+}
+
+// TODO class olb::ParticleSystem3D<double, olb::Particle3D>’ has no member named ‘radiusDistribution
+/*
+ template<typename T, template<typename U> class PARTICLETYPE>
+ std::map<T, int> SuperParticleSystem3D<T, PARTICLETYPE>::rankRadiusDistribution()
+ {
+ std::map<T, int> distrAll;
+ typename std::map<T, int>::iterator ita;
+ for (auto pS : _pSystems) {
+ std::map<T, int> distrPs;
+ distrPs = pS->radiusDistribution();
+ for (typename std::map<T, int>::iterator itp = distrPs.begin();
+ itp != distrPs.end(); ++itp) {
+ T radKeyP = itp->first;
+ int countP = itp->second;
+ if (distrAll.count(radKeyP) > 0) {
+ ita = distrAll.find(radKeyP);
+ ita->second = ita->second + countP;
+ } else {
+ distrAll.insert(std::pair<T, int>(radKeyP, countP));
+ }
+ }
+ }
+ return distrAll;
+ }
+ */
+
+template<typename T, template<typename U> class PARTICLETYPE>
+std::vector<int> SuperParticleSystem3D<T, PARTICLETYPE>::numOfForces()
+{
+  std::vector<int> dummy;
+  for (auto pS : _pSystems) {
+    dummy.push_back(pS->numOfForces());
+  }
+  return dummy;
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+template<typename DESCRIPTOR>
+void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToFluidVel(
+  SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR>& fVel)
+{
+  for (auto pS : _pSystems) {
+    pS->setVelToFluidVel(fVel);
+  }
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToAnalyticalVel(
+  AnalyticalConst3D<T, T>& aVel)
+{
+  for (auto pS : _pSystems) {
+    pS->setVelToAnalyticalVel(aVel);
+  }
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p)
+{
+  return findCuboid(p, int(_overlap));
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p,
+    int overlap)
+{
+  int C = this->_cuboidGeometry.get_iC(p.getPos()[0], p.getPos()[1],
+                                       p.getPos()[2], overlap);
+  if (C != this->_cuboidGeometry.getNc()) {
+    p.setCuboid(C);
+    return 1;
+  } else {
+    clout << "Lost Particle! Pos: " << p.getPos()[0] << " " << p.getPos()[1]
+          << " " << p.getPos()[2] << " Vel: " << p.getVel()[0] << " "
+          << p.getVel()[1] << " " << p.getVel()[2] << " Cuboid: " << C << std::endl;
+    p.setActive(false);
+    return 0;
+  }
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p,
+    T overlap)
+{
+  return checkCuboid(p, overlap, p.getCuboid());
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p,
+    T overlap, int iC)
+{
+  // Checks whether particle is contained in the cuboid
+  // extended with an layer of size overlap*delta
+  return this->_cuboidGeometry.get(iC).physCheckPoint(p.getPos()[0],
+         p.getPos()[1],
+         p.getPos()[2], overlap);
+}
+
+template<typename T, template<typename U> class PARTICLETYPE>
+void SuperParticleSystem3D<T, PARTICLETYPE>::updateParticleDistribution()
+{
+  // Find particles on wrong cuboid, store in relocate and delete
+  // maps particles to their new rank
+  // std::multimap<int, PARTICLETYPE<T> > relocate;
+  _relocate.clear();
+#ifdef shadows
+  _relocateShadow.clear();
+#endif
+  for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
+
+    _pSystems[pS]->_shadowParticles.clear();
+    auto par = _pSystems[pS]->_particles.begin();
+
+    while (par != _pSystems[pS]->_particles.end()) {
+
+      //Check if particle is still in its cuboid
+      if (checkCuboid(*par, 0)) {
+#ifdef shadows
+        // Check if its inside boundary layer of its cuboid
+        if (!(checkCuboid(*par, -_overlap))) {
+
+          std::set<int> sendTo;
+          // Run through all cuboids to search the one that
+          // shares its boundary with the cuboid containing the particle
+          for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) {
+            // Check if iC is neighbor cuboid in which overlap the particle is
+            // located
+            if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) {
+              // rank is the processing thread of the global cuboid number iC
+              int rank = this->_loadBalancer.rank(iC);
+              // insert rank into sendTo, if sendTo does not already contain
+              // rank for the actual particle
+              if (!sendTo.count(rank)) {
+                _relocateShadow.insert(std::make_pair(rank, (*par)));
+                sendTo.insert(rank);
+              }
+            }
+          }
+        }
+#endif
+        //If not --> find new cuboid
+      } else {
+        // check to which cuboid particle is located
+        // and give new cuboid to par
+        findCuboid(*par, 0);
+        // gives particle the new cuboid number where it is now located to
+        _relocate.insert(
+          std::make_pair(this->_loadBalancer.rank(par->getCuboid()), (*par)));
+#ifdef shadows
+        // If the new particle position is in boundary layer
+        // If yes --> check if its inside boundary layer
+        if (!(checkCuboid(*par, -_overlap))) {
+          // yes, particle is in boundary layer
+          std::set<int> sendTo;
+          for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) {
+            // checks if iC is neighbor cuboid in whose overlap the particle
+            // is located
+            if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) {
+              int rank = this->_loadBalancer.rank(iC);
+              if (!sendTo.count(rank)) {
+                // shadow particle (copy of original) is inserted for iC
+                // but just with information of rank, not of cuboid!
+                _relocateShadow.insert(std::make_pair(rank, (*par)));
+                sendTo.insert(rank);
+              }
+            }
+          }
+        }
+#endif
+        par = _pSystems[pS]->_particles.erase(par);
+        par--;
+      }
+      par++;
+    }
+  }
+  /* Communicate number of Particles per cuboid*/
+#ifdef PARALLEL_MODE_MPI
+  singleton::MpiNonBlockingHelper mpiNbHelper;
+  //mpiNbHelper.allocate(_rankNeighbours.size());
+
+  int k = 0;
+
+  /* Serialize particles */
+  // it is: std::map<int, std::vector<double> > _send_buffer
+  _send_buffer.clear();
+
+  // fill std::map<int, std::vector<double> > _send_buffer with
+  // particle properties of particles in _relocate
+
+  // buffer size equals number of particle properties
+  T buffer[PARTICLETYPE<T>::serialPartSize];
+  // insert all regular particles of _relocate into _send_buffer
+  for (auto rN : _relocate) {
+    // second element in rN is (*par)
+    // write particle properties into buffer
+    rN.second.serialize(buffer);
+    // first element in rN is rank
+    // it is: std::map<int, std::vector<double> > _send_buffer;
+    // enlarge the _send_buffer element (rank), which is of vector type, with
+    // new elements of buffer (begin: buffer, end: buffer + serialPartSize)
+    _send_buffer[rN.first].insert(_send_buffer[rN.first].end(),
+                                  buffer, buffer + PARTICLETYPE<T>::serialPartSize);
+  }
+
+  /*Send Particles */
+  k = 0;
+  // noSends contains number of neighboring cuboids
+  int noSends = 0;
+  // it is: std::list<int> _rankNeighbours, rank of neighboring cuboids
+  for (auto rN : _rankNeighbours) {
+    // if there are particles contained in _send_buffer, increase noSends
+    if (_send_buf