From c5901b0fd5db799e3e365a43657f29dc650f2089 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Fri, 21 Feb 2020 18:12:47 +0100
Subject: [PATCH] [BUGFIX] checking of simulation is working again

---
 apps/benchmarks/GranularGas/GranularGas.cfg   |  2 ++
 .../GranularGas/MESA_PD_GranularGas.cpp       | 22 +++++++++----------
 .../GranularGas/MESA_PD_KernelBenchmark.cpp   | 19 ++++++++--------
 .../MESA_PD_KernelLoadBalancing.cpp           |  2 +-
 .../GranularGas/MESA_PD_LoadBalancing.cpp     |  2 +-
 apps/benchmarks/GranularGas/check.h           |  4 ++--
 6 files changed, 26 insertions(+), 25 deletions(-)

diff --git a/apps/benchmarks/GranularGas/GranularGas.cfg b/apps/benchmarks/GranularGas/GranularGas.cfg
index 69be192fa..647eca2ef 100644
--- a/apps/benchmarks/GranularGas/GranularGas.cfg
+++ b/apps/benchmarks/GranularGas/GranularGas.cfg
@@ -8,6 +8,8 @@ GranularGas
    sorting linear;
    shift < 0.01, 0.01, 0.01 >;
 
+   checkSimulation 0;
+
    LBAlgorithm Morton;
    baseWeight 1;
 
diff --git a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
index cb19bf5cd..59682aa03 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
@@ -74,13 +74,12 @@ namespace mesa_pd {
 class DEM
 {
 public:
-   DEM(const std::shared_ptr<domain::BlockForestDomain>& domain)
+   DEM(const std::shared_ptr<domain::BlockForestDomain>& domain, real_t dt, real_t mass)
    : domain_(domain)
    {
-      dem_.setStiffness(0, 0, real_t(0));
-      dem_.setDampingN(0, 0, real_t(0));
       dem_.setDampingT(0, 0, real_t(0));
       dem_.setFriction(0, 0, real_t(0));
+      dem_.setParametersFromCOR(0, 0, real_t(0.9), dt*real_t(20), mass * real_t(0.5));
    }
 
    inline
@@ -208,20 +207,20 @@ int main( int argc, char ** argv )
 
    if (!forest->isPeriodic(0))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(+1,0,0));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(-1,0,0));
+      createPlane(*ps, *ss, confiningDomain.minCorner()+params.shift, Vec3(+1,0,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+params.shift, Vec3(-1,0,0));
    }
 
    if (!forest->isPeriodic(1))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,+1,0));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,-1,0));
+      createPlane(*ps, *ss, confiningDomain.minCorner()+params.shift, Vec3(0,+1,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+params.shift, Vec3(0,-1,0));
    }
 
    if (!forest->isPeriodic(2))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,0,+1));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,0,-1));
+      createPlane(*ps, *ss, confiningDomain.minCorner()+params.shift, Vec3(0,0,+1));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+params.shift, Vec3(0,0,-1));
    }
 
    WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
@@ -237,7 +236,7 @@ int main( int argc, char ** argv )
    WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
    // Init kernels
    kernel::ExplicitEulerWithShape        explicitEulerWithShape( params.dt );
-   DEM dem(domain);
+   DEM dem(domain, params.dt, ss->shapes[smallSphere]->getMass());
    kernel::InsertParticleIntoLinkedCells ipilc;
    kernel::AssocToBlock                  assoc(forest);
    mpi::ReduceProperty                   RP;
@@ -356,7 +355,8 @@ int main( int argc, char ** argv )
 
       if (params.checkSimulation)
       {
-         check(*ps, *forest, params.spacing);
+         //if you want to activate checking you have to deactivate sorting
+         check(*ps, *forest, params.spacing, params.shift);
       }
 
       WALBERLA_LOG_INFO_ON_ROOT("*** SQL OUTPUT - START ***");
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
index 0ae1b0e67..d2c1452c4 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
@@ -92,7 +92,6 @@ public:
    inline
    void operator()(const size_t idx1, const size_t idx2, T& ac)
    {
-
       ++contactsChecked_;
       if (double_cast_(idx1, idx2, ac, acd_, ac))
       {
@@ -233,20 +232,20 @@ int main( int argc, char ** argv )
 
    if (!forest->isPeriodic(0))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(+1,0,0));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(-1,0,0));
+      createPlane(*ps, *ss, confiningDomain.minCorner() + params.shift, Vec3(+1,0,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner() + params.shift, Vec3(-1,0,0));
    }
 
    if (!forest->isPeriodic(1))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,+1,0));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,-1,0));
+      createPlane(*ps, *ss, confiningDomain.minCorner() + params.shift, Vec3(0,+1,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner() + params.shift, Vec3(0,-1,0));
    }
 
    if (!forest->isPeriodic(2))
    {
-      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,0,+1));
-      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,0,-1));
+      createPlane(*ps, *ss, confiningDomain.minCorner() + params.shift, Vec3(0,0,+1));
+      createPlane(*ps, *ss, confiningDomain.maxCorner() + params.shift, Vec3(0,0,-1));
    }
 
    WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
@@ -266,10 +265,9 @@ int main( int argc, char ** argv )
    kernel::AssocToBlock                  assoc(forest);
    kernel::InsertParticleIntoLinkedCells ipilc;
    kernel::SpringDashpot                 dem(1);
-   dem.setStiffness(0, 0, real_t(0));
-   dem.setDampingN (0, 0, real_t(0));
    dem.setDampingT (0, 0, real_t(0));
    dem.setFriction (0, 0, real_t(0));
+   dem.setParametersFromCOR(0, 0, real_t(0.9), params.dt*real_t(20), ss->shapes[smallSphere]->getMass() * real_t(0.5));
    mpi::ReduceProperty                   RP;
    mpi::SyncNextNeighbors                SNN;
    ContactDetection                      CD(domain);
@@ -376,7 +374,8 @@ int main( int argc, char ** argv )
 
       if (params.checkSimulation)
       {
-         check(*ps, *forest, params.spacing);
+         //if you want to activate checking you have to deactivate sorting
+         check(*ps, *forest, params.spacing, params.shift);
       }
 
       WALBERLA_LOG_INFO_ON_ROOT("*** SQL OUTPUT - START ***");
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
index 44375e3f5..dfff27199 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
@@ -495,7 +495,7 @@ int main( int argc, char ** argv )
 
    if (params.checkSimulation)
    {
-      check(*ps, *forest, params.spacing);
+      check(*ps, *forest, params.spacing, params.shift);
    }
 
 
diff --git a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
index 6048d5b7b..5c658a301 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
@@ -483,7 +483,7 @@ int main( int argc, char ** argv )
 
    if (params.checkSimulation)
    {
-      check(*ps, *forest, params.spacing);
+      check(*ps, *forest, params.spacing, params.shift);
    }
 
    WALBERLA_LOG_INFO_ON_ROOT("*** SQL OUTPUT - START ***");
diff --git a/apps/benchmarks/GranularGas/check.h b/apps/benchmarks/GranularGas/check.h
index 31099e229..73c457592 100644
--- a/apps/benchmarks/GranularGas/check.h
+++ b/apps/benchmarks/GranularGas/check.h
@@ -27,13 +27,13 @@
 namespace walberla {
 namespace mesa_pd {
 
-void check( data::ParticleStorage& ps, blockforest::BlockForest& forest, real_t spacing )
+void check( data::ParticleStorage& ps, blockforest::BlockForest& forest, real_t spacing, const Vec3& shift )
 {
    WALBERLA_LOG_INFO_ON_ROOT("*** CHECKING RESULT - START ***");
    auto pIt = ps.begin();
    for (auto& iBlk : forest)
    {
-      for (auto it = grid_generator::SCIterator(iBlk.getAABB(), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing);
+      for (auto it = grid_generator::SCIterator(iBlk.getAABB(), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5) + shift, spacing);
            it != grid_generator::SCIterator();
            ++it, ++pIt)
       {
-- 
GitLab