From 46d893d581621d65b17e70a9a08e15aeb62a3613 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Sat, 27 Jan 2018 15:11:48 +0100
Subject: [PATCH] Remove boost random, math::special_functions

Also remove some unneeded boost includes
---
 .../MeshDistance/MeshDistanceBenchmark.cpp    |  4 +-
 .../MotionSingleHeavySphere.cpp               |  2 +-
 src/core/debug/CheckFunctions.h               |  1 -
 src/core/math/DistributedSample.cpp           | 11 ++-
 src/core/math/FPClassify.h                    |  8 +-
 src/core/math/GenericAABB.h                   |  2 +-
 src/core/math/GenericAABB.impl.h              |  8 +-
 src/core/math/Random.cpp                      |  6 +-
 src/core/math/Random.h                        | 80 +++++++++++++++----
 src/core/math/Sample.cpp                      |  6 +-
 src/core/math/equation_system/Equation.cpp    |  7 +-
 src/core/math/equation_system/Variable.cpp    |  4 +-
 src/domain_decomposition/IBlock.h             |  1 -
 src/gather/CurveGatherPackInfo.impl.h         |  1 -
 tests/blockforest/BlockDataIOTest.cpp         |  2 +-
 .../blockforest/StructuredBlockForestTest.cpp |  4 +-
 tests/core/cell/CellIntervalTest.cpp          | 23 +++---
 tests/core/cell/CellTest.cpp                  |  8 +-
 tests/core/load_balancing/MetisTest.cpp       |  1 -
 tests/core/load_balancing/ParMetisTest.cpp    |  1 -
 tests/core/math/GenericAABBTest.cpp           | 12 +--
 tests/core/math/PlaneTest.cpp                 | 18 ++---
 tests/core/math/PrimesTest.cpp                |  8 +-
 tests/core/mpi/BufferSystemTest.cpp           | 11 +--
 tests/core/mpi/BufferTest.cpp                 | 39 +++++----
 .../GPUPackInfoCommunicationTest.cpp          |  2 +-
 tests/field/AccuracyEvaluationTest.cpp        |  4 +-
 tests/field/FieldFileIOTest.cpp               |  2 +-
 tests/field/FieldMPIDatatypesTest.cpp         | 10 +--
 tests/geometry/RandomPointGenerator.h         | 20 ++---
 tests/geometry/VoxelFileTest.cpp              | 67 +++++++++++++---
 tests/lbm/geometry/IntersectionRatioTest.cpp  |  4 +-
 tests/mesh/MeshAABBIntersectionTest.cpp       |  4 +-
 tests/mesh/MeshAABBSelectionTest.cpp          |  4 +-
 tests/mesh/MeshContainmentOctreeTest.cpp      |  4 +-
 tests/mesh/MeshDistanceCompareTest.cpp        |  4 +-
 tests/mesh/MeshDistanceOctreeTest.cpp         |  4 +-
 tests/mesh/MeshDistancePlausibilityTest.cpp   |  1 -
 tests/mesh/MeshVTKTest.cpp                    |  2 -
 tests/mesh/PeVTKMeshWriterTest.cpp            | 16 ++--
 tests/mesh/QHullTest.cpp                      | 18 +++--
 tests/pe/SyncEquivalence.cpp                  |  4 +-
 42 files changed, 268 insertions(+), 170 deletions(-)

diff --git a/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp b/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
index 64e9001fc..0e9585346 100644
--- a/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
+++ b/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
@@ -34,7 +34,7 @@
 #include "mesh/distance_octree/DistanceOctree.h"
 #include "mesh/MeshIO.h"
 
-#include <boost/random.hpp>
+#include <random>
 
 #include <algorithm>
 #include <vector>
@@ -64,7 +64,7 @@ void runBenchmark( const std::string & meshFile, const uint_t numPoints, const u
    timer.end();
    WALBERLA_LOG_INFO( "Mesh preparation took " << timer.last() << "s" );
 
-   boost::random::mt19937 rng;  
+   std::mt19937 rng;
 
    std::vector< typename MeshType::Point > points( numPoints );
    for(auto it = points.begin(); it != points.end(); ++it)
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index 10c97a06f..c54228f17 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -831,7 +831,7 @@ int main( int argc, char **argv )
       else if( int(Galileo) == 250 )
       {
          // add random perturbance for chaotic regime
-         walberla::math::seedRandomGenerator( boost::mt19937::result_type(std::time(0)) );
+         walberla::math::seedRandomGenerator( std::mt19937::result_type(std::time(0)) );
          xParticle = real_c( xlength ) * real_t(0.5) + walberla::math::realRandom( real_t(-0.5), real_t(0.5) );
          yParticle = real_c( ylength ) * real_t(0.5) + walberla::math::realRandom( real_t(-0.5), real_t(0.5) );
 
diff --git a/src/core/debug/CheckFunctions.h b/src/core/debug/CheckFunctions.h
index 84130a8d1..4e64bbf35 100644
--- a/src/core/debug/CheckFunctions.h
+++ b/src/core/debug/CheckFunctions.h
@@ -28,7 +28,6 @@
 #include "core/math/Utility.h"
 #include "core/VectorTrait.h"
 
-//#include <boost/math/special_functions/next.hpp>
 #include <boost/type_traits/is_arithmetic.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/type_traits/has_left_shift.hpp>
diff --git a/src/core/math/DistributedSample.cpp b/src/core/math/DistributedSample.cpp
index 8d7a531b5..509af72af 100644
--- a/src/core/math/DistributedSample.cpp
+++ b/src/core/math/DistributedSample.cpp
@@ -26,7 +26,6 @@
 
 #include <boost/algorithm/string/replace.hpp>
 #include <boost/lexical_cast.hpp>
-#include <boost/math/special_functions/pow.hpp>
 
 
 
@@ -64,7 +63,10 @@ void DistributedSample::mpiAllGather()
    variance_ = real_t(0);
 
    for( auto it = data_.begin(); it != data_.end(); ++it )
-      variance_ += boost::math::pow<2>( *it - mean_ );
+   {
+      real_t val = *it - mean_;
+      variance_ += val*val;
+   }
 
    WALBERLA_MPI_SECTION()
    {
@@ -120,7 +122,10 @@ void DistributedSample::mpiGather( int rank )
    variance_ = real_t(0);
 
    for( auto it = data_.begin(); it != data_.end(); ++it )
-      variance_ += boost::math::pow<2>( *it - mean_ );
+   {
+      real_t val = *it - mean_;
+      variance_ += val*val;
+   }
 
    WALBERLA_MPI_SECTION()
    {
diff --git a/src/core/math/FPClassify.h b/src/core/math/FPClassify.h
index 2c87b8e2b..abd8aee27 100644
--- a/src/core/math/FPClassify.h
+++ b/src/core/math/FPClassify.h
@@ -22,7 +22,7 @@
 
 #pragma once
 
-#include <boost/math/special_functions/fpclassify.hpp>
+#include <cmath>
 
 
 namespace walberla {
@@ -39,7 +39,7 @@ namespace math {
 template<typename T>
 inline bool isnan(T x)
 {
-   return (boost::math::isnan)(x);
+   return (std::isnan)(x);
 }
 
 
@@ -53,7 +53,7 @@ inline bool isnan(T x)
 template<typename T>
 inline bool isinf(T x)
 {
-   return (boost::math::isinf)(x);
+   return (std::isinf)(x);
 }
 
 
@@ -67,7 +67,7 @@ inline bool isinf(T x)
 template<typename T>
 bool finite(T x)
 {
-    return (boost::math::isfinite)(x);
+    return (std::isfinite)(x);
 }
 
 
diff --git a/src/core/math/GenericAABB.h b/src/core/math/GenericAABB.h
index 99e1ac647..aa8e15ac5 100644
--- a/src/core/math/GenericAABB.h
+++ b/src/core/math/GenericAABB.h
@@ -30,7 +30,7 @@
 #include <boost/array.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
 
-#include <boost/random/uniform_real_distribution.hpp>
+#include <random>
 
 
 namespace walberla {
diff --git a/src/core/math/GenericAABB.impl.h b/src/core/math/GenericAABB.impl.h
index e5db1a23c..5709d770e 100644
--- a/src/core/math/GenericAABB.impl.h
+++ b/src/core/math/GenericAABB.impl.h
@@ -1717,7 +1717,7 @@ void GenericAABB< T >::intersect( const GenericAABB & other )
  * The point is in ( [ xMin(), xMax() ), [ yMin(), yMax() ), [ zMin(), zMax() ) )
  *
  * \pre !empty()
- * \param engine  An Uniform Random Number Generator (e.g. boost::random::mt19937)
+ * \param engine  An Uniform Random Number Generator (e.g. std::mt19937)
  * \returns Random point within *this
  */
 template< typename T >
@@ -1725,9 +1725,9 @@ template< typename Engine >
 typename GenericAABB< T >::vector_type GenericAABB< T >::randomPoint( Engine & engine ) const
 {
    WALBERLA_ASSERT( !empty() );
-   boost::random::uniform_real_distribution< T > randX( xMin(), xMax() );
-   boost::random::uniform_real_distribution< T > randY( yMin(), yMax() );
-   boost::random::uniform_real_distribution< T > randZ( zMin(), zMax() );
+   std::uniform_real_distribution< T > randX( xMin(), xMax() );
+   std::uniform_real_distribution< T > randY( yMin(), yMax() );
+   std::uniform_real_distribution< T > randZ( zMin(), zMax() );
 
    return vector_type( randX( engine ), randY( engine ), randZ( engine ) );
 }
diff --git a/src/core/math/Random.cpp b/src/core/math/Random.cpp
index 0236e1012..72707f208 100644
--- a/src/core/math/Random.cpp
+++ b/src/core/math/Random.cpp
@@ -29,9 +29,9 @@ namespace math {
 
 namespace internal {
 
-static boost::mt19937 generator; // static boost::random::mt19937_64 generator;
+static std::mt19937 generator; // static std::mt19937_64 generator;
 
-boost::mt19937 & getGenerator() // boost::random::mt19937_64
+std::mt19937 & getGenerator() // std::mt19937_64
 {
    return generator;
 }
@@ -40,7 +40,7 @@ boost::mt19937 & getGenerator() // boost::random::mt19937_64
 
 
 
-void seedRandomGenerator( const boost::mt19937::result_type & seed )
+void seedRandomGenerator( const std::mt19937::result_type & seed )
 {
 #ifdef _OPENMP
    #pragma omp critical (random)
diff --git a/src/core/math/Random.h b/src/core/math/Random.h
index 774764615..517c4d4be 100644
--- a/src/core/math/Random.h
+++ b/src/core/math/Random.h
@@ -24,9 +24,7 @@
 #include "core/DataTypes.h"
 #include "core/debug/Debug.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>  // #include <boost/random/uniform_int_distribution.hpp>
-#include <boost/random/uniform_real.hpp> // #include <boost/random/uniform_real_distribution.hpp>
+#include <random>
 #include <limits>
 
 
@@ -36,12 +34,12 @@ namespace math {
 
 
 namespace internal {
-boost::mt19937 & getGenerator(); // boost::random::mt19937_64
+std::mt19937 & getGenerator(); // std::mt19937_64
 }
 
 
 
-void seedRandomGenerator( const boost::mt19937::result_type & seed ); // boost::random::mt19937_64
+void seedRandomGenerator( const std::mt19937::result_type & seed ); // std::mt19937_64
 
 
 
@@ -52,12 +50,13 @@ void seedRandomGenerator( const boost::mt19937::result_type & seed ); // boost::
 //**********************************************************************************************************************
 template< typename INT >
 INT intRandom( const INT min = std::numeric_limits<INT>::min(), const INT max = std::numeric_limits<INT>::max(),
-               boost::mt19937 & generator = internal::getGenerator() )
+               std::mt19937 & generator = internal::getGenerator() )
 {
    static_assert_int_t< INT >();
+   static_assert(sizeof(INT) > sizeof(char), "cannot use char");
    WALBERLA_ASSERT_LESS_EQUAL( min, max );
 
-   boost::uniform_int< INT > distribution( min, max ); // boost::random::uniform_int_distribution< INT > distribution;
+   std::uniform_int_distribution< INT > distribution( min, max );
 
    INT value;
 #ifdef _OPENMP
@@ -68,19 +67,70 @@ INT intRandom( const INT min = std::numeric_limits<INT>::min(), const INT max =
    return value;
 }
 
+template<>
+inline char intRandom<char>( const char min, const char max, std::mt19937 & generator )
+{
+   static_assert_int_t< char >();
+   WALBERLA_ASSERT_LESS_EQUAL( min, max );
+
+   std::uniform_int_distribution< int16_t > distribution( min, max );
+
+   char value;
+#ifdef _OPENMP
+   #pragma omp critical (random)
+#endif
+   { value = static_cast<char>( distribution( generator ) ); }
+
+   return value;
+}
+
+template<>
+inline unsigned char intRandom<unsigned char>( const unsigned char min, const unsigned char max, std::mt19937 & generator )
+{
+   static_assert_int_t< unsigned char >();
+   WALBERLA_ASSERT_LESS_EQUAL( min, max );
+
+   std::uniform_int_distribution< int16_t > distribution( min, max );
+
+   unsigned char value;
+#ifdef _OPENMP
+   #pragma omp critical (random)
+#endif
+   { value = static_cast<unsigned char>( distribution( generator ) ); }
+
+   return value;
+}
+
+template<>
+inline signed char intRandom<signed char>( const signed char min, const signed char max, std::mt19937 & generator )
+{
+   static_assert_int_t< signed char >();
+   WALBERLA_ASSERT_LESS_EQUAL( min, max );
+
+   std::uniform_int_distribution< int16_t > distribution( min, max );
+
+   signed char value;
+#ifdef _OPENMP
+   #pragma omp critical (random)
+#endif
+   { value = static_cast<signed char>( distribution( generator ) ); }
+
+   return value;
+}
+
 
 
 template< typename INT >
 class IntRandom
 {
 public:
-   IntRandom( const boost::mt19937::result_type & seed = boost::mt19937::result_type() ) { generator_.seed( seed ); }
+   IntRandom( const std::mt19937::result_type & seed = std::mt19937::result_type() ) { generator_.seed( seed ); }
    INT operator()( const INT min = std::numeric_limits<INT>::min(), const INT max = std::numeric_limits<INT>::max() )
    {
       return intRandom( min, max, generator_ );
    }
 private:
-   boost::mt19937 generator_;
+   std::mt19937 generator_;
 };
 
 
@@ -91,12 +141,12 @@ private:
 */
 //**********************************************************************************************************************
 template< typename REAL >
-REAL realRandom( const REAL min = REAL(0), const REAL max = REAL(1), boost::mt19937 & generator = internal::getGenerator() )
+REAL realRandom( const REAL min = REAL(0), const REAL max = REAL(1), std::mt19937 & generator = internal::getGenerator() )
 {
    static_assert( std::numeric_limits<REAL>::is_specialized && !std::numeric_limits<REAL>::is_integer, "Floating point type required/expected!" );
    WALBERLA_ASSERT_LESS( min, max );
 
-   boost::uniform_real< REAL > distribution( min, max ); // boost::uniform_real_distribution< REAL > distribution( min, max );
+   std::uniform_real_distribution< REAL > distribution( min, max );
 
    REAL value;
 #ifdef _OPENMP
@@ -113,13 +163,13 @@ template< typename REAL >
 class RealRandom
 {
 public:
-   RealRandom( const boost::mt19937::result_type & seed = boost::mt19937::result_type() ) { generator_.seed( seed ); }
+   RealRandom( const std::mt19937::result_type & seed = std::mt19937::result_type() ) { generator_.seed( seed ); }
    REAL operator()( const REAL min = REAL(0), const REAL max = REAL(1) )
    {
       return realRandom( min, max, generator_ );
    }
 private:
-   boost::mt19937 generator_;
+   std::mt19937 generator_;
 };
 
 
@@ -133,10 +183,10 @@ inline bool boolRandom() { ///< Randomly returns 'true' or 'false'
 class BoolRandom
 {
 public:
-   BoolRandom( const boost::mt19937::result_type & seed = boost::mt19937::result_type() ) { generator_.seed( seed ); }
+   BoolRandom( const std::mt19937::result_type & seed = std::mt19937::result_type() ) { generator_.seed( seed ); }
    bool operator()() { return boolRandom(); }
 private:
-   boost::mt19937 generator_;
+   std::mt19937 generator_;
 };
 
 
diff --git a/src/core/math/Sample.cpp b/src/core/math/Sample.cpp
index 6db1ad17a..eabe897cd 100644
--- a/src/core/math/Sample.cpp
+++ b/src/core/math/Sample.cpp
@@ -26,7 +26,6 @@
 
 #include <boost/algorithm/string/replace.hpp>
 #include <boost/lexical_cast.hpp>
-#include <boost/math/special_functions/pow.hpp>
 
 #include <functional>
 #include <iterator>
@@ -128,7 +127,10 @@ real_t Sample::variance( real_t theMean ) const
 
    KahanAccumulator< real_t > acc;
    for(auto it = begin(); it != end(); ++it)
-      acc += boost::math::pow<2>(*it - theMean);
+   {
+      real_t val = *it - theMean;
+      acc += val*val;
+   }
       
    return acc.get() / real_c(size());
 }
diff --git a/src/core/math/equation_system/Equation.cpp b/src/core/math/equation_system/Equation.cpp
index 6c355b86f..53bef41aa 100644
--- a/src/core/math/equation_system/Equation.cpp
+++ b/src/core/math/equation_system/Equation.cpp
@@ -23,7 +23,8 @@
 #include "Operator.h"
 #include "Variable.h"
 
-#include <boost/math/special_functions/fpclassify.hpp>
+#include <cmath>
+#include <algorithm>
 
 
 namespace walberla {
@@ -128,10 +129,10 @@ namespace math {
       double left  = root_->left_->compute();
       double right = root_->right_->compute();
 
-      if ( boost::math::isnan(left) && boost::math::isnan(right) ){
+      if ( std::isnan(left) && std::isnan(right) ){
          //WALBERLA_LOG_WARNING( "WARNING: Both values are NAN -> return true" );
          return true;
-      } else if ( boost::math::isinf(left) && boost::math::isinf(right) ){
+      } else if ( std::isinf(left) && std::isinf(right) ){
     	 //WALBERLA_LOG_WARNING( "WARNING: Both values are INF -> return true" );
          return true;
       }
diff --git a/src/core/math/equation_system/Variable.cpp b/src/core/math/equation_system/Variable.cpp
index 12295cb9c..04d9b7517 100644
--- a/src/core/math/equation_system/Variable.cpp
+++ b/src/core/math/equation_system/Variable.cpp
@@ -21,7 +21,7 @@
 
 #include "Variable.h"
 
-#include <boost/math/special_functions/fpclassify.hpp>
+#include <cmath>
 #include <sstream>
 
 
@@ -36,7 +36,7 @@ namespace math {
 
    void Var::setValue( const double value ){
       value_ = value;
-      valid_ = !boost::math::isnan( value );
+      valid_ = !std::isnan( value );
    }
 
    bool Var::operator==( const Var& var) const {
diff --git a/src/domain_decomposition/IBlock.h b/src/domain_decomposition/IBlock.h
index e9b4ec6fb..215f8a51d 100644
--- a/src/domain_decomposition/IBlock.h
+++ b/src/domain_decomposition/IBlock.h
@@ -30,7 +30,6 @@
 #include "core/math/AABB.h"
 #include "core/uid/SUID.h"
 
-#include <boost/cast.hpp>
 #include <typeinfo>
 #include <vector>
 
diff --git a/src/gather/CurveGatherPackInfo.impl.h b/src/gather/CurveGatherPackInfo.impl.h
index c18ff74d5..0039b2a05 100644
--- a/src/gather/CurveGatherPackInfo.impl.h
+++ b/src/gather/CurveGatherPackInfo.impl.h
@@ -20,7 +20,6 @@
 //
 //======================================================================================================================
 
-#include "boost/lambda/bind.hpp"
 #include "core/debug/Debug.h"
 #include "core/logging/Logging.h"
 #include "core/math/Parser.h"
diff --git a/tests/blockforest/BlockDataIOTest.cpp b/tests/blockforest/BlockDataIOTest.cpp
index 6a550599a..291b3322a 100644
--- a/tests/blockforest/BlockDataIOTest.cpp
+++ b/tests/blockforest/BlockDataIOTest.cpp
@@ -91,7 +91,7 @@ int main( int argc, char* argv[] )
    auto dataHandling = make_shared< field::DefaultBlockDataHandling< FieldType > >( sbf, uint_t(3), 0.0, field::zyxf );
    auto originalFieldId = sbf->addBlockData( dataHandling, "OriginalField", None, Empty );
    
-   math::seedRandomGenerator( numeric_cast<boost::mt19937::result_type>( MPIManager::instance()->rank() ) );
+   math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( MPIManager::instance()->rank() ) );
 
    for( auto it = sbf->begin( None, Empty ); it != sbf->end(); ++it )
    {
diff --git a/tests/blockforest/StructuredBlockForestTest.cpp b/tests/blockforest/StructuredBlockForestTest.cpp
index 84032f809..fc9ce847c 100644
--- a/tests/blockforest/StructuredBlockForestTest.cpp
+++ b/tests/blockforest/StructuredBlockForestTest.cpp
@@ -25,7 +25,7 @@
 #include "core/debug/TestSubsystem.h"
 #include "core/mpi/Environment.h"
 
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 
 namespace walberla {
 namespace blockforest {
@@ -418,7 +418,7 @@ static void test() {
    bb = *( block->getData< CellInterval >( bid ) );
    WALBERLA_CHECK_EQUAL( bb, CellInterval(Cell(150,180,210),Cell(199,239,279)) );
 
-   boost::random::mt19937 eng( 23 );
+   std::mt19937 eng( 23 );
    for( int i = 0; i < 23; ++i )
    {
       const Vector3<real_t> globalPoint = forest.getDomain().getScaled( real_t(1.25) ).randomPoint( eng );
diff --git a/tests/core/cell/CellIntervalTest.cpp b/tests/core/cell/CellIntervalTest.cpp
index bf4f42da5..5ce0c57f2 100644
--- a/tests/core/cell/CellIntervalTest.cpp
+++ b/tests/core/cell/CellIntervalTest.cpp
@@ -23,19 +23,20 @@
 #include "core/cell/CellInterval.h"
 #include "core/debug/TestSubsystem.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>
+#include <random>
 #include <iostream>
 #include <iterator>
 
 
 using namespace walberla;
 
+typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+
 CellInterval makeRandomInterval(uint_t maxSize)
 {
-   static boost::mt11213b rng;
-   boost::uniform_int<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() - cell_idx_c( maxSize ) );
-   boost::uniform_int<uint_t> dist2( uint_t(0), maxSize );
+   static mt11213b rng;
+   std::uniform_int_distribution<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() - cell_idx_c( maxSize ) );
+   std::uniform_int_distribution<uint_t> dist2( uint_t(0), maxSize );
 
    cell_idx_t xMin = dist(rng);
    cell_idx_t yMin = dist(rng);
@@ -50,10 +51,10 @@ CellInterval makeRandomInterval(uint_t maxSize)
 
 CellInterval makeRandomEmptyInterval(uint_t maxSize)
 {
-   static boost::mt11213b rng;
-   boost::uniform_int<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min() + cell_idx_c( maxSize ), std::numeric_limits<cell_idx_t>::max() - cell_idx_c( maxSize ) );
-   boost::uniform_int<uint_t> dist2( uint_t(1), maxSize );
-   boost::uniform_int<uint_t> dist3( uint_t(0), uint_t(1) );
+   static mt11213b rng;
+   std::uniform_int_distribution<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min() + cell_idx_c( maxSize ), std::numeric_limits<cell_idx_t>::max() - cell_idx_c( maxSize ) );
+   std::uniform_int_distribution<uint_t> dist2( uint_t(1), maxSize );
+   std::uniform_int_distribution<uint_t> dist3( uint_t(0), uint_t(1) );
 
    cell_idx_t xMin = dist(rng);
    cell_idx_t yMin = dist(rng);
@@ -81,8 +82,8 @@ CellInterval makeRandomEmptyInterval(uint_t maxSize)
 
 Cell makeRandomCell()
 {
-   static boost::mt11213b rng;
-   boost::uniform_int<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() );
+   static mt11213b rng;
+   std::uniform_int_distribution<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() );
    return Cell( dist(rng), dist(rng), dist(rng) );
 }
 
diff --git a/tests/core/cell/CellTest.cpp b/tests/core/cell/CellTest.cpp
index 1cdc57474..55c6097ba 100644
--- a/tests/core/cell/CellTest.cpp
+++ b/tests/core/cell/CellTest.cpp
@@ -23,8 +23,7 @@
 #include "core/cell/Cell.h"
 #include "core/debug/TestSubsystem.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>
+#include <random>
 
 
 using namespace walberla;
@@ -135,8 +134,9 @@ int main( int /*argc*/, char** /*argv*/ ) {
 
    debug::enterTestMode();
 
-   boost::mt11213b rng;
-   boost::uniform_int<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() );
+   typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+   mt11213b rng;
+   std::uniform_int_distribution<cell_idx_t> dist( std::numeric_limits<cell_idx_t>::min(), std::numeric_limits<cell_idx_t>::max() );
 
    for(int i = 0; i < 100000; ++i)
    {
diff --git a/tests/core/load_balancing/MetisTest.cpp b/tests/core/load_balancing/MetisTest.cpp
index a69cffb99..a460077c2 100644
--- a/tests/core/load_balancing/MetisTest.cpp
+++ b/tests/core/load_balancing/MetisTest.cpp
@@ -40,7 +40,6 @@
 #include "vtk/BlockCellDataWriter.h"
 
 #include <boost/lexical_cast.hpp>
-#include <boost/random.hpp>
 
 int main( int argc, char * argv[] )
 {
diff --git a/tests/core/load_balancing/ParMetisTest.cpp b/tests/core/load_balancing/ParMetisTest.cpp
index ce43ce97b..320315e8d 100644
--- a/tests/core/load_balancing/ParMetisTest.cpp
+++ b/tests/core/load_balancing/ParMetisTest.cpp
@@ -42,7 +42,6 @@
 #include "vtk/BlockCellDataWriter.h"
 
 #include <boost/lexical_cast.hpp>
-#include <boost/random.hpp>
 
 int main( int argc, char * argv[] )
 {
diff --git a/tests/core/math/GenericAABBTest.cpp b/tests/core/math/GenericAABBTest.cpp
index 6a5d4ce9a..9194c0ba8 100644
--- a/tests/core/math/GenericAABBTest.cpp
+++ b/tests/core/math/GenericAABBTest.cpp
@@ -27,8 +27,7 @@
 
 #include "stencil/D3CornerStencil.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_real.hpp>
+#include <random>
 
 
 using namespace walberla;
@@ -112,7 +111,7 @@ void testNonEmptyAABB( const GenericAABB< T > & aabb )
    WALBERLA_CHECK_FLOAT_EQUAL( tmpAABB.volume(), T(0) );
    WALBERLA_CHECK_IDENTICAL( tmpAABB.volume(), aabb.intersectionVolume( intersectingBox ) );
 
-   boost::random::mt19937 urng;
+   std::mt19937 urng;
    for( int i = 0; i < 100; ++i )
    {
       auto p = aabb.randomPoint( urng );
@@ -530,9 +529,10 @@ void testConstructors( const T x0, const T y0, const T z0, const T x1, const T y
 template< typename T >
 void randomTest()
 {
-   boost::mt11213b rng;
+   typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+   mt11213b rng;
 
-   boost::uniform_real<T> dist( -T(10), T(10) );
+   std::uniform_real_distribution<T> dist( -T(10), T(10) );
 
    for( int i = 0; i < 1000; ++i )
    {
@@ -600,7 +600,7 @@ void testAABBDistancesRandom( const GenericAABB< T > & baseAABB )
 {
    static const uint_t NUM_BOXES  = 100;
    static const uint_t NUM_POINTS = 1000;
-   boost::random::mt19937 rng;  
+   std::mt19937 rng;
 
    for( uint_t i = 0; i < NUM_BOXES; ++i )
    {
diff --git a/tests/core/math/PlaneTest.cpp b/tests/core/math/PlaneTest.cpp
index b74f2c30f..b8f2d49b7 100644
--- a/tests/core/math/PlaneTest.cpp
+++ b/tests/core/math/PlaneTest.cpp
@@ -25,9 +25,7 @@
 #include "core/math/Vector3.h"
 #include "core/mpi/Environment.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 #include <cmath>
 #include <vector>
 
@@ -39,9 +37,9 @@ template < typename scalar_t >
 struct RandomPointGenerator
 {
    typedef walberla::Vector3<scalar_t> vector_t;
-   typedef boost::mt11213b RandomNumberEngine;
-   typedef boost::normal_distribution<scalar_t> NormalDistribution;
-   typedef boost::variate_generator< RandomNumberEngine, NormalDistribution > Generator;
+   typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+   typedef mt11213b RandomNumberEngine;
+   typedef std::normal_distribution<scalar_t> NormalDistribution;
 
    RandomPointGenerator( const vector_t & mu, const vector_t & sigma )
    {
@@ -49,17 +47,19 @@ struct RandomPointGenerator
       {
          RandomNumberEngine rne;
          rne.seed( numeric_cast< RandomNumberEngine::result_type >( i * 642573 ) );
-         normalDistributions.push_back( Generator( rne, NormalDistribution( mu[i], sigma[i] ) ) );
+         engines.push_back( rne );
+         normalDistributions.push_back( NormalDistribution( mu[i], sigma[i] ) );
       }
    }
 
    vector_t operator()()
    {
-      return vector_t( normalDistributions[0](), normalDistributions[1](), normalDistributions[2]() );
+      return vector_t( normalDistributions[0](engines[0]), normalDistributions[1](engines[1]), normalDistributions[2](engines[2]) );
    }
 
 private:
-   std::vector< Generator > normalDistributions;
+   std::vector< RandomNumberEngine > engines;
+   std::vector< NormalDistribution > normalDistributions;
 };
 
 void testIOStream( const Plane & p )
diff --git a/tests/core/math/PrimesTest.cpp b/tests/core/math/PrimesTest.cpp
index 9e8a776a0..141a9704b 100644
--- a/tests/core/math/PrimesTest.cpp
+++ b/tests/core/math/PrimesTest.cpp
@@ -27,8 +27,7 @@
 #include "core/math/Primes.h"
 #include "core/mpi/MPIManager.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>
+#include <random>
 #include <algorithm>
 #include <numeric>
 
@@ -85,8 +84,9 @@ int main(int argc, char * argv[])
       runTests( n );
    }
 
-   boost::mt11213b rng;
-   boost::uniform_int<uint_t> dist( 100, 10000 );
+   typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+   mt11213b rng;
+   std::uniform_int_distribution<uint_t> dist( 100, 10000 );
    for(int i = 0; i < 100; ++i)
    {
       const uint_t n = dist(rng);
diff --git a/tests/core/mpi/BufferSystemTest.cpp b/tests/core/mpi/BufferSystemTest.cpp
index 312de8799..a6104a87e 100644
--- a/tests/core/mpi/BufferSystemTest.cpp
+++ b/tests/core/mpi/BufferSystemTest.cpp
@@ -26,9 +26,7 @@
 #include "core/mpi/BufferSystem.h"
 #include "core/mpi/Environment.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 #include <cmath>
 #include <iostream>
@@ -43,7 +41,7 @@ using namespace std::literals::chrono_literals;
 
 
 
-typedef boost::mt19937 base_generator_type;
+typedef std::mt19937 base_generator_type;
 
 /**
  * Utility function for sleeping a random time
@@ -60,10 +58,9 @@ void randomSleep( int maxTimeInMs = 20 )
    unsigned int seed = static_cast<unsigned int>(std::time(0)) + static_cast<unsigned int>(rank*1000) + counter;
    generator.seed(seed);
 
-   boost::uniform_int<> uni_dist(0,maxTimeInMs);
-   boost::variate_generator<base_generator_type&, boost::uniform_int<> > uni(generator, uni_dist);
+   std::uniform_int_distribution<> uni_dist(0,maxTimeInMs);
 
-   int sleepTime = uni();
+   int sleepTime = uni_dist(generator);
    std::this_thread::sleep_for( sleepTime * 1ms );
 }
 
diff --git a/tests/core/mpi/BufferTest.cpp b/tests/core/mpi/BufferTest.cpp
index 4d1819dcd..34f961f1a 100644
--- a/tests/core/mpi/BufferTest.cpp
+++ b/tests/core/mpi/BufferTest.cpp
@@ -36,8 +36,7 @@
 #include "core/mpi/RecvBuffer.h"
 #include "core/mpi/SendBuffer.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>  // #include <boost/random/uniform_int_distribution.hpp>
+#include <random>
 #include <cstring>
 #include <iostream>
 
@@ -48,9 +47,9 @@ using namespace mpi;
 template<typename T>
 void initIntegerContainer( T & container )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<typename T::value_type> dist;
-   boost::uniform_int<size_t> size_dist(10,1000);
+   static std::mt19937 rng;
+   std::uniform_int_distribution<typename T::value_type> dist;
+   std::uniform_int_distribution<size_t> size_dist(10,1000);
 
    size_t size = size_dist(rng);
    container.clear();
@@ -61,9 +60,9 @@ void initIntegerContainer( T & container )
 template<typename T>
 void initIntegerAssocContainer( T & container )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<typename T::value_type> dist;
-   boost::uniform_int<size_t> size_dist(10,1000);
+   static std::mt19937 rng;
+   std::uniform_int_distribution<typename T::value_type> dist;
+   std::uniform_int_distribution<size_t> size_dist(10,1000);
 
    size_t size = size_dist(rng);
    container.clear();
@@ -74,9 +73,9 @@ void initIntegerAssocContainer( T & container )
 
 void initVecBool( std::vector<bool> & vecBool )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<walberla::uint32_t> dist;
-   boost::uniform_int<size_t> size_dist(10,1000);
+   static std::mt19937 rng;
+   std::uniform_int_distribution<walberla::uint32_t> dist;
+   std::uniform_int_distribution<size_t> size_dist(10,1000);
 
    size_t size = size_dist(rng);
    vecBool.clear();
@@ -88,10 +87,10 @@ void initVecBool( std::vector<bool> & vecBool )
 template<typename T>
 void initIntegerMap( T & container )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<typename T::mapped_type> mapped_dist;
-   boost::uniform_int<typename T::key_type> key_dist;
-   boost::uniform_int<size_t> size_dist(10,1000);
+   static std::mt19937 rng;
+   std::uniform_int_distribution<typename T::mapped_type> mapped_dist;
+   std::uniform_int_distribution<typename T::key_type> key_dist;
+   std::uniform_int_distribution<size_t> size_dist(10,1000);
 
    size_t size = size_dist(rng);
    container.clear();
@@ -103,9 +102,9 @@ void initIntegerMap( T & container )
 template<typename T>
 void initCellContainer( T & container )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<cell_idx_t> dist;
-   boost::uniform_int<size_t> size_dist(10,1000);
+   static std::mt19937 rng;
+   std::uniform_int_distribution<cell_idx_t> dist;
+   std::uniform_int_distribution<size_t> size_dist(10,1000);
 
    size_t size = size_dist(rng);
    container.clear();
@@ -116,8 +115,8 @@ void initCellContainer( T & container )
 template<typename T, std::size_t N>
 void initBoostArray( boost::array< T, N > & array )
 {
-   static boost::mt19937 rng;
-   boost::uniform_int<T> dist;
+   static std::mt19937 rng;
+   std::uniform_int_distribution<T> dist;
 
    for( auto it = array.begin(); it != array.end(); ++it )
       *it = dist( rng );
diff --git a/tests/cuda/communication/GPUPackInfoCommunicationTest.cpp b/tests/cuda/communication/GPUPackInfoCommunicationTest.cpp
index 0d118f1f3..8d0d66c8e 100644
--- a/tests/cuda/communication/GPUPackInfoCommunicationTest.cpp
+++ b/tests/cuda/communication/GPUPackInfoCommunicationTest.cpp
@@ -131,7 +131,7 @@ int main( int argc, char ** argv )
       BlockDataID asyncGPUFieldId = blocks->addStructuredBlockData< GPUFieldType >( &createGPUField,
                                                                                     "asyncGPUField" );
 
-      math::seedRandomGenerator( numeric_cast<boost::mt19937::result_type>( MPIManager::instance()->rank() ) );
+      math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( MPIManager::instance()->rank() ) );
       // Initialize CPU field with random values
       initFields( blocks, sourceFieldId );
 
diff --git a/tests/field/AccuracyEvaluationTest.cpp b/tests/field/AccuracyEvaluationTest.cpp
index 8fedfe317..c7656509f 100644
--- a/tests/field/AccuracyEvaluationTest.cpp
+++ b/tests/field/AccuracyEvaluationTest.cpp
@@ -67,8 +67,8 @@ int main( int argc, char* argv[] )
                                                       real_t(1), // dx
                                                       uint_t( 2), uint_t( 1), uint_t( 2) ); // number of processes
 
-   //math::seedRandomGenerator( numeric_cast<boost::mt19937::result_type>( std::time(0) ) );
-   math::seedRandomGenerator( numeric_cast<boost::mt19937::result_type>( MPIManager::instance()->rank() ) );
+   //math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( std::time(0) ) );
+   math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( MPIManager::instance()->rank() ) );
 
    auto sId = field::addToStorage< ScalarField_T >( blocks, "scalar field" );
    auto vId = field::addToStorage< VectorField_T >( blocks, "vector field" );
diff --git a/tests/field/FieldFileIOTest.cpp b/tests/field/FieldFileIOTest.cpp
index 2e33e0dce..6e7a95a2d 100644
--- a/tests/field/FieldFileIOTest.cpp
+++ b/tests/field/FieldFileIOTest.cpp
@@ -114,7 +114,7 @@ int main( int argc, char* argv[] )
    auto originalFieldId = field::addToStorage< FieldType >( sbf, "OriginalField" );
    auto readFieldId     = field::addToStorage< FieldType >( sbf, "ReadField" );
 
-   math::seedRandomGenerator( numeric_cast<boost::mt19937::result_type>( MPIManager::instance()->rank() ) );
+   math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( MPIManager::instance()->rank() ) );
 
    for( auto it = sbf->begin(); it != sbf->end(); ++it )
    {
diff --git a/tests/field/FieldMPIDatatypesTest.cpp b/tests/field/FieldMPIDatatypesTest.cpp
index f5411aeb5..3328c1bd7 100644
--- a/tests/field/FieldMPIDatatypesTest.cpp
+++ b/tests/field/FieldMPIDatatypesTest.cpp
@@ -33,9 +33,7 @@
 
 #include "stencil/D3Q27.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int_distribution.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
+#include <random>
 
 namespace mpi_datatypes_test {
 
@@ -51,7 +49,7 @@ public:
    template< typename T, uint_t fSize >
    void operator()( Field< T, fSize > & field )
    {
-      boost::random::uniform_real_distribution< T > distribution;
+      std::uniform_real_distribution< T > distribution;
 
       for( auto it = field.begin(); it != field.end(); ++it )
       {
@@ -62,7 +60,7 @@ public:
    template< typename T, uint_t fSize >
    void operator()( GhostLayerField< T, fSize > & field )
    {
-      boost::random::uniform_real_distribution< T > distribution;
+      std::uniform_real_distribution< T > distribution;
 
       for( auto it = field.beginWithGhostLayer(); it != field.end(); ++it )
       {
@@ -71,7 +69,7 @@ public:
    }
 
 private:
-   boost::random::mt19937 generator_;
+   std::mt19937 generator_;
 };
 
 
diff --git a/tests/geometry/RandomPointGenerator.h b/tests/geometry/RandomPointGenerator.h
index 8275946f4..097f610c4 100644
--- a/tests/geometry/RandomPointGenerator.h
+++ b/tests/geometry/RandomPointGenerator.h
@@ -23,9 +23,7 @@
 
 #include "core/math/Vector3.h"
 
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 #include <vector>
 
 
@@ -33,21 +31,25 @@ template < typename scalar_t >
 struct RandomPointGenerator
 {
    typedef walberla::Vector3<scalar_t> vector_t;
-   typedef boost::mt11213b RandomNumberEngine;
-   typedef boost::normal_distribution<scalar_t> NormalDistribution;
-   typedef boost::variate_generator< RandomNumberEngine, NormalDistribution > Generator;
+   typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+   typedef mt11213b RandomNumberEngine;
+   typedef std::normal_distribution<scalar_t> NormalDistribution;
 
    RandomPointGenerator( const vector_t & mu, const vector_t & sigma )
    {
       for( walberla::uint_t i = 0; i < 3; ++i )
-         normalDistributions.push_back( Generator( RandomNumberEngine(), NormalDistribution( mu[i], sigma[i] ) ) );
+      {
+         engines.push_back( RandomNumberEngine() );
+         normalDistributions.push_back( NormalDistribution( mu[i], sigma[i] ) );
+      }
    }
 
    vector_t operator()()
    {
-      return vector_t( normalDistributions[0](), normalDistributions[1](), normalDistributions[2]() );
+      return vector_t( normalDistributions[0](engines[0]), normalDistributions[1](engines[1]), normalDistributions[2](engines[2]) );
    }
 
 private:
-   std::vector< Generator > normalDistributions;
+   std::vector< RandomNumberEngine > engines;
+   std::vector< NormalDistribution > normalDistributions;
 };
diff --git a/tests/geometry/VoxelFileTest.cpp b/tests/geometry/VoxelFileTest.cpp
index 3918b08de..8de8c2156 100644
--- a/tests/geometry/VoxelFileTest.cpp
+++ b/tests/geometry/VoxelFileTest.cpp
@@ -54,8 +54,7 @@
 
 #include <boost/filesystem.hpp>
 #include <boost/foreach.hpp>
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_int.hpp>
+#include <random>
 
 #ifdef _MSC_VER
 #  pragma warning(push)
@@ -69,15 +68,17 @@
 #  pragma warning(pop)
 #endif //_MSC_VER
 
+typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
+
 
 /// randomize the memory underlying the vector up the maximal size (== capacity)
 template<typename T>
 void randomizeVector( std::vector<T> & v )
 {
-   //boost::random::mt11213b rng;
-   boost::mt11213b rng;
-   //boost::random::uniform_int_distribution<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
-   boost::uniform_int<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
+   static_assert(sizeof(T) > sizeof(char), "cannot use char");
+
+   mt11213b rng;
+   std::uniform_int_distribution<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
 
    size_t oldSize = v.size();
    v.resize( v.capacity() );
@@ -86,18 +87,64 @@ void randomizeVector( std::vector<T> & v )
    v.resize(oldSize);
 }
 
+template<>
+void randomizeVector( std::vector<unsigned char> & v )
+{
+   mt11213b rng;
+   std::uniform_int_distribution<walberla::int16_t> dist( std::numeric_limits<unsigned char>::min(), std::numeric_limits<unsigned char>::max() );
+
+   size_t oldSize = v.size();
+   v.resize( v.capacity() );
+   for(typename std::vector<unsigned char>::iterator it = v.begin(); it != v.end(); ++it)
+      *it = static_cast<unsigned char>( dist(rng) );
+   v.resize(oldSize);
+}
+
+template<>
+void randomizeVector( std::vector<char> & v )
+{
+   mt11213b rng;
+   std::uniform_int_distribution<int16_t> dist( std::numeric_limits<char>::min(), std::numeric_limits<char>::max() );
+
+   size_t oldSize = v.size();
+   v.resize( v.capacity() );
+   for(typename std::vector<char>::iterator it = v.begin(); it != v.end(); ++it)
+      *it = static_cast<char>( dist(rng) );
+   v.resize(oldSize);
+}
+
 template<typename T>
 void makeRandomMultiArray( boost::multi_array<T, 3> & ma)
 {
-   //boost::random::mt11213b rng;
-   boost::mt11213b rng;
-   //boost::random::uniform_int_distribution<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
-   boost::uniform_int<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
+   static_assert(sizeof(T) > sizeof(char), "cannot use char");
+
+   mt11213b rng;
+   std::uniform_int_distribution<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
 
    for(T* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
       *it = dist(rng);
 }
 
+template<>
+void makeRandomMultiArray( boost::multi_array<unsigned char, 3> & ma)
+{
+   mt11213b rng;
+   std::uniform_int_distribution<walberla::int16_t> dist( std::numeric_limits<unsigned char>::min(), std::numeric_limits<unsigned char>::max() );
+
+   for(unsigned char* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
+      *it = static_cast<unsigned char>( dist(rng) );
+}
+
+template<>
+void makeRandomMultiArray( boost::multi_array<char, 3> & ma)
+{
+   mt11213b rng;
+   std::uniform_int_distribution<int16_t> dist( std::numeric_limits<char>::min(), std::numeric_limits<char>::max() );
+
+   for(char* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
+      *it = static_cast<char>( dist(rng) );
+}
+
 template<typename T>
 void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t zSize);
 
diff --git a/tests/lbm/geometry/IntersectionRatioTest.cpp b/tests/lbm/geometry/IntersectionRatioTest.cpp
index d02960184..ae7a4a159 100644
--- a/tests/lbm/geometry/IntersectionRatioTest.cpp
+++ b/tests/lbm/geometry/IntersectionRatioTest.cpp
@@ -26,7 +26,7 @@
 #include "geometry/bodies/Sphere.h"
 #include "geometry/bodies/AABBBody.h"
 
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 
 
 namespace walberla {
@@ -77,7 +77,7 @@ void testAABB()
    static const math::Vector3<real_t> UNIT( real_t( 1 ), real_t( 1 ), real_t( 1 ) );
    static const real_t EPSILON = real_t(1e-4);
 
-   boost::random::mt19937 randomEngine;
+   std::mt19937 randomEngine;
 
    std::vector<math::AABB> testAABBs;
    testAABBs.push_back( math::AABB( -UNIT, UNIT ) );
diff --git a/tests/mesh/MeshAABBIntersectionTest.cpp b/tests/mesh/MeshAABBIntersectionTest.cpp
index f950e7609..fc4711134 100644
--- a/tests/mesh/MeshAABBIntersectionTest.cpp
+++ b/tests/mesh/MeshAABBIntersectionTest.cpp
@@ -29,7 +29,7 @@
 #include "mesh/DistanceComputations.h"
 #include "mesh/TriangleMeshes.h"
 
-#include <boost/random.hpp>
+#include <random>
 #include <boost/lexical_cast.hpp>
 
 #include <vector>
@@ -53,7 +53,7 @@ void runTests( const uint_t numAABBs )
 
    WALBERLA_CHECK( isIntersecting( triDist, meshAABB, real_t(0) ) );
 
-   boost::random::mt19937 rng( uint32_t(42) );
+   std::mt19937 rng( uint32_t(42) );
 
    for(uint_t i = 0; i < numAABBs; ++i)
    {
diff --git a/tests/mesh/MeshAABBSelectionTest.cpp b/tests/mesh/MeshAABBSelectionTest.cpp
index 7a95cb4ea..37d421bfe 100644
--- a/tests/mesh/MeshAABBSelectionTest.cpp
+++ b/tests/mesh/MeshAABBSelectionTest.cpp
@@ -31,7 +31,7 @@
 #include "mesh/MeshIO.h"
 #include "mesh/MeshOperations.h"
 
-#include <boost/random.hpp>
+#include <random>
 #include <boost/lexical_cast.hpp>
 
 #include <algorithm>
@@ -61,7 +61,7 @@ int main( int argc, char * argv[] )
 
    auto aabb = computeAABB( *mesh );
 
-   boost::random::mt19937 rng;  
+   std::mt19937 rng;
 
    for( uint_t i = 0; i < numBoxes; ++i )
    {
diff --git a/tests/mesh/MeshContainmentOctreeTest.cpp b/tests/mesh/MeshContainmentOctreeTest.cpp
index 8858790eb..4e922e17a 100644
--- a/tests/mesh/MeshContainmentOctreeTest.cpp
+++ b/tests/mesh/MeshContainmentOctreeTest.cpp
@@ -35,7 +35,7 @@
 #include "mesh/vtk/VTKMeshWriter.h"
 #include "mesh/vtk/CommonDataSources.h"
 
-#include <boost/random.hpp>
+#include <random>
 
 #include <vector>
 #include <string>
@@ -104,7 +104,7 @@ int main( int argc, char * argv[] )
   if( writeVtk )
      containmentOctree.writeVTKOutput( "containment_octree" );
   
-  boost::random::mt19937 rng;  
+  std::mt19937 rng;
   
   for( int i = 0; i < 10000; ++i )
   {
diff --git a/tests/mesh/MeshDistanceCompareTest.cpp b/tests/mesh/MeshDistanceCompareTest.cpp
index ee89527ec..4067a28e8 100644
--- a/tests/mesh/MeshDistanceCompareTest.cpp
+++ b/tests/mesh/MeshDistanceCompareTest.cpp
@@ -32,7 +32,7 @@
 #include "mesh/DistanceComputations.h"
 #include "mesh/MeshIO.h"
 
-#include <boost/random.hpp>
+#include <random>
 
 #include <vector>
 #include <string>
@@ -54,7 +54,7 @@ void testAABBDistance( const Vector3<real_t> & translationVector = Vector3<real_
 
    TriangleDistance<MeshType> triDist( mesh );
 
-   boost::random::mt19937 rng;
+   std::mt19937 rng;
 
    for( int i = 0; i < 10000; ++i )
    {
diff --git a/tests/mesh/MeshDistanceOctreeTest.cpp b/tests/mesh/MeshDistanceOctreeTest.cpp
index 24601416f..3d849c3f6 100644
--- a/tests/mesh/MeshDistanceOctreeTest.cpp
+++ b/tests/mesh/MeshDistanceOctreeTest.cpp
@@ -32,7 +32,7 @@
 #include "mesh/distance_octree/DistanceOctree.h"
 #include "mesh/MeshIO.h"
 
-#include <boost/random.hpp>
+#include <random>
 
 #include <vector>
 #include <string>
@@ -55,7 +55,7 @@ void test( const std::string & meshFile )
 
    //distanceOctree.writeVTKOutput( "distance_octree" );
 
-   boost::random::mt19937 rng;
+   std::mt19937 rng;
 
    for( int i = 0; i < 1000; ++i )
    {
diff --git a/tests/mesh/MeshDistancePlausibilityTest.cpp b/tests/mesh/MeshDistancePlausibilityTest.cpp
index ff99381ff..5ca33efc1 100644
--- a/tests/mesh/MeshDistancePlausibilityTest.cpp
+++ b/tests/mesh/MeshDistancePlausibilityTest.cpp
@@ -45,7 +45,6 @@
 #include "stencil/D3Q27.h"
 
 #include <boost/lexical_cast.hpp>
-#include <boost/random.hpp>
 
 #include <vector>
 #include <string>
diff --git a/tests/mesh/MeshVTKTest.cpp b/tests/mesh/MeshVTKTest.cpp
index 211bffd11..ac6ecd307 100644
--- a/tests/mesh/MeshVTKTest.cpp
+++ b/tests/mesh/MeshVTKTest.cpp
@@ -30,8 +30,6 @@
 #include "mesh/vtk/CommonDataSources.h"
 #include "mesh/vtk/CommonFilters.h"
 
-#include <boost/random.hpp>
-
 #include <vector>
 #include <string>
 
diff --git a/tests/mesh/PeVTKMeshWriterTest.cpp b/tests/mesh/PeVTKMeshWriterTest.cpp
index 0737b841f..a765f30c4 100644
--- a/tests/mesh/PeVTKMeshWriterTest.cpp
+++ b/tests/mesh/PeVTKMeshWriterTest.cpp
@@ -38,10 +38,11 @@
 #include <core/logging/Logging.h>
 #include <core/timing/TimingTree.h>
 #include <core/waLBerlaBuildInfo.h>
+#include <core/math/Utility.h>
 #include <postprocessing/sqlite/SQLite.h>
 #include <vtk/VTKOutput.h>
 
-#include <boost/random.hpp>
+#include <random>
 
 using namespace walberla;
 using namespace walberla::pe;
@@ -86,15 +87,16 @@ std::vector<Vector3<real_t>> generatePointCloudDodecahedron()
 template< typename RNG >
 std::vector<Vector3<real_t>> generatPointCloudOnSphere( const real_t radius, const uint_t numPoints, RNG & rng )
 {
-   boost::uniform_on_sphere<real_t> distribution(3);
+   std::uniform_real_distribution<real_t> distribution;
 
    std::vector<Vector3<real_t>> pointCloud( numPoints );
    for( auto & p : pointCloud )
    {
-      auto v = distribution(rng);
-      p[0] = v[0] * radius;
-      p[1] = v[1] * radius;
-      p[2] = v[2] * radius;
+      real_t theta = 2 * math::PI * distribution(rng);
+      real_t phi = std::acos( real_t(1.0) - real_t(2.0) * distribution(rng) );
+      p[0] = std::sin(phi) * std::cos(theta) * radius;
+      p[1] = std::sin(phi) * std::sin(theta) * radius;
+      p[2] = std::cos(phi) * radius;
    }
 
    return pointCloud;
@@ -177,7 +179,7 @@ int main( int argc, char ** argv )
    //for( auto & p: pointCloud)
    //   p = p.getNormalized() * radius;
 
-   boost::random::mt19937 rng(42);
+   std::mt19937 rng(42);
    for (auto blkIt = forest->begin(); blkIt != forest->end(); ++blkIt)
    {
       IBlock & currentBlock = *blkIt;
diff --git a/tests/mesh/QHullTest.cpp b/tests/mesh/QHullTest.cpp
index 1183325e0..c755a46ec 100644
--- a/tests/mesh/QHullTest.cpp
+++ b/tests/mesh/QHullTest.cpp
@@ -23,6 +23,7 @@
 #include "core/logging/Logging.h"
 #include "core/mpi/Environment.h"
 #include "core/timing/Timer.h"
+#include "core/math/Utility.h"
 
 #include "mesh/TriangleMeshes.h"
 #include "mesh/QHull.h"
@@ -32,7 +33,7 @@
 #include "vtk/VTKOutput.h"
 #include "vtk/PointDataSource.h"
 
-#include <boost/random.hpp>
+#include <random>
 
 #include <vector>
 #include <string>
@@ -190,7 +191,7 @@ std::vector<Vector3<real_t>> generatePointCloudDodecahedron()
 
 std::vector<Vector3<real_t>> generatePointCloudInAABB( const math::AABB & aabb, const uint_t numPoints )
 {
-   boost::random::mt19937 rng(42);
+   std::mt19937 rng(42);
    
    std::vector<Vector3<real_t>> pointCloud( numPoints );
    for( auto & p : pointCloud )
@@ -202,16 +203,17 @@ std::vector<Vector3<real_t>> generatePointCloudInAABB( const math::AABB & aabb,
 
 std::vector<Vector3<real_t>> generatPointCloudOnSphere( const real_t radius, const uint_t numPoints )
 {
-   boost::random::mt19937 rng(42);
-   boost::uniform_on_sphere<real_t> distribution(3);
+   std::mt19937 rng(42);
+   std::uniform_real_distribution<real_t> distribution;
 
    std::vector<Vector3<real_t>> pointCloud( numPoints );
    for( auto & p : pointCloud )
    {
-      auto v = distribution(rng);
-      p[0] = v[0] * radius;
-      p[1] = v[1] * radius;
-      p[2] = v[2] * radius;
+      real_t theta = 2 * math::PI * distribution(rng);
+      real_t phi = std::acos( real_t(1.0) - real_t(2.0) * distribution(rng) );
+      p[0] = std::sin(phi) * std::cos(theta) * radius;
+      p[1] = std::sin(phi) * std::sin(theta) * radius;
+      p[2] = std::cos(phi) * radius;
    }
 
    return pointCloud;
diff --git a/tests/pe/SyncEquivalence.cpp b/tests/pe/SyncEquivalence.cpp
index 9d5cb83b5..c311e1d67 100644
--- a/tests/pe/SyncEquivalence.cpp
+++ b/tests/pe/SyncEquivalence.cpp
@@ -95,9 +95,9 @@ void createSimulation(math::AABB& simulationDomain,
                       real_t vMax,
                       SimInfo& info)
 {
-    boost::mt19937 generator;
+    std::mt19937 generator;
     generator.seed(1337);
-    //math::seedRandomGenerator( 1337 ); //static_cast<boost::mt19937::result_type>(1337 * mpi::MPIManager::instance()->worldRank()) );
+    //math::seedRandomGenerator( 1337 ); //static_cast<std::mt19937::result_type>(1337 * mpi::MPIManager::instance()->worldRank()) );
 
     info.globalBodyStorage = make_shared<BodyStorage>();
 
-- 
GitLab