diff --git a/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp b/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
index 8a18cbccdbc2b9ae94280468cfd60349d85eff94..ce7a492f76cb6aaef9965f7f11286cab878ee134 100644
--- a/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
+++ b/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
@@ -78,8 +78,6 @@
 
 #include "core/timing/RemainingTimeLogger.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <cmath>
 #include <vector>
 #include <string>
diff --git a/apps/benchmarks/DEM/DEM.cpp b/apps/benchmarks/DEM/DEM.cpp
index 52b6e856a4fb47403d3fb0e9db448c3b1756e460..fb9de46118aea30d7cebd0aa9455ed9b2b46b4e6 100644
--- a/apps/benchmarks/DEM/DEM.cpp
+++ b/apps/benchmarks/DEM/DEM.cpp
@@ -32,13 +32,13 @@ namespace dem {
 real_t calcCoefficientOfRestitution(const real_t k, const real_t gamma, const real_t meff)
 {
    auto a = real_t(0.5) * gamma / meff;
-   return std::exp(-a * math::PI / std::sqrt(k / meff - a*a));
+   return std::exp(-a * math::M_PI / std::sqrt(k / meff - a*a));
 }
 
 real_t calcCollisionTime(const real_t k, const real_t gamma, const real_t meff)
 {
    auto a = real_t(0.5) * gamma / meff;
-   return math::PI / std::sqrt( k/meff - a*a);
+   return math::M_PI / std::sqrt( k/meff - a*a);
 }
 }
 
diff --git a/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp b/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
index 0e95853463bb92a51bff94e4606427416dc7d994..ec2212242422eade37e6a5c22171e6852b610e83 100644
--- a/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
+++ b/apps/benchmarks/MeshDistance/MeshDistanceBenchmark.cpp
@@ -164,8 +164,8 @@ int main( int argc, char * argv[] )
       WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args[0] << " [--no-brute-force] [--force-float] MESH_FILE NUM_POINTS NUM_REPETITIONS" );
 
    const std::string & meshFile = args[1];
-   const uint_t numPoints       = boost::lexical_cast<uint_t>( args[2] );
-   const uint_t numRepetitions  = boost::lexical_cast<uint_t>( args[3] );
+   const uint_t numPoints       = string_to_num<uint_t>( args[2] );
+   const uint_t numRepetitions  = string_to_num<uint_t>( args[3] );
 
    if(forceFloat)
    {
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index 2e6f5af9f3ad4f510ad7d3088f6ae46b10318dd8..fb355ddde53c979f3d608f69b524e087df878963 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -1081,7 +1081,7 @@ int main( int argc, char **argv )
       WALBERLA_LOG_INFO_ON_ROOT("Initial simulation has ended.")
 
       //evaluate the gravitational force necessary to keep the sphere at a approximately fixed position
-      gravity = forceEval->getForce() / ( (densityRatio - real_t(1) ) * diameter * diameter * diameter * math::PI / real_t(6) );
+      gravity = forceEval->getForce() / ( (densityRatio - real_t(1) ) * diameter * diameter * diameter * math::M_PI / real_t(6) );
       GalileoSim = std::sqrt( ( densityRatio - real_t(1) ) * gravity * diameter * diameter * diameter ) / viscosity;
       ReynoldsSim = uIn * diameter / viscosity;
       u_ref = std::sqrt( std::fabs(densityRatio - real_t(1)) * gravity * diameter );
@@ -1235,7 +1235,7 @@ int main( int argc, char **argv )
    }
 
    // add gravity
-   Vector3<real_t> extForcesOnSphere( real_t(0), real_t(0), - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::PI / real_t(6));
+   Vector3<real_t> extForcesOnSphere( real_t(0), real_t(0), - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::M_PI / real_t(6));
    timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, extForcesOnSphere ), "Add external forces (gravity and buoyancy)" );
 
    // evaluate the sphere properties
diff --git a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
index dfb8eefd3ab839cf40924e6f709818bae2f0ea68..bd9debc51ed74524a914302789bfcb262299d790 100644
--- a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
+++ b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
@@ -34,6 +34,7 @@
 #include "core/debug/Debug.h"
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "core/math/Sample.h"
 #include "core/math/Vector3.h"
 #include "core/mpi/Environment.h"
@@ -762,7 +763,7 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
    {
       setup.maxVelocity_L = ( setup.acceleration_L * setup.radius_L * setup.radius_L ) / ( real_t(4) * setup.viscosity_L );
       setup.meanVelocity_L = ( setup.acceleration_L * setup.radius_L * setup.radius_L ) / ( real_t(8) * setup.viscosity_L );
-      setup.flowRate_L = setup.meanVelocity_L * math::PI * setup.radius_L * setup.radius_L;
+      setup.flowRate_L = setup.meanVelocity_L * math::M_PI * setup.radius_L * setup.radius_L;
    }
    else
    {
diff --git a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
index d83409c4c6e4b49ebdd745904c7b74fce00fdda2..42885f9587e0dce2ed2f2e27e334ac312d43da7e 100644
--- a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
+++ b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
@@ -38,6 +38,7 @@
 #include "core/debug/Debug.h"
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "core/math/Sample.h"
 #include "core/math/Vector3.h"
 #include "core/mpi/Environment.h"
@@ -746,8 +747,7 @@ public:
 
    void operator()( const real_t t )
    {
-      const real_t PI = real_t( 3.141592653589793238462643383279502884 );
-      tConstTerm_ = ( sinPeriod_ > real_t(0) ) ? ( std::abs( std::sin( PI * t / sinPeriod_ ) ) ) : real_t(1);
+      tConstTerm_ = ( sinPeriod_ > real_t(0) ) ? ( std::abs( std::sin( math::M_PI * t / sinPeriod_ ) ) ) : real_t(1);
       tConstTerm_ *= uTerm_ * HTerm_;
       tConstTerm_ *= ( raisingTime_ > real_t(0) ) ? std::min( t / raisingTime_, real_t(1) ) : real_t(1);
    }
diff --git a/apps/tutorials/pde/01_SolvingPDE.cpp b/apps/tutorials/pde/01_SolvingPDE.cpp
index 55eb90fb22a099a27fefaf27d82c57307f76d0c3..11b47c0f0c989122816afa4060f7d69aab306856 100644
--- a/apps/tutorials/pde/01_SolvingPDE.cpp
+++ b/apps/tutorials/pde/01_SolvingPDE.cpp
@@ -22,6 +22,7 @@
 #include "blockforest/communication/UniformBufferedScheme.h"
 
 #include "core/Environment.h"
+#include "core/math/Constants.h"
 
 #include "field/Field.h"
 #include "field/AddToStorage.h"
@@ -69,9 +70,9 @@ void initBC( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDat
             // obtain the physical coordinate of the center of the current cell
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-            // set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*PI*x)*sinh(2*PI) in the source and destination field
-            src->get( *cell ) = std::sin( real_c(2) * math::PI * p[0] ) * std::sinh( real_c(2) * math::PI );
-            dst->get( *cell ) = std::sin( real_c(2) * math::PI * p[0] ) * std::sinh( real_c(2) * math::PI );
+            // set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*M_PI*x)*sinh(2*M_PI) in the source and destination field
+            src->get( *cell ) = std::sin( real_c(2) * math::M_PI * p[0] ) * std::sinh( real_c(2) * math::M_PI );
+            dst->get( *cell ) = std::sin( real_c(2) * math::M_PI * p[0] ) * std::sinh( real_c(2) * math::M_PI );
          }
       }
    }
@@ -97,9 +98,9 @@ void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDa
          // obtain the physical coordinate of the center of the current cell
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-         // set the right-hand side, given by the function f(x,y) = 4*PI*PI*sin(2*PI*x)*sinh(2*PI*y)
-         f->get( *cell ) = real_c(4) * math::PI * math::PI * std::sin( real_c(2) * math::PI * p[0] ) *
-                           std::sinh( real_c(2) * math::PI * p[1] );
+         // set the right-hand side, given by the function f(x,y) = 4*M_PI*PI*sin(2*M_PI*x)*sinh(2*M_PI*y)
+         f->get( *cell ) = real_c(4) * math::M_PI * math::M_PI * std::sin( real_c(2) * math::M_PI * p[0] ) *
+                           std::sinh( real_c(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -146,7 +147,7 @@ void JacobiSweep::operator()( IBlock * const block )
       dst->get(x,y,z) += ( real_c(1) / (dx_ * dx_) ) * src->get( x-1,  y , z );
       dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y+1, z );
       dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y-1, z );
-      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::PI * math::PI );
+      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::M_PI * math::M_PI );
    )
 
    // swap source and destination fields
@@ -281,7 +282,7 @@ int main( int argc, char ** argv )
    // ...or the variant using the stencil concept
    // set up the stencil weights
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::E ] ] = real_c(-1) / ( dx * dx );
    weights[ Stencil_T::idx[ stencil::W ] ] = real_c(-1) / ( dx * dx );
    weights[ Stencil_T::idx[ stencil::N ] ] = real_c(-1) / ( dy * dy );
diff --git a/apps/tutorials/pde/01_SolvingPDE.dox b/apps/tutorials/pde/01_SolvingPDE.dox
index 9298f9c687f462a6a2a96d89bffccdfb20f7b621..8a7e33dce7da53943f4fbdffd4a5f387bfd1165f 100644
--- a/apps/tutorials/pde/01_SolvingPDE.dox
+++ b/apps/tutorials/pde/01_SolvingPDE.dox
@@ -71,7 +71,7 @@ void JacobiSweep::operator()( IBlock * const block )
       dst->get(x,y,z) += ( real_c(1) / (dx_ * dx_) ) * src->get( x-1,  y , z );
       dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y+1, z );
       dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y-1, z );
-      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::PI * math::PI );
+      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::M_PI * math::M_PI );
    )
 
    // swap source and destination fields
@@ -90,7 +90,7 @@ typedef stencil::D2Q5 Stencil_T;
 
 // set up the stencil weights
 std::vector< real_t > weights( Stencil_T::Size );
-weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::PI * math::PI;
+weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::M_PI * math::M_PI;
 weights[ Stencil_T::idx[ stencil::E ] ] = real_c(-1) / ( dx * dx );
 weights[ Stencil_T::idx[ stencil::W ] ] = real_c(-1) / ( dx * dx );
 weights[ Stencil_T::idx[ stencil::N ] ] = real_c(-1) / ( dy * dy );
@@ -185,9 +185,9 @@ void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDa
          // obtain the physical coordinate of the center of the current cell
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-         // set the right-hand side, given by the function f(x,y) = 4*PI*PI*sin(2*PI*x)*sinh(2*PI*y)
-         f->get( *cell ) = real_c(4) * math::PI * math::PI * std::sin( real_c(2) * math::PI * p[0] ) *
-                           std::sinh( real_c(2) * math::PI * p[1] );
+         // set the right-hand side, given by the function f(x,y) = 4*M_PI*PI*sin(2*M_PI*x)*sinh(2*M_PI*y)
+         f->get( *cell ) = real_c(4) * math::M_PI * math::M_PI * std::sin( real_c(2) * math::M_PI * p[0] ) *
+                           std::sinh( real_c(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -236,10 +236,10 @@ void initBC( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDat
             // obtain the physical coordinate of the center of the current cell
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-            // set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*PI*x)*sinh(2*PI)
+            // set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*M_PI*x)*sinh(2*M_PI)
             // in the source and destination field
-            src->get( *cell ) = std::sin( real_c(2) * math::PI * p[0] ) * std::sinh( real_c(2) * math::PI );
-            dst->get( *cell ) = std::sin( real_c(2) * math::PI * p[0] ) * std::sinh( real_c(2) * math::PI );
+            src->get( *cell ) = std::sin( real_c(2) * math::M_PI * p[0] ) * std::sinh( real_c(2) * math::M_PI );
+            dst->get( *cell ) = std::sin( real_c(2) * math::M_PI * p[0] ) * std::sinh( real_c(2) * math::M_PI );
          }
       }
    }
diff --git a/apps/tutorials/pde/02_HeatEquation.cpp b/apps/tutorials/pde/02_HeatEquation.cpp
index 931644c344b2770deec2a63ce653c7cc82c6d481..c4ec151ea024d15505e13429e64c82e642111700 100644
--- a/apps/tutorials/pde/02_HeatEquation.cpp
+++ b/apps/tutorials/pde/02_HeatEquation.cpp
@@ -22,6 +22,7 @@
 #include "blockforest/communication/UniformBufferedScheme.h"
 
 #include "core/Environment.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Reduce.h"
 
 #include "field/Field.h"
@@ -65,8 +66,8 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          // obtain the physical coordinate of the center of the current cell
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-         // set the initial condition, given by the function u(x,y,0) = sin(PI*x)*sin(PI*y)
-         u->get( *cell ) = std::sin( math::PI * p[0] ) * std::sin( math::PI * p[1] );
+         // set the initial condition, given by the function u(x,y,0) = sin(M_PI*x)*sin(M_PI*y)
+         u->get( *cell ) = std::sin( math::M_PI * p[0] ) * std::sin( math::M_PI * p[1] );
       }
    }
 }
@@ -255,4 +256,4 @@ int main( int argc, char ** argv )
 int main( int argc, char ** argv )
 {
    walberla::main(argc, argv);
-}
\ No newline at end of file
+}
diff --git a/apps/tutorials/pde/03_HeatEquation_Extensions.cpp b/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
index 15d8859c1997fb064f0f9715ea7171f1b248fed8..92f3b79acf19310c2baa6e5c192ad1409dfc3248 100644
--- a/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
+++ b/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
@@ -22,6 +22,7 @@
 #include "blockforest/communication/UniformBufferedScheme.h"
 
 #include "core/Environment.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Reduce.h"
 
 #include "field/Field.h"
@@ -66,8 +67,8 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          // obtain the physical coordinate of the center of the current cell
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
 
-         // set the initial condition, given by the function u(x,y,0) = sin(PI*x)*sin(PI*y)
-         u->get( *cell ) = std::sin( math::PI * p[0] ) * std::sin( math::PI * p[1] );
+         // set the initial condition, given by the function u(x,y,0) = sin(M_PI*x)*sin(M_PI*y)
+         u->get( *cell ) = std::sin( math::M_PI * p[0] ) * std::sin( math::M_PI * p[1] );
       }
    }
 }
@@ -352,4 +353,4 @@ int main( int argc, char ** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/src/core/DataTypes.h b/src/core/DataTypes.h
index e4c755d104988335a71158aaf2f9955e12d523a9..fb7d16af20a19eaad7a98baa8bdf1971f48c5d3b 100644
--- a/src/core/DataTypes.h
+++ b/src/core/DataTypes.h
@@ -24,14 +24,17 @@
 #include "waLBerlaDefinitions.h"
 
 #include <cstdint>
-#include <boost/numeric/conversion/cast.hpp>
 #include <memory>
 #include <type_traits>
-#include <boost/units/detail/utility.hpp>
-
+#include <string>
 #include <cmath>
 #include <limits>
 
+#if (defined( __has_include ) && __has_include(<cxxabi.h>)) || defined( __GLIBCXX__ )
+#define HAVE_CXXABI_H
+#include <cxxabi.h>
+#endif
+
 
 namespace walberla {
 
@@ -41,6 +44,8 @@ namespace walberla {
 
 template <typename> struct never_true : std::false_type {};
 
+template< typename T > bool isIdentical( const T a, const T b );
+
 
 // shared ptr
 
@@ -54,15 +59,52 @@ using std::dynamic_pointer_cast;
 
 template< typename S, typename T >
 inline S numeric_cast( T t ) {
-#ifdef NDEBUG
-   return static_cast< S >(t);
-#else
-   return boost::numeric_cast< S >(t);
+#ifndef NDEBUG
+   if( std::is_integral<S>::value && std::is_integral<T>::value && !std::is_same<S,T>::value )
+        // integer to different integer: check that forward and back conversion does not change value
+   {
+      if( !isIdentical( static_cast<T>( static_cast<S>(t) ), t ) )
+      {
+         throw std::range_error("out of range");
+      }
+   }
+   else if( !std::is_integral<S>::value && !std::is_integral<T>::value && sizeof(S) < sizeof(T) )
+       // float to shorter float: check that value within limits of shorter type
+   {
+      using H = typename std::conditional< !std::is_integral<S>::value && !std::is_integral<T>::value && (sizeof(S) < sizeof(T)), T, long double >::type; // always true, but makes Intel's overflow check happy
+      H h = static_cast<H>(t);
+      if( h < static_cast<H>(std::numeric_limits<S>::lowest()) || h > static_cast<H>(std::numeric_limits<S>::max()) ) {
+         throw std::range_error("out of range");
+      }
+   }
+   else if( std::is_integral<S>::value && !std::is_integral<T>::value )
+       // float to integer: check that value within limits of integer
+   {
+      using H = typename std::conditional< std::is_integral<S>::value && !std::is_integral<T>::value, T, long double >::type; // always true, but makes Intel's overflow check happy
+      H h = static_cast<H>(t);
+      if( h < static_cast<H>(std::numeric_limits<S>::lowest()) || h > static_cast<H>(std::numeric_limits<S>::max()) ) {
+         throw std::range_error("out of range");
+      }
+   }
 #endif
+   return static_cast< S >(t);
 }
 
 
 
+template<typename S>
+inline S string_to_num( std::string & t );
+template <> inline float              string_to_num( std::string & t ) { return std::stof(t); }
+template <> inline double             string_to_num( std::string & t ) { return std::stod(t); }
+template <> inline long double        string_to_num( std::string & t ) { return std::stold(t); }
+template <> inline int                string_to_num( std::string & t ) { return std::stoi(t); }
+template <> inline long               string_to_num( std::string & t ) { return std::stol(t); }
+template <> inline long long          string_to_num( std::string & t ) { return std::stoll(t); }
+template <> inline unsigned long      string_to_num( std::string & t ) { return std::stoul(t); }
+template <> inline unsigned long long string_to_num( std::string & t ) { return std::stoull(t); }
+
+
+
 // fixed size signed integral types
 typedef std::int8_t   int8_t;    ///<  8 bit signed integer
 typedef std::int16_t  int16_t;   ///< 16 bit signed integer
@@ -286,7 +328,20 @@ template< typename T > inline const char* typeToString( T ) { return typeToStrin
 
 inline std::string demangle( const std::string & name )
 {
-   return boost::units::detail::demangle( name.c_str() );
+#ifdef HAVE_CXXABI_H
+   int status = 0;
+   std::size_t size = 0;
+   const char * demangled = abi::__cxa_demangle( name.c_str(), NULL, &size, &status );
+   if( demangled == nullptr )
+   {
+      return name;
+   }
+   std::string demangled_str(demangled);
+   std::free( const_cast<char*>(demangled) );
+   return demangled_str;
+#else
+   return name;
+#endif
 }
 
 
diff --git a/src/core/cell/Cell.h b/src/core/cell/Cell.h
index cf54ccaae269df7c6c4114c0eb5ed60941a2ccdf..11fa7d1d1470084c6a6701bb3e9011b771270b22 100644
--- a/src/core/cell/Cell.h
+++ b/src/core/cell/Cell.h
@@ -28,10 +28,9 @@
 #include "core/mpi/RecvBuffer.h"
 #include "core/mpi/SendBuffer.h"
 
-#include <boost/functional/hash.hpp>
 #include <algorithm>
 #include <iterator>
-#include <ostream>
+#include <iostream>
 
 
 namespace walberla {
@@ -336,9 +335,12 @@ inline std::istream & operator>>( std::istream & is, Cell & cell )
 inline std::size_t hash_value( const Cell & cell )
 {
   std::size_t seed = 0;
-  boost::hash_combine( seed, cell.x() );
-  boost::hash_combine( seed, cell.y() );
-  boost::hash_combine( seed, cell.z() );
+  std::hash<cell_idx_t> hasher;
+
+  seed ^= hasher(cell.x()) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+  seed ^= hasher(cell.y()) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+  seed ^= hasher(cell.z()) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+
   return seed;
 }
 
diff --git a/src/core/debug/CheckFunctions.impl.h b/src/core/debug/CheckFunctions.impl.h
index fac44221d516a0a8be73690787eecff21a672e88..71f0895654f7c16d2e3f49a0b76de941c196679a 100644
--- a/src/core/debug/CheckFunctions.impl.h
+++ b/src/core/debug/CheckFunctions.impl.h
@@ -20,6 +20,8 @@
 //======================================================================================================================
 
 
+#include <algorithm>
+
 #include <boost/type_traits/integral_constant.hpp>
 
 /// \cond internal
diff --git a/src/core/debug/PrintStacktrace.cpp b/src/core/debug/PrintStacktrace.cpp
index b4a8e6888c86784861f9985a358345011a5c300b..a83b13d23ad7c33e0e2b4b6d6d3b5caa075bc175 100644
--- a/src/core/debug/PrintStacktrace.cpp
+++ b/src/core/debug/PrintStacktrace.cpp
@@ -65,6 +65,24 @@ namespace debug {
       for (i = 0; i < size; i++)
       {
          std::string line ( strings[i] );
+#ifdef __APPLE__
+         // one line might look like this:
+         //0   PrintStacktraceTest                 0x0000000000408c6b _ZN8walberla4core15printStacktraceEv + 75
+
+         size_t plusPos       = line.find_last_of('+');
+         size_t funcPos       = line.find_last_of(' ', plusPos-2)+1;
+
+         string functionName = line.substr( funcPos, plusPos-funcPos-1 );
+
+         size_t addrPos       = line.find_last_of('x', funcPos-2)-2;
+         size_t afterLevelPos = line.find_first_of(' ');
+         size_t beforeAppPos  = line.find_first_not_of(" ", afterLevelPos);
+
+         string appName      = line.substr( beforeAppPos, addrPos-beforeAppPos );
+
+         size_t afterAppPos   = appName.find_last_not_of(" ")+1;
+         appName             = appName.substr( 0, afterAppPos );
+#else
          // one line might look like this:
          //./PrintStacktraceTest(_ZN8walberla4core15printStacktraceEv+0x4b) [0x408c6b]
 
@@ -80,6 +98,7 @@ namespace debug {
          size_t plusPos = bracketPart.find_first_of('+');
          string functionName = bracketPart.substr(0, plusPos );
          string offset       = bracketPart.substr( plusPos+1 );
+#endif
 
          // try to demangle -> no return code if successfull
          // but the returned string starts with "demangle" if demangling failed
diff --git a/src/core/math/DistributedSample.cpp b/src/core/math/DistributedSample.cpp
index 509af72af40ccd817e29f183e0d1b8925574a414..3a1d40a7f191fc304bfd567a1532227248057e84 100644
--- a/src/core/math/DistributedSample.cpp
+++ b/src/core/math/DistributedSample.cpp
@@ -25,7 +25,6 @@
 #include "core/mpi/Reduce.h"
 
 #include <boost/algorithm/string/replace.hpp>
-#include <boost/lexical_cast.hpp>
 
 
 
@@ -172,28 +171,27 @@ void DistributedSample::mpiGatherRoot()
 std::string DistributedSample::format( const std::string & formatString ) const
 {
    using boost::algorithm::replace_all;
-   using boost::lexical_cast;
 
    std::string result = formatString;
 
    if( size_ > uint_t(0) )
    {
       if( result.find( "%min" ) != std::string::npos )
-         replace_all(result, "%min", lexical_cast<std::string>( min_ ) );
+         replace_all(result, "%min", std::to_string( min_ ) );
       if( result.find( "%max" ) != std::string::npos )
-         replace_all(result, "%max", lexical_cast<std::string>( max_ ) );
+         replace_all(result, "%max", std::to_string( max_ ) );
       if( result.find( "%sum" ) != std::string::npos )
-         replace_all(result, "%sum",lexical_cast<std::string>( sum_ ) );
+         replace_all(result, "%sum",std::to_string( sum_ ) );
       if( result.find( "%mean" ) != std::string::npos )
-         replace_all(result, "%mean", lexical_cast<std::string>( mean_ ) );
+         replace_all(result, "%mean", std::to_string( mean_ ) );
       if( result.find( "%var" ) != std::string::npos )
-         replace_all(result, "%var", lexical_cast<std::string>( variance_ ) );
+         replace_all(result, "%var", std::to_string( variance_ ) );
       if( result.find( "%stddev" ) != std::string::npos )
-         replace_all(result, "%stddev", lexical_cast<std::string>( stdDeviation() ) );
+         replace_all(result, "%stddev", std::to_string( stdDeviation() ) );
       if( result.find( "%relstddev" ) != std::string::npos )
-         replace_all(result, "%relstddev", lexical_cast<std::string>( relativeStdDeviation() ) );
+         replace_all(result, "%relstddev", std::to_string( relativeStdDeviation() ) );
       if( result.find( "%size" ) != std::string::npos )
-         replace_all(result, "%size", lexical_cast<std::string>( size_ ) );
+         replace_all(result, "%size", std::to_string( size_ ) );
    }
    else // empty()
    {
diff --git a/src/core/math/PhysicalCheck.cpp b/src/core/math/PhysicalCheck.cpp
index d87a0ad18b8abd855a340a03e748c1743e8c1ae3..7c3342449e13ab888e8186c3e385f204b25ffa7b 100644
--- a/src/core/math/PhysicalCheck.cpp
+++ b/src/core/math/PhysicalCheck.cpp
@@ -25,8 +25,6 @@
 #include "core/math/Parser.h"
 #include "equation_system/EquationParser.h"
 
-#include <boost/lexical_cast.hpp>
-
 
 namespace walberla {
 namespace math {
@@ -96,7 +94,7 @@ namespace math {
 
       for (size_t i=0; i<equations.size(); ++i){
          index = 0;
-         es_.add( boost::lexical_cast<std::string>(es_.getNumberOfEquations()+1), ep.parseEquation( equations[i], index ) );
+         es_.add( std::to_string(es_.getNumberOfEquations()+1), ep.parseEquation( equations[i], index ) );
       }
    }
 
@@ -107,7 +105,7 @@ namespace math {
 
       solved_ = false;
 
-      es_.add( boost::lexical_cast<std::string>(es_.getNumberOfEquations()+1), ep.parseEquation( equation, index ) );
+      es_.add( std::to_string(es_.getNumberOfEquations()+1), ep.parseEquation( equation, index ) );
    }
 
    void PhysicalCheck::addUnitParameterRelations( const std::map< std::string, std::string >& unitParameterRelations )
@@ -134,7 +132,7 @@ namespace math {
                if( j < unitString.size() && unitString[j] == '^' )
                {
                   ++j;
-                  int expo = factor * boost::lexical_cast<int>( unitString[j] );
+                  int expo = factor * std::atoi( &unitString[j] );
                   if( !setVarUnit( i->first, unitString.substr(index,1), expo ) )
                      WALBERLA_ABORT( "Error in PhysicalCheck::addUnitParameterRelations(). Non-unique description for unit '" << unitString[index] << "' for parameter '" << i->first << "'." );
                }
@@ -153,7 +151,7 @@ namespace math {
                if( j < unitString.size() && unitString[j] == '^' )
                {
                   ++j;
-                  int expo = factor * boost::lexical_cast<int>( unitString[j] );
+                  int expo = factor * std::atoi( &unitString[j] );
                   if( !setVarUnit( i->first, unitString.substr(index,2), expo ) )
                      WALBERLA_ABORT( "Error in PhysicalCheck::addUnitParameterRelations(). Non-unique description for unit 'kg' for parameter '" << i->first << "'.");
                }
@@ -277,7 +275,9 @@ namespace math {
             }
 
             // set the current parameter to the evaluated result in the current block
-            configBlock.setParameter( param->first, boost::lexical_cast<std::string>(result) );
+            std::ostringstream os;
+            os << std::setprecision(16) << result;
+            configBlock.setParameter( param->first, os.str() );
 
             // check for the equality of the parameter and the result
             double val = configBlock.getParameter<double>( param->first );
diff --git a/src/core/math/Sample.cpp b/src/core/math/Sample.cpp
index eabe897cd81c45ba07445f168e87841eecea99fa..8e376edf74f71b5f8bdfb68f97ea8515a33c34c7 100644
--- a/src/core/math/Sample.cpp
+++ b/src/core/math/Sample.cpp
@@ -25,7 +25,6 @@
 #include "core/mpi/MPIManager.h"
 
 #include <boost/algorithm/string/replace.hpp>
-#include <boost/lexical_cast.hpp>
 
 #include <functional>
 #include <iterator>
@@ -271,32 +270,31 @@ real_t Sample::quantile(const real_t p) const
 std::string Sample::format(const std::string & formatString) const
 {
    using boost::algorithm::replace_all;
-   using boost::lexical_cast;
 
    std::string result = formatString;
 
    if( !empty() )
    {
       if( result.find( "%min" ) != std::string::npos )
-         replace_all(result, "%min", lexical_cast<std::string>( min() ) );
+         replace_all(result, "%min", std::to_string( min() ) );
       if( result.find( "%max" ) != std::string::npos )
-         replace_all(result, "%max", lexical_cast<std::string>( max() ) );
+         replace_all(result, "%max", std::to_string( max() ) );
       if( result.find( "%sum" ) != std::string::npos )
-         replace_all(result, "%sum",lexical_cast<std::string>( sum() ) );
+         replace_all(result, "%sum",std::to_string( sum() ) );
       if( result.find( "%mean" ) != std::string::npos )
-         replace_all(result, "%mean", lexical_cast<std::string>( mean() ) );
+         replace_all(result, "%mean", std::to_string( mean() ) );
       if( result.find( "%med" ) != std::string::npos )
-         replace_all(result, "%med", lexical_cast<std::string>( median() ) );
+         replace_all(result, "%med", std::to_string( median() ) );
       if( result.find( "%var" ) != std::string::npos )
-         replace_all(result, "%var", lexical_cast<std::string>( variance() ) );
+         replace_all(result, "%var", std::to_string( variance() ) );
       if( result.find( "%stddev" ) != std::string::npos )
-         replace_all(result, "%stddev", lexical_cast<std::string>( stdDeviation() ) );
+         replace_all(result, "%stddev", std::to_string( stdDeviation() ) );
       if( result.find( "%relstddev" ) != std::string::npos )
-         replace_all(result, "%relstddev", lexical_cast<std::string>( relativeStdDeviation() ) );
+         replace_all(result, "%relstddev", std::to_string( relativeStdDeviation() ) );
       if( result.find( "%mad" ) != std::string::npos )
-         replace_all(result, "%mad", lexical_cast<std::string>( mad() ) );
+         replace_all(result, "%mad", std::to_string( mad() ) );
       if( result.find( "%size" ) != std::string::npos )
-         replace_all(result, "%size", lexical_cast<std::string>( size() ) );
+         replace_all(result, "%size", std::to_string( size() ) );
    }
    else // empty()
    {
diff --git a/src/core/math/Utility.h b/src/core/math/Utility.h
index 198cd4c73dd9e311e3f8fa2766d2a827f9eb0dd2..ba1f62c108f555806fe92f9d62effe2187cf6e16 100644
--- a/src/core/math/Utility.h
+++ b/src/core/math/Utility.h
@@ -25,26 +25,15 @@
 #include "MathTrait.h"
 #include "core/DataTypes.h"
 
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <cstddef>
 #include <limits>
+#include <type_traits>
 
 
 namespace walberla {
 namespace math {
 
-//======================================================================================================================
-//
-//  MATHEMATICAL CONSTANTS
-//
-//======================================================================================================================
-
-//! Definition of the mathematical value \f$ \pi \f$.
-const real_t PI = boost::math::constants::pi<real_t>();
-
-
-
 //======================================================================================================================
 //
 //  MATHEMATICAL UTILITY FUNCTIONS
@@ -58,9 +47,9 @@ template< typename T >
 inline const T sign( T a );
 
 template< typename T >
-inline const typename boost::enable_if_c<  boost::is_unsigned<T>::value, T >::type abs( T a );
+inline const typename std::enable_if<  std::is_unsigned<T>::value, T >::type abs( T a );
 template< typename T >
-inline const typename boost::enable_if_c< !boost::is_unsigned<T>::value, T >::type abs( T a );
+inline const typename std::enable_if< !std::is_unsigned<T>::value, T >::type abs( T a );
 
 
 template< typename T1, typename T2 >
@@ -113,13 +102,13 @@ inline const T sign( T a )
 // \return The value if it is greater than or equal to zero, -1 times the value if the value is smaller than zero.
  */
 template< typename T >
-inline const typename boost::enable_if_c<  boost::is_unsigned<T>::value, T >::type abs( T a )
+inline const typename std::enable_if<  std::is_unsigned<T>::value, T >::type abs( T a )
 {
    return a;
 }
 
 template< typename T >
-inline const typename boost::enable_if_c< !boost::is_unsigned<T>::value, T >::type abs( T a )
+inline const typename std::enable_if< !std::is_unsigned<T>::value, T >::type abs( T a )
 {
    return std::abs( a );
 }
diff --git a/src/core/math/Vector2.h b/src/core/math/Vector2.h
index a6ad5c30091e2de9c14245e62590659bb4a9ec45..220bb21b060b4f81a7d0a48c244b34732bfd72bd 100644
--- a/src/core/math/Vector2.h
+++ b/src/core/math/Vector2.h
@@ -37,7 +37,6 @@
 #include "core/mpi/RecvBuffer.h"
 #include "core/mpi/SendBuffer.h"
 
-#include <boost/functional/hash.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/type_traits/is_fundamental.hpp>
 #include <boost/type_traits/is_same.hpp>
@@ -1579,9 +1578,10 @@ template< typename T >
 size_t hash_value( const Vector2<T> & v )
 {
    size_t seed = 0;
+   std::hash<T> hasher;
 
-   boost::hash_combine( seed, v[0] );
-   boost::hash_combine( seed, v[1] );
+   seed ^= hasher(v[0]) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+   seed ^= hasher(v[1]) + 0x9e3779b9 + (seed<<6) + (seed>>2);
 
    return seed;
 
diff --git a/src/core/math/Vector3.h b/src/core/math/Vector3.h
index 24485fef5a005aaf95dcfdc8f1fe913d62d4ab60..e0c2965795ab2a9c00f06fff7af96d348c7e6ef6 100644
--- a/src/core/math/Vector3.h
+++ b/src/core/math/Vector3.h
@@ -38,7 +38,6 @@
 #include "core/mpi/SendBuffer.h"
 #include "core/debug/CheckFunctions.h"
 
-#include <boost/functional/hash.hpp>
 #include <boost/type_traits.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/type_traits/is_fundamental.hpp>
@@ -1837,10 +1836,11 @@ template< typename T >
 size_t hash_value( const Vector3<T> & v )
 {
    size_t seed = 0;
+   std::hash<T> hasher;
 
-   boost::hash_combine( seed, v[0] );
-   boost::hash_combine( seed, v[1] );
-   boost::hash_combine( seed, v[2] );
+   seed ^= hasher(v[0]) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+   seed ^= hasher(v[1]) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+   seed ^= hasher(v[2]) + 0x9e3779b9 + (seed<<6) + (seed>>2);
 
    return seed;
 
diff --git a/src/core/math/all.h b/src/core/math/all.h
index 162414f7bf491709eb9a96eaa3db3b73b0922914..29295fc4f159d9540e7e21b4e54ccb9270c222ec 100644
--- a/src/core/math/all.h
+++ b/src/core/math/all.h
@@ -23,6 +23,7 @@
 #pragma once
 
 #include "AABB.h"
+#include "Constants.h"
 #include "DistributedSample.h"
 #include "FPClassify.h"
 #include "FastInvSqrt.h"
diff --git a/src/core/math/equation_system/EquationParser.cpp b/src/core/math/equation_system/EquationParser.cpp
index a58be4bd796d1a9eb33c8e8f7f37c9aa1c3e4a3a..cba3f45b1f7cde60118bca90b8723d54f80d92c4 100644
--- a/src/core/math/equation_system/EquationParser.cpp
+++ b/src/core/math/equation_system/EquationParser.cpp
@@ -26,7 +26,6 @@
 #include "core/math/Constants.h"
 
 #include <boost/algorithm/string/trim.hpp>
-#include <boost/lexical_cast.hpp>
 #include <memory>
 
 
@@ -78,10 +77,10 @@ NodePtr EquationParser::parseNumber( const std::string& str, size_t& index ) con
          THROW( "Number ends with 'e'", str, index );
       while( isdigit(str[++index]) != int(0) );
 
-      value =  boost::lexical_cast< double >( str.substr(start, estart-start-1) ) *
-            pow(10.0, boost::lexical_cast< int >( str.substr(estart, index-estart) ) );
+      value =  std::stod( str.substr(start, estart-start-1) ) *
+            pow(10.0, std::stoi( str.substr(estart, index-estart) ) );
    } else {
-      value = boost::lexical_cast< double >( str.substr(start, index-start) );
+      value = std::stod( str.substr(start, index-start) );
    }
 
    return std::make_shared<Node>( value );
diff --git a/src/core/mpi/BufferSystem.h b/src/core/mpi/BufferSystem.h
index 2a04730f7e2cbade9c4bf2e8addc612786945234..6a0048b2c1055992c413fd82f6045fd7edcf1b41 100644
--- a/src/core/mpi/BufferSystem.h
+++ b/src/core/mpi/BufferSystem.h
@@ -26,7 +26,6 @@
 #include "core/DataTypes.h"
 #include "core/debug/Debug.h"
 
-#include <boost/range/counting_range.hpp>
 #include <map>
 #include <set>
 #include <vector>
@@ -209,8 +208,7 @@ public:
    //* Rank Ranges     *************************************************************************************************
    /*! \name Rank Ranges  */
    //@{
-   typedef boost::counting_iterator<MPIRank> RankCountIter;
-   typedef boost::iterator_range< RankCountIter > RankRange;
+   typedef std::set<MPIRank> RankRange;
    static RankRange noRanks();
    static RankRange allRanks();
    static RankRange allRanksButRoot();
diff --git a/src/core/mpi/BufferSystem.impl.h b/src/core/mpi/BufferSystem.impl.h
index 852c7a1fac9a14f490917a4340d4fa180d3525c8..483013860cfb65df8c0892ac770ec21057b89953 100644
--- a/src/core/mpi/BufferSystem.impl.h
+++ b/src/core/mpi/BufferSystem.impl.h
@@ -525,36 +525,35 @@ void GenericBufferSystem<Rb, Sb>::setCommunicationType( const bool knownSize )
 template< typename Rb, typename Sb>
 typename GenericBufferSystem<Rb, Sb>::RankRange GenericBufferSystem<Rb,Sb>::noRanks()
 {
-   return RankRange ( RankCountIter( 0 ),
-                      RankCountIter( 0 ) );
+   return RankRange();
 }
 template< typename Rb, typename Sb>
 typename GenericBufferSystem<Rb, Sb>::RankRange GenericBufferSystem<Rb,Sb>::allRanks()
 {
    int numProcesses = MPIManager::instance()->numProcesses();
-   return RankRange ( RankCountIter( 0            ),
-                      RankCountIter( numProcesses ) );
+   RankRange r;
+   std::generate_n( std::inserter(r, r.end()), numProcesses, [&]{ return MPIRank(r.size()); } );
+   return r;
 }
 template< typename Rb, typename Sb>
 typename GenericBufferSystem<Rb, Sb>::RankRange GenericBufferSystem<Rb,Sb>::allRanksButRoot()
 {
    int numProcesses = MPIManager::instance()->numProcesses();
-   return RankRange ( RankCountIter( 1            ),
-                      RankCountIter( numProcesses ) );
+   RankRange r;
+   std::generate_n( std::inserter(r, r.end()), numProcesses-1, [&]{ return MPIRank(r.size()+size_t(1)); } );
+   return r;
 }
 template< typename Rb, typename Sb>
 typename GenericBufferSystem<Rb, Sb>::RankRange GenericBufferSystem<Rb,Sb>::onlyRank( int includedRank )
 {
    WALBERLA_ASSERT_LESS( includedRank, MPIManager::instance()->numProcesses() );
-   return RankRange ( RankCountIter( includedRank   ),
-                      RankCountIter( includedRank+1 ) );
+   return RankRange { includedRank };
 }
 
 template< typename Rb, typename Sb>
 typename GenericBufferSystem<Rb, Sb>::RankRange GenericBufferSystem<Rb,Sb>::onlyRoot()
 {
-   return RankRange ( RankCountIter( 0 ),
-                      RankCountIter( 1 ) );
+   return RankRange { 0 };
 }
 
 
diff --git a/src/core/mpi/OpenMPBufferSystem.impl.h b/src/core/mpi/OpenMPBufferSystem.impl.h
index c79b961ba0b0deaf0a8af506fc508e80715900cb..dc042e8f44d870ff02d513938006c0ad981579dd 100644
--- a/src/core/mpi/OpenMPBufferSystem.impl.h
+++ b/src/core/mpi/OpenMPBufferSystem.impl.h
@@ -19,8 +19,6 @@
 //
 //======================================================================================================================
 
-#include <boost/range/adaptor/map.hpp>
-
 
 namespace walberla {
 namespace mpi {
@@ -69,8 +67,10 @@ void GenericOpenMPBufferSystem<Rb, Sb>::setupBufferSystem()
    if ( ! dirty_ )
       return;
 
-   using boost::adaptors::map_keys;
-   bs_.setReceiverInfo( recvFunctions_ | map_keys, sizeChangesEverytime_ );
+   std::set<MPIRank> recvRanks;
+   std::transform( recvFunctions_.begin(), recvFunctions_.end(), std::inserter(recvRanks, recvRanks.end()),
+                   []( typename decltype(recvFunctions_)::value_type r ){ return r.first; });
+   bs_.setReceiverInfo( recvRanks, sizeChangesEverytime_ );
 
    for( auto sendRank = sendRanks_.begin(); sendRank != sendRanks_.end(); ++sendRank ) // Do NOT delete this for loop! This loop is needed ...
       bs_.sendBuffer( *sendRank ); // ... so that the "sendBuffer(rank)" call in startCommunicationOpenMP is thread-safe!
diff --git a/src/core/mpi/Reduce.h b/src/core/mpi/Reduce.h
index cb4afcfaab0339d7a2329a1e2ef9deec069289c5..cc1ba91c2fd564b077458249f3740a9a94716850 100644
--- a/src/core/mpi/Reduce.h
+++ b/src/core/mpi/Reduce.h
@@ -30,6 +30,7 @@
 #include "core/mpi/MPIWrapper.h"
 
 #include <boost/type_traits/is_arithmetic.hpp>
+#include <boost/type_traits/is_same.hpp>
 #include <vector>
 
 
diff --git a/src/core/uid/UID.h b/src/core/uid/UID.h
index 7739f7a39fe5139d27f9dd6614111d98d3e2c731..f90f7f9de7e8bf9cb5a65d1df5341e806f73d96e 100644
--- a/src/core/uid/UID.h
+++ b/src/core/uid/UID.h
@@ -27,8 +27,6 @@
 #include "core/debug/Debug.h"
 #include "core/mpi/BufferDataTypeExtensions.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <iostream>
 #include <map>
 #include <sstream>
@@ -201,7 +199,7 @@ void UID<T>::init( const std::string& identifier, const bool newUid, const bool
 
       uid_ = T::generateUID();
 
-      std::string idString = ( appendUIDtoIdentifier ? ( identifier + " [" + boost::lexical_cast< std::string >( uid_ ) + "]" ) : identifier );
+      std::string idString = ( appendUIDtoIdentifier ? ( identifier + " [" + std::to_string( uid_ ) + "]" ) : identifier );
 
       WALBERLA_ASSERT( stringToUid().find( idString ) == stringToUid().end() ); // 'idString' must not exist
 
@@ -250,7 +248,7 @@ UID<T>::UID( const bool createAnonymousUID ) : uid_( 0 ) {
 
       uid_ = T::generateUID();
 
-      std::string identifier = std::string("[anonymous #") + boost::lexical_cast< std::string >( T::toIndex( uid_ ) ) + std::string("]");
+      std::string identifier = std::string("[anonymous #") + std::to_string( T::toIndex( uid_ ) ) + std::string("]");
 
       WALBERLA_ASSERT( stringToUid().find( identifier ) == stringToUid().end() ); // 'identifier' must not exist
 
diff --git a/src/cuda/AddGPUFieldToStorage.h b/src/cuda/AddGPUFieldToStorage.h
index 358a3397cdf7fe635fb30d07a2d24e3d27147695..3803ff8c814bb4df9e21312cbd00bee562100492 100644
--- a/src/cuda/AddGPUFieldToStorage.h
+++ b/src/cuda/AddGPUFieldToStorage.h
@@ -25,8 +25,6 @@
 
 #include "domain_decomposition/StructuredBlockStorage.h"
 
-#include <boost/bind.hpp>
-
 
 namespace walberla {
 namespace cuda {
diff --git a/src/gather/FileGatherScheme.cpp b/src/gather/FileGatherScheme.cpp
index 0aa2e40d7121db5eacbda388a5f55afcabe434c9..238fe341034895fbce46a485381762805c4ed5ae 100644
--- a/src/gather/FileGatherScheme.cpp
+++ b/src/gather/FileGatherScheme.cpp
@@ -28,7 +28,6 @@
 #include "core/mpi/SendBuffer.h"
 
 #include "core/Filesystem.h"
-#include <boost/lexical_cast.hpp>
 
 
 namespace walberla {
@@ -57,12 +56,11 @@ FileGatherScheme::~FileGatherScheme()
 void FileGatherScheme::deleteTemporaryFiles()
 {
    using namespace filesystem;
-   using std::string;
 
    for(int rank=0; rank < MPIManager::instance()->numProcesses(); ++rank)
    {
-      string fileName = "tmp_collection_" + boost::lexical_cast<string>(myId_) + "_"
-                                          + boost::lexical_cast<string>(rank);
+      std::string fileName = "tmp_collection_" + std::to_string(myId_) + "_"
+                                               + std::to_string(rank);
 
       if( exists(fileName) )
          remove(fileName);
@@ -73,8 +71,6 @@ void FileGatherScheme::deleteTemporaryFiles()
 
 void FileGatherScheme::writeToFile()
 {
-   using std::string;
-
    static mpi::SendBuffer buffer;
 
    buffer.clear();
@@ -88,8 +84,8 @@ void FileGatherScheme::writeToFile()
       // not important if rank() or worldRank() we need only a distinct number
       // for each process
       int rank = MPIManager::instance()->worldRank();
-      string fileName = "tmp_collection_" + boost::lexical_cast<string>(myId_) + "_"
-                                          + boost::lexical_cast<string>(rank);
+      std::string fileName = "tmp_collection_" + std::to_string(myId_) + "_"
+                                               + std::to_string(rank);
 
       fileStream_.open( fileName.c_str(), std::ios::binary );
    }
@@ -126,8 +122,8 @@ void FileGatherScheme::collectFromFiles()
    //look for files from all processes
    for(int rank=0; rank < MPIManager::instance()->numProcesses(); ++rank)
    {
-      string fileName = "tmp_collection_" + boost::lexical_cast<string>(myId_) + "_"
-                                          + boost::lexical_cast<string>(rank);
+      string fileName = "tmp_collection_" + to_string(myId_) + "_"
+                                          + to_string(rank);
 
       if(exists(fileName)) {
 		  std::ifstream * str = new std::ifstream(fileName.c_str(),ios::in|ios::binary);
diff --git a/src/geometry/initializer/ScalarFieldFromBody.impl.h b/src/geometry/initializer/ScalarFieldFromBody.impl.h
index 98b6d2a167ad7c64b309da26a0e20d7db6c388a8..0a166fd0ef58eec0a610370415a5aee9ff6c968b 100644
--- a/src/geometry/initializer/ScalarFieldFromBody.impl.h
+++ b/src/geometry/initializer/ScalarFieldFromBody.impl.h
@@ -66,11 +66,11 @@ namespace initializer {
       
       try
       {
-         Value_T value = boost::lexical_cast<Value_T>(expression);
+         Value_T value = string_to_num<Value_T>(expression);
          init ( *bodyFromConfig ( subBlock ), value, addOrSet, id );
       }
       
-      catch(boost::bad_lexical_cast&)
+      catch(std::invalid_argument&)
       {
          math::FunctionParser p;
          p.parse(expression);         
diff --git a/src/geometry/initializer/ScalarFieldFromCellInterval.impl.h b/src/geometry/initializer/ScalarFieldFromCellInterval.impl.h
index d9b97544938b8553cc361b53728b999ba9612921..8e65dcfb475790f585984252b79d0e5a2aa6b896 100644
--- a/src/geometry/initializer/ScalarFieldFromCellInterval.impl.h
+++ b/src/geometry/initializer/ScalarFieldFromCellInterval.impl.h
@@ -62,11 +62,11 @@ void ScalarFieldFromCellInterval<Field_T>::init( const Config::BlockHandle & blo
    
    try
    {
-      Value_T value = boost::lexical_cast<Value_T>(expression);
+      Value_T value = string_to_num<Value_T>(expression);
       init(globalCellInterval, value, id);
    }
    
-   catch(boost::bad_lexical_cast&)
+   catch(std::invalid_argument&)
    {
       math::FunctionParser p;
       p.parse(expression);
diff --git a/src/pde/iterations/VCycles.impl.h b/src/pde/iterations/VCycles.impl.h
index f202d89e20b45e8c3cd26fb9398b0acbff6fb1cc..fea9a4031534af18fd812000f657f65c8e41033c 100644
--- a/src/pde/iterations/VCycles.impl.h
+++ b/src/pde/iterations/VCycles.impl.h
@@ -91,9 +91,9 @@ VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::
    for ( uint_t lvl = 1; lvl < numLvl; ++lvl )
    {
       auto getSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
-      uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
-      fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
-      rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
 
       for ( auto & w: weights_[lvl] )
       {
@@ -197,10 +197,10 @@ VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::
    for ( uint_t lvl = 1; lvl < numLvl; ++lvl )
    {
       auto getSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
-      uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
-      fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
-      rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
-      stencilId_.push_back( field::addToStorage< StencilField_T >( blocks, "w_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
+      stencilId_.push_back( field::addToStorage< StencilField_T >( blocks, "w_"+std::to_string(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
    }
 
    // CoarsenStencilFieldsDCA<Stencil_T>( blocks, stencilId_, numLvl, uint_t(2)) ();  // scaling by ( 1/h^2 )^lvl
diff --git a/src/pe/cr/HCSITS.impl.h b/src/pe/cr/HCSITS.impl.h
index 314b3d1c1dae7313d706053d66011ebf62b7e540..a03c3107fdc7e6162a21af9fd353ec19358e6d3e 100644
--- a/src/pe/cr/HCSITS.impl.h
+++ b/src/pe/cr/HCSITS.impl.h
@@ -1807,7 +1807,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::integratePositions( Body
             v = v * (edge * getSpeedLimitFactor() / dt / speed );
          }
 
-         const real_t maxPhi = real_t(2) * math::PI * getSpeedLimitFactor();
+         const real_t maxPhi = real_t(2) * math::M_PI * getSpeedLimitFactor();
          const real_t phi    = w.length() * dt;
          if (phi > maxPhi)
          {
diff --git a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
index 91ae0da40a0a12448f30e0c996924f6dc76ae6e1..8d7fc55c0fa014f6eddd1189618b5a985de1dc67 100644
--- a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
+++ b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
@@ -54,7 +54,7 @@ real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
 }
 
 // Coefficient from Stokes' law for drag, only valid for Stokes regime (low Reynolds numbers)
-// = 3 * PI * mu * D * fluidVolumeFraction
+// = 3 * M_PI * mu * D * fluidVolumeFraction
 real_t dragCoeffStokes ( real_t fluidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity )
 {
    return real_t(3) * math::M_PI * diameter * fluidDynamicViscosity * fluidVolumeFraction;
diff --git a/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h b/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
index b9bf10293b45ccbf7831871e2f1014cf783d3abd..cd036bfc064357692fa645c7ad99df529f5910f3 100644
--- a/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
+++ b/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
@@ -225,9 +225,9 @@ pe::Vec3 LubricationForceEvaluator::compLubricationSphrSphr( real_t gap, const p
    real_t d = real_t(2) * diameterSphereI * diameterSphereJ / ( diameterSphereI + diameterSphereJ );
    real_t h = gap;
    real_t r = d + h;
-   real_t a_sq = ( real_t(3) * dynamicViscosity_ * walberla::math::PI * d / real_t(2) ) * ( d / ( real_t(4) * h ) + ( real_t(18) / real_t(40) ) * std::log( d / ( real_t(2) * h ) ) +
+   real_t a_sq = ( real_t(3) * dynamicViscosity_ * walberla::math::M_PI * d / real_t(2) ) * ( d / ( real_t(4) * h ) + ( real_t(18) / real_t(40) ) * std::log( d / ( real_t(2) * h ) ) +
                                                                                             ( real_t(9)/real_t(84) ) * ( h / d ) * std::log( d/( real_t(2)*h ) ) );
-   real_t a_sh = ( dynamicViscosity_ * walberla::math::PI * d / real_t(2) ) * std::log( d / ( real_t(2) * h ) ) * ( d + h ) * ( d + h ) / real_t(4);
+   real_t a_sh = ( dynamicViscosity_ * walberla::math::M_PI * d / real_t(2) ) * std::log( d / ( real_t(2) * h ) ) * ( d + h ) * ( d + h ) / real_t(4);
    pe::Vec3 fLub( - a_sq * length * rIJ - a_sh * ( real_t(2) / r ) * ( real_t(2) / r ) * ( velDiff - length * rIJ ) );
 
    WALBERLA_LOG_DETAIL_SECTION()
@@ -274,9 +274,9 @@ pe::Vec3 LubricationForceEvaluator::compLubricationSphrPlane( real_t gap, const
    real_t d = real_t(4) * radiusSphereI;
    real_t h = gap;
    real_t r = d + h;
-   real_t a_sq = ( real_t(3) * dynamicViscosity_ * walberla::math::PI * d / real_t(2) ) * ( d / ( real_t(4) * h ) + ( real_t(18) / real_t(40) ) * std::log( d / ( real_t(2) * h ) ) +
+   real_t a_sq = ( real_t(3) * dynamicViscosity_ * walberla::math::M_PI * d / real_t(2) ) * ( d / ( real_t(4) * h ) + ( real_t(18) / real_t(40) ) * std::log( d / ( real_t(2) * h ) ) +
                                                                                             ( real_t(9)/real_t(84) ) * ( h / d ) * std::log( d/( real_t(2)*h ) ) );
-   real_t a_sh = ( dynamicViscosity_ * walberla::math::PI * d / real_t(2) ) * std::log( d / ( real_t(2) * h ) ) * ( d + h ) * ( d + h ) / real_t(4);
+   real_t a_sh = ( dynamicViscosity_ * walberla::math::M_PI * d / real_t(2) ) * std::log( d / ( real_t(2) * h ) ) * ( d + h ) * ( d + h ) / real_t(4);
    pe::Vec3 fLub( - a_sq * length * rIJ - a_sh * ( real_t(2) / r ) * ( real_t(2) / r ) * ( v1 - length * rIJ ) );
 
    WALBERLA_LOG_DETAIL_SECTION() {
diff --git a/src/pe_coupling/utility/LubricationCorrection.cpp b/src/pe_coupling/utility/LubricationCorrection.cpp
index eaf1952149fc4c35d7b6b45b23549c3516557f76..6cee7eb439df47af3c8051e9a8fb5d2c882e175d 100644
--- a/src/pe_coupling/utility/LubricationCorrection.cpp
+++ b/src/pe_coupling/utility/LubricationCorrection.cpp
@@ -189,7 +189,7 @@ pe::Vec3 LubricationCorrection::compLubricationSphrSphr( real_t gap, const pe::S
    real_t radiiSQR    = ( radiusSphereI * radiusSphereJ ) * ( radiusSphereI * radiusSphereJ );
    real_t radiiSumSQR = ( radiusSphereI + radiusSphereJ ) * ( radiusSphereI + radiusSphereJ );
 
-   pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+   pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::M_PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
 
    WALBERLA_LOG_DETAIL_SECTION()
    {
@@ -247,7 +247,7 @@ pe::Vec3 LubricationCorrection::compLubricationSphrPlane( real_t gap, const pe::
 
    real_t radiiSQR = radiusSphereI * radiusSphereI;
 
-   pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+   pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::M_PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
 
    WALBERLA_LOG_DETAIL_SECTION() {
       std::stringstream ss;
diff --git a/src/vtk/Initialization.cpp b/src/vtk/Initialization.cpp
index d6df0a4fb137d35fb8e4c8834d137bad99029b35..7109d9e294c23cd310f88b0cd51297ed8d595eed 100644
--- a/src/vtk/Initialization.cpp
+++ b/src/vtk/Initialization.cpp
@@ -49,9 +49,9 @@ static void splitVector( T& x, T& y, T& z, const Config::BlockHandle& bb, const
    if( coordinates.size() != 3 )
       WALBERLA_ABORT( errorMsg );
 
-   x = boost::lexical_cast< T >( coordinates[0] );
-   y = boost::lexical_cast< T >( coordinates[1] );
-   z = boost::lexical_cast< T >( coordinates[2] );
+   x = string_to_num< T >( coordinates[0] );
+   y = string_to_num< T >( coordinates[1] );
+   z = string_to_num< T >( coordinates[2] );
 }
 
 
diff --git a/tests/core/debug/CheckMacroTest.cpp b/tests/core/debug/CheckMacroTest.cpp
index d2e4cf7a8ed94de05f876af83cc3e6ac67280fa3..2b11a23c0e599b743f747478e9ca6b36297f90c0 100644
--- a/tests/core/debug/CheckMacroTest.cpp
+++ b/tests/core/debug/CheckMacroTest.cpp
@@ -27,6 +27,7 @@
 #include <iostream>
 #include <vector>
 #include <map>
+#include <algorithm>
 
 
 // Own data structure with comparison operators but without operator<<
diff --git a/tests/core/load_balancing/MetisTest.cpp b/tests/core/load_balancing/MetisTest.cpp
index 22038220a93d192e4ac06c88a5fb220e5c415f8d..5db60e5ae786e14811fab22cd87cfccb903fa57b 100644
--- a/tests/core/load_balancing/MetisTest.cpp
+++ b/tests/core/load_balancing/MetisTest.cpp
@@ -39,8 +39,6 @@
 #include "vtk/VTKOutput.h"
 #include "vtk/BlockCellDataWriter.h"
 
-#include <boost/lexical_cast.hpp>
-
 int main( int argc, char * argv[] )
 {
    using namespace walberla;
@@ -62,8 +60,8 @@ int main( int argc, char * argv[] )
    }
 
    try {
-      fieldSize.set( boost::lexical_cast< uint_t >( args.at(1) ), boost::lexical_cast< uint_t >( args.at(2) ) );
-      partitions = boost::lexical_cast< uint_t >( args.at(3) );
+      fieldSize.set( string_to_num< uint_t >( args.at(1) ), string_to_num< uint_t >( args.at(2) ) );
+      partitions = string_to_num< uint_t >( args.at(3) );
    }
    catch( std::exception & e )
    {
diff --git a/tests/core/load_balancing/ParMetisTest.cpp b/tests/core/load_balancing/ParMetisTest.cpp
index fd72061009c5e4b26648e38473bd0f3185800ecc..a9e79d3848dd327abc2d3682b12ed633b1e386fd 100644
--- a/tests/core/load_balancing/ParMetisTest.cpp
+++ b/tests/core/load_balancing/ParMetisTest.cpp
@@ -41,8 +41,6 @@
 #include "vtk/VTKOutput.h"
 #include "vtk/BlockCellDataWriter.h"
 
-#include <boost/lexical_cast.hpp>
-
 int main( int argc, char * argv[] )
 {
    using namespace walberla;
@@ -57,8 +55,8 @@ int main( int argc, char * argv[] )
    bool            vtk = true;
 
    try {
-      fieldSize.set( boost::lexical_cast< uint_t >( args.at(1) ), boost::lexical_cast< uint_t >( args.at(2) ) );
-      partitions = boost::lexical_cast< uint_t >( args.at(3) );
+      fieldSize.set( string_to_num< uint_t >( args.at(1) ), string_to_num< uint_t >( args.at(2) ) );
+      partitions = string_to_num< uint_t >( args.at(3) );
 
       auto it = std::find( args.begin(), args.end(), "--no-vtk" );
       if(it != args.end())
diff --git a/tests/core/math/GenericAABBTest.cpp b/tests/core/math/GenericAABBTest.cpp
index d574ba15343cb932633661886ee2cc6c8cd582dd..ee6db55766e9528b00de3a7a6bb323ca0677ee3b 100644
--- a/tests/core/math/GenericAABBTest.cpp
+++ b/tests/core/math/GenericAABBTest.cpp
@@ -28,6 +28,7 @@
 #include "stencil/D3CornerStencil.h"
 
 #include <random>
+#include <array>
 
 
 using namespace walberla;
diff --git a/tests/core/math/PhysicalCheckTest.cpp b/tests/core/math/PhysicalCheckTest.cpp
index 2ad726467cc94b6d72a5fc5ef0e8d44a29010565..abb019c4ad4dad6de41989817da5ce7bd05ef4e2 100644
--- a/tests/core/math/PhysicalCheckTest.cpp
+++ b/tests/core/math/PhysicalCheckTest.cpp
@@ -26,7 +26,6 @@
 #include "core/logging/Logging.h"
 #include "core/math/PhysicalCheck.h"
 
-#include <boost/lexical_cast.hpp>
 #include <iostream>
 #include <map>
 #include <string>
diff --git a/tests/core/math/PlaneTest.cpp b/tests/core/math/PlaneTest.cpp
index 98cfdb45eca908c2a137488089a0a03314626ff8..2f895119184e14b303536a3a84a44cdd707a8d68 100644
--- a/tests/core/math/PlaneTest.cpp
+++ b/tests/core/math/PlaneTest.cpp
@@ -21,6 +21,7 @@
 
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "core/math/Plane.h"
 #include "core/math/Vector3.h"
 #include "core/mpi/Environment.h"
@@ -137,7 +138,7 @@ int main(int argc, char * argv[])
 
       real_t angle = std::acos( (p1-p0) * (p2-p0) / std::sqrt( (p1-p0).sqrLength() * (p2-p0).sqrLength() ) );
 
-      if( (p0 - p1).sqrLength() < 1e-6 || (p0 - p2).sqrLength() < 1e-6 || (p2 - p1).sqrLength() < 1e-6 || angle < math::PI / real_t(180) )
+      if( (p0 - p1).sqrLength() < 1e-6 || (p0 - p2).sqrLength() < 1e-6 || (p2 - p1).sqrLength() < 1e-6 || angle < math::M_PI / real_t(180) )
       {
          --i;
          continue;
diff --git a/tests/core/math/equation_system/EquationSolverTest.cpp b/tests/core/math/equation_system/EquationSolverTest.cpp
index 4319769cfc8790368de21be1d261d6b5f13d4e04..e711d04f8300e746ce14d922d7103c3afe3b0670 100644
--- a/tests/core/math/equation_system/EquationSolverTest.cpp
+++ b/tests/core/math/equation_system/EquationSolverTest.cpp
@@ -26,7 +26,6 @@
 #include "core/math/equation_system/EquationParser.h"
 #include "core/math/equation_system/EquationSystem.h"
 
-#include <boost/lexical_cast.hpp>
 #include <iostream>
 #include <string>
 #include <vector>
@@ -128,7 +127,7 @@ int equationInput(){
 
    for (size_t i=0; i<eqStringList.size(); ++i){
       index = 0;
-      es.add( boost::lexical_cast<std::string>(++number), ep.parseEquation( eqStringList[i], index ) );
+      es.add( std::to_string(++number), ep.parseEquation( eqStringList[i], index ) );
    }
 
    WALBERLA_CHECK( es.solve() );
@@ -139,7 +138,7 @@ int equationInput(){
 
 int unitTest(double v)
 {
-   std::string s = boost::lexical_cast<std::string>( v );
+   std::string s = std::to_string( v );
 
    std::vector<std::string> eqStringList;
    eqStringList.push_back(       "a = " + s );
diff --git a/tests/core/mpi/MPITextFileTest.cpp b/tests/core/mpi/MPITextFileTest.cpp
index 1500d730f0fdbada5b765286d26b27e5727fa438..a0ed6e61cc4429ca9f5369b71ec908a8fd4a5bb9 100644
--- a/tests/core/mpi/MPITextFileTest.cpp
+++ b/tests/core/mpi/MPITextFileTest.cpp
@@ -120,7 +120,7 @@ int main( int argc, char * argv[] )
    std::string filename;
    try
    {
-      chunkSize = boost::lexical_cast<size_t>( args.at(2) );
+      chunkSize = walberla::string_to_num<size_t>( args.at(2) );
       filename  = args.at( 1 );
    }
    catch( ... )
diff --git a/tests/core/timing/TimerTest.cpp b/tests/core/timing/TimerTest.cpp
index 4df41d6a156d5f214c59d13b79cde8659af21931..a6ef83f9ee0fe8680089813d429c575ebf997acf 100644
--- a/tests/core/timing/TimerTest.cpp
+++ b/tests/core/timing/TimerTest.cpp
@@ -21,6 +21,7 @@
 
 
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/timing/Timer.h"
 
 #include <cmath>
@@ -30,7 +31,7 @@ namespace walberla {
 static double burnTime()
 {
    double sum = 0.0;
-   for( double d = 0.0; d < math::PI; d += 0.000001 )
+   for( double d = 0.0; d < math::M_PI; d += 0.000001 )
    {
       sum += std::atan( std::tan( d ) );
       sum += std::asin( std::sin( d ) );
@@ -80,4 +81,4 @@ int main()
 {
   walberla::debug::enterTestMode();
   walberla::mergeTest();
-}
\ No newline at end of file
+}
diff --git a/tests/core/timing/TimingPoolTest.cpp b/tests/core/timing/TimingPoolTest.cpp
index 11f5f5b88562e3ecc883330e1602e1aae1cc3d27..65bc2ca5ac679fe703bc3721aa301498619f0b81 100644
--- a/tests/core/timing/TimingPoolTest.cpp
+++ b/tests/core/timing/TimingPoolTest.cpp
@@ -21,6 +21,7 @@
 
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "core/mpi/MPIManager.h"
 #include "core/timing/Timer.h"
 #include "core/timing/TimingPool.h"
@@ -61,7 +62,7 @@ void scopedTimer()
       auto scopedTimer = pool.getScopeTimer( "scope timer" );
 
       double sum = 0.0;
-      for( double d = 0.0; d < math::PI; d += 0.00001 )
+      for( double d = 0.0; d < math::M_PI; d += 0.00001 )
       {
          sum += std::atan( std::tan( d ) );
          sum += std::asin( std::sin( d ) );
diff --git a/tests/fft/GreensTest.cpp b/tests/fft/GreensTest.cpp
index 720640d7cd0eb0db8d1c2aeaaecd947be6dd5888..aaa669c2ac0f79048e379e09919afdba22412ec6 100644
--- a/tests/fft/GreensTest.cpp
+++ b/tests/fft/GreensTest.cpp
@@ -3,6 +3,7 @@
 #include "core/Environment.h"
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "field/AddToStorage.h"
 #include "blockforest/communication/UniformBufferedScheme.h"
 #include "stencil/D3Q7.h"
@@ -55,9 +56,9 @@ int main (int argc, char** argv)
    auto greens = [&dim] (uint_t x, uint_t y, uint_t z) -> real_t {
       if (x == 0 && y == 0 && z == 0)
          return 0;
-      return real_c(0.5) / ( std::cos( real_c(2) * real_c(M_PI) * real_c(x) / real_c(dim[0])) +
-                             std::cos( real_c(2) * real_c(M_PI) * real_c(y) / real_c(dim[1])) +
-                             std::cos( real_c(2) * real_c(M_PI) * real_c(z) / real_c(dim[2])) -
+      return real_c(0.5) / ( std::cos( real_c(2) * real_c(math::M_PI) * real_c(x) / real_c(dim[0])) +
+                             std::cos( real_c(2) * real_c(math::M_PI) * real_c(y) / real_c(dim[1])) +
+                             std::cos( real_c(2) * real_c(math::M_PI) * real_c(z) / real_c(dim[2])) -
                              real_c(3) ) / real_c(dim[0]*dim[1]*dim[2]);
    };
    
diff --git a/tests/field/FieldFileIOTest.cpp b/tests/field/FieldFileIOTest.cpp
index 6e7a95a2dec44fb408c17e1039df77867e8fdb0b..1214f01992c75d6a1226ea2306e9386848fe7faf 100644
--- a/tests/field/FieldFileIOTest.cpp
+++ b/tests/field/FieldFileIOTest.cpp
@@ -33,8 +33,6 @@
 #include "field/Field.h"
 #include "field/FileIO.h"
 
-#include <boost/lexical_cast.hpp>
-
 
 namespace mpi_file_io_test {
    
@@ -86,10 +84,10 @@ int main( int argc, char* argv[] )
    
    if( args.size() == 5 )
    {
-      numBlocks  = boost::lexical_cast<uint_t>( args[1] );
-      xBlockSize = boost::lexical_cast<uint_t>( args[2] );
-      yBlockSize = boost::lexical_cast<uint_t>( args[3] );
-      zBlockSize = boost::lexical_cast<uint_t>( args[4] );
+      numBlocks  = string_to_num<uint_t>( args[1] );
+      xBlockSize = string_to_num<uint_t>( args[2] );
+      yBlockSize = string_to_num<uint_t>( args[3] );
+      zBlockSize = string_to_num<uint_t>( args[4] );
    }
    else if( args.size() > 5 )
    {
diff --git a/tests/geometry/ScalarFieldFromBodyTest.cpp b/tests/geometry/ScalarFieldFromBodyTest.cpp
index d6c29f3d1b4ddd97b36155b78e9a8087ea426867..0aa382715a8d664570635d958a86fd3f23bc58cc 100644
--- a/tests/geometry/ScalarFieldFromBodyTest.cpp
+++ b/tests/geometry/ScalarFieldFromBodyTest.cpp
@@ -29,6 +29,7 @@
 
 #include "core/debug/TestSubsystem.h"
 #include "core/logging/Logging.h"
+#include "core/math/Constants.h"
 #include "core/math/Vector3.h"
 
 #include "field/AddToStorage.h"
@@ -47,7 +48,6 @@ using namespace geometry;
 const uint_t confBlockCount []      = { 1, 1, 1 };
 const uint_t confCells []           = { 30, 30, 30 };
 
-const real_t PI = real_t(3.1415927);
 const bool   useGui = false;
 
 
@@ -156,7 +156,7 @@ void ellipsoidTest( StructuredBlockStorage & storage,
    const Vector3<real_t> axis1 (  real_t(1), real_t(1), real_t(0) );
    const Vector3<real_t> axis2 ( -real_t(1), real_t(1), real_t(0) );
 
-   const real_t expectedVolume = real_t(4) / real_t(3) * radii[0] * radii[1] * radii[2] * PI;
+   const real_t expectedVolume = real_t(4) / real_t(3) * radii[0] * radii[1] * radii[2] * math::M_PI;
 
    resetField( storage, fieldID );
 
@@ -209,7 +209,7 @@ void sphereTest( StructuredBlockStorage & storage,
 {
    const Vector3<real_t> midpoint ( real_t(15), real_t(15), real_t(15) );
    const real_t radius = real_t(5);
-   const real_t expectedVolume = real_t(4) / real_t(3) * radius * radius * radius * PI;
+   const real_t expectedVolume = real_t(4) / real_t(3) * radius * radius * radius * math::M_PI;
 
    resetField( storage, fieldID );
 
@@ -294,4 +294,4 @@ int main( int argc, char ** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/geometry/VoxelFileTest.cpp b/tests/geometry/VoxelFileTest.cpp
index 726616dc11a6c68ef962ee5ac9bdfe462b95ce81..d8ab6ddb679cff372a72fea03d056b9d0e2bd523 100644
--- a/tests/geometry/VoxelFileTest.cpp
+++ b/tests/geometry/VoxelFileTest.cpp
@@ -21,18 +21,6 @@
 
 
 
-#ifdef _MSC_VER
-#  pragma warning(push)
-// disable warning boost/concept/detail/msvc.hpp(22): warning C4100: 'x' : unreferenced formal parameter
-#  pragma warning( disable : 4100 )
-// disable warning bboost/concept/detail/msvc.hpp(76): warning C4067 : unexpected tokens following preprocessor directive - expected a newline
-#  pragma warning( disable : 4067 )
-
-#include <boost/concept/detail/msvc.hpp>
-
-#  pragma warning(pop)
-#endif //_MSC_VER
-
 #ifdef _MSC_VER
 #  pragma warning(push)
 // disable warning boost/multi_array/base.hpp(475): warning C4189: 'bound_adjustment' : local variable is initialized but not referenced
diff --git a/tests/lbm/DiffusionTest.cpp b/tests/lbm/DiffusionTest.cpp
index c1aa845088b935ef751deb37c31df158dc5bab13..bdcb97ff68a2ce57a9a32a15671d5986a63024c5 100644
--- a/tests/lbm/DiffusionTest.cpp
+++ b/tests/lbm/DiffusionTest.cpp
@@ -59,6 +59,7 @@
 #include "core/cell/CellVector.h"
 #include "core/config/Config.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/math/IntegerFactorization.h"
 #include "core/mpi/Environment.h"
 
@@ -77,8 +78,6 @@
 
 #include "vtk/VTKOutput.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <functional>
 
 
@@ -198,14 +197,14 @@ int run( int argc, char **argv )
    if( argc > 1 ) {
       std::vector<std::string> args( argv, argv + argc );
       for( uint_t i = 1; i < uint_c(argc); ++i ) {
-              if( boost::equals(argv[i], "-d"      ) )   d      = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-dim"    ) )   dim    = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-dx"     ) )   dx     = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-dt"     ) )   dt     = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-dv"     ) )   dv     = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-v"      ) )   u_in   = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-t"      ) )   time   = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-err"    ) )   err    = boost::lexical_cast<real_t>( args[++i] );
+              if( boost::equals(argv[i], "-d"      ) )   d      = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-dim"    ) )   dim    = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-dx"     ) )   dx     = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-dt"     ) )   dt     = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-dv"     ) )   dv     = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-v"      ) )   u_in   = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-t"      ) )   time   = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-err"    ) )   err    = string_to_num<real_t>( args[++i] );
          else if( boost::equals(argv[i], "--gui"   ) )   useGui = true;
          else if( boost::equals(argv[i], "--quiet" ) )   quiet  = true;
          else if( boost::equals(argv[i], "--vtk"   ) )   useVTK = true;
@@ -236,8 +235,8 @@ int run( int argc, char **argv )
    if(!quiet) WALBERLA_LOG_RESULT( " -> u   = " << u   );
    if(!quiet) WALBERLA_LOG_RESULT( " -> tau = " << tau );
 
-   const real_t tperiod = real_t(2) * math::PI / real_c( timesteps  );
-   const real_t cperiod = real_t(2) * math::PI / real_c( cells[dim] );
+   const real_t tperiod = real_t(2) * math::M_PI / real_c( timesteps  );
+   const real_t cperiod = real_t(2) * math::M_PI / real_c( cells[dim] );
 
    // --- create blockstorage --- //
    auto blockStorage = blockforest::createUniformBlockGrid(
@@ -347,7 +346,7 @@ int main(int argc, char **argv)
    bool corr = true;
    for( int i=0; i<argc; ++i ){
       if( boost::equals(argv[i], "-c" ) ){
-         corr = boost::lexical_cast<bool>( argv[++i] );
+         corr = std::atoi( argv[++i] ) != 0;
          break;
       }
    }
diff --git a/tests/lbm/boundary/BoundaryForce.cpp b/tests/lbm/boundary/BoundaryForce.cpp
index 08b41401ca9641f3daab0dd30ff656a473b165f5..c36db2ec17579532697f6ea9897f372d5da1ecd0 100644
--- a/tests/lbm/boundary/BoundaryForce.cpp
+++ b/tests/lbm/boundary/BoundaryForce.cpp
@@ -25,6 +25,7 @@
 
 #include "boundary/BoundaryHandling.h"
 
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 #include "core/mpi/Reduce.h"
@@ -167,7 +168,7 @@ int main( int argc, char ** argv )
    
    uint_t timeSteps = uint_c(200);
    if( argc > 1 )
-      timeSteps = boost::lexical_cast<uint_t>( argv[1] );
+      timeSteps = uint_c(std::stoul( argv[1] ));
    SweepTimeloop timeloop( blocks->getBlockStorage(), timeSteps );
  
    blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
@@ -211,7 +212,7 @@ int main( int argc, char ** argv )
    mpi::allReduceInplace( force[2], mpi::SUM );
    
    real_t visc = lbm::collision_model::viscosityFromOmega( omega );
-   Vector3<real_t> stokes = 6 * math::PI * visc * R * velocity;
+   Vector3<real_t> stokes = 6 * math::M_PI * visc * R * velocity;
 
    WALBERLA_LOG_RESULT_ON_ROOT("Expected force: " << stokes);
    WALBERLA_LOG_RESULT_ON_ROOT("Actual force: " << force);
diff --git a/tests/lbm/boundary/BoundaryForceCouette.cpp b/tests/lbm/boundary/BoundaryForceCouette.cpp
index 1acf8f8c729f002b90eb575f2f5b2f654d4fea11..974015b9d5c0d64e705694f28cff9182b182d4e1 100644
--- a/tests/lbm/boundary/BoundaryForceCouette.cpp
+++ b/tests/lbm/boundary/BoundaryForceCouette.cpp
@@ -162,7 +162,7 @@ int main( int argc, char ** argv )
    
    uint_t timeSteps = uint_c(2000);
    if( argc > 1 )
-      timeSteps = boost::lexical_cast<uint_t>( argv[1] );
+      timeSteps = uint_c(std::stoul( argv[1] ));
    SweepTimeloop timeloop( blocks->getBlockStorage(), timeSteps );
  
    blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
diff --git a/tests/lbm/boundary/DiffusionDirichlet.cpp b/tests/lbm/boundary/DiffusionDirichlet.cpp
index 79cee7c428ced16bbf99006d0011cee3be53e4a5..9661765230757b846acc2aa41b4e45e9a9ad8514 100644
--- a/tests/lbm/boundary/DiffusionDirichlet.cpp
+++ b/tests/lbm/boundary/DiffusionDirichlet.cpp
@@ -56,6 +56,7 @@
 #include "core/cell/CellVector.h"
 #include "core/config/Config.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/math/IntegerFactorization.h"
 #include "core/mpi/Environment.h"
 
@@ -79,8 +80,6 @@
 
 #include "vtk/VTKOutput.h"
 
-#include <boost/lexical_cast.hpp>
-
 
 namespace walberla {
 
@@ -107,7 +106,7 @@ public:
    using cplx_t = std::complex<real_t>;
 
    PlugFlow( real_t L, real_t H, real_t u, real_t k ) :
-      period_( real_t(2)*math::PI/L ),
+      period_( real_t(2)*math::M_PI/L ),
       lambda_( period_*sqrt( cplx_t(real_t(1), u/k/period_) ) ),
       emH_   ( real_t(1) - exp(-lambda_*H) ),
       epH_   ( real_t(1) - exp(+lambda_*H) ),
@@ -220,12 +219,12 @@ int main( int argc, char **argv )
    if( argc > 1 ) {
       std::vector<std::string> args( argv, argv + argc );
       for( uint_t i = 1; i < uint_c(argc); ++i ) {
-              if( boost::equals(argv[i], "-o"    ) ) omega  = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-l"    ) ) length = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-w"    ) ) width  = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-v"    ) ) velx   = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-t"    ) ) time   = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-e"    ) ) error  = boost::lexical_cast<real_t>( args[++i] );
+              if( boost::equals(argv[i], "-o"    ) ) omega  = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-l"    ) ) length = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-w"    ) ) width  = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-v"    ) ) velx   = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-t"    ) ) time   = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-e"    ) ) error  = string_to_num<real_t>( args[++i] );
          else if( boost::equals(argv[i], "--gui" ) ) useGui = true;
          else if( boost::equals(argv[i], "--vtk" ) ) useVTK = true;
          else if( argv[i][0] != '-' ){
@@ -252,7 +251,7 @@ int main( int argc, char **argv )
 
    BlockDataID boundaryHandling = MyBoundaryHandling::addDefaultDiffusionBoundaryHandlingToStorage( blockStorage, "BoundaryHandling", flagFieldID, getFluidFlag(), srcFieldID );
 
-   auto cbc = make_shared<CosBoundaryConfiguration>( real_t(2)*math::PI/real_c(length) );
+   auto cbc = make_shared<CosBoundaryConfiguration>( real_t(2)*math::M_PI/real_c(length) );
    geometry::initializer::BoundaryFromDomainBorder<MyBoundaryHandling::BoundaryHandling_T> bfdb( *blockStorage, boundaryHandling );
    bfdb.init( MyBoundaryHandling::getDiffusionDirichletBoundaryUID(), stencil::N, cbc, -1, 1 );
    bfdb.init( MyBoundaryHandling::getDiffusionDirichletBoundaryUID(), stencil::S, cbc, -1, 1 );
diff --git a/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp b/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp
index 30c95b82abe998b1ced956ea650fd6ec7c9f7e8a..bb86d552a3262adfeb99a121f94ed480d7ac9002 100644
--- a/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp
+++ b/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp
@@ -47,6 +47,7 @@
 #include "core/cell/CellVector.h"
 #include "core/config/Config.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/math/IntegerFactorization.h"
 #include "core/math/Limits.h"
 #include "core/mpi/Environment.h"
@@ -72,8 +73,6 @@
 
 #include "timeloop/SweepTimeloop.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <stdexcept>
 #include <array>
 #include <functional>
@@ -237,8 +236,8 @@ public:
          delta_( maxValue - minValue),
          length_(real_c(length)),
          lengthInv_(real_t(1)/real_c(length)),
-         pi_(math::PI),
-         piInv_(real_t(1)/math::PI),
+         pi_(math::M_PI),
+         piInv_(real_t(1)/math::M_PI),
          valid_(uint_c(std::ceil(omega*omega*omega*real_t(10)))),
          time_( time ),
          expArray(),
@@ -430,13 +429,13 @@ int main( int argc, char **argv )
    if( argc > 1 ) {
       std::vector<std::string> args( argv, argv + argc );
       for( uint_t i = 1; i < uint_c(argc); ++i ) {
-              if( boost::equals(argv[i], "-l"    ) )   length  = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-w"    ) )   width   = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-t"    ) )   time    = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-d"    ) )   dv      = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-o"    ) )   omega   = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-c"    ) )   closed  = boost::lexical_cast< bool >( args[++i] );
-         else if( boost::equals(argv[i], "-r"    ) )   levels += boost::lexical_cast<uint_t>( args[++i] );
+              if( boost::equals(argv[i], "-l"    ) )   length  = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-w"    ) )   width   = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-t"    ) )   time    = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-d"    ) )   dv      = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-o"    ) )   omega   = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-c"    ) )   closed  = string_to_num<int>( args[++i] ) != 0;
+         else if( boost::equals(argv[i], "-r"    ) )   levels += string_to_num<uint_t>( args[++i] );
          else if( boost::equals(argv[i], "--vtk" ) )   useVTK  = true;
          else if( argv[i][0] != '-' ){
             std::cerr << "Usage: -option value" << std::endl; return EXIT_FAILURE;
diff --git a/tests/lbm/boundary/SimplePABTest.cpp b/tests/lbm/boundary/SimplePABTest.cpp
index 31bc21a747617028e54cd2be91bfa537a60ca498..6d5b8c478e11b09873636ee68040b24b91b204a6 100644
--- a/tests/lbm/boundary/SimplePABTest.cpp
+++ b/tests/lbm/boundary/SimplePABTest.cpp
@@ -56,7 +56,6 @@
 
 #include "vtk/VTKOutput.h"
 
-#include <boost/lexical_cast.hpp>
 #include <stdexcept>
 
 
@@ -218,11 +217,11 @@ int main( int argc, char **argv )
       if( args.size() != 6 )
          throw std::invalid_argument( "Wrong number of command line arguments!" );
 
-      channelLength = boost::lexical_cast<uint_t>( args[1] );
-      channelWidth  = boost::lexical_cast<uint_t>( args[2] );
-      omega         = boost::lexical_cast<real_t>( args[3] );
-      deltaDensity  = boost::lexical_cast<real_t>( args[4] );
-      numTimesteps  = boost::lexical_cast<uint_t>( args[5] );
+      channelLength = string_to_num<uint_t>( args[1] );
+      channelWidth  = string_to_num<uint_t>( args[2] );
+      omega         = string_to_num<real_t>( args[3] );
+      deltaDensity  = string_to_num<real_t>( args[4] );
+      numTimesteps  = string_to_num<uint_t>( args[5] );
    }
    catch( std::exception & )
    {
diff --git a/tests/lbm/evaluations/PermeabilityTest.cpp b/tests/lbm/evaluations/PermeabilityTest.cpp
index 022ad27889a260167a7d900f50ee64b41efbe308..7826539357a236d6cc3976631e9cf6163fcd65fc 100644
--- a/tests/lbm/evaluations/PermeabilityTest.cpp
+++ b/tests/lbm/evaluations/PermeabilityTest.cpp
@@ -107,7 +107,7 @@ real_t permeability( Setup setup )
    for( uint_t i = 0; i < 31; i++ )
       drag += real_c(qs[i]) * real_c(std::pow( setup.kappa, real_c(i) ));
 
-   return ( L * L * L ) / ( real_t(6) * math::PI * r * real_t(2) * drag );
+   return ( L * L * L ) / ( real_t(6) * math::M_PI * r * real_t(2) * drag );
 }
 
 
diff --git a/tests/lbm/refinement/NonConstantDiffusion.cpp b/tests/lbm/refinement/NonConstantDiffusion.cpp
index 83ac2f989663821ac7f2db15a49c904d5cef42ca..592d23ce7fbe9e8535e0076d95960b0b40d23778 100644
--- a/tests/lbm/refinement/NonConstantDiffusion.cpp
+++ b/tests/lbm/refinement/NonConstantDiffusion.cpp
@@ -72,8 +72,6 @@
 
 #include "timeloop/SweepTimeloop.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <stdexcept>
 #include <functional>
 
@@ -275,15 +273,15 @@ int main( int argc, char **argv )
    if( argc > 1 ) {
       std::vector<std::string> args( argv, argv + argc );
       for( uint_t i = 1; i < uint_c(argc); ++i ) {
-              if( boost::equals(argv[i], "-l"    ) )   length  = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-w"    ) )   width   = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-t"    ) )   time    = boost::lexical_cast<uint_t>( args[++i] );
-         else if( boost::equals(argv[i], "-dv"   ) )   dv      = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-v"    ) )   v       = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-do"   ) )   domega  = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-o"    ) )   omega   = boost::lexical_cast<real_t>( args[++i] );
-         else if( boost::equals(argv[i], "-c"    ) )   closed  = boost::lexical_cast< bool >( args[++i] );
-         else if( boost::equals(argv[i], "-r"    ) )   levels += boost::lexical_cast<uint_t>( args[++i] );
+              if( boost::equals(argv[i], "-l"    ) )   length  = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-w"    ) )   width   = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-t"    ) )   time    = string_to_num<uint_t>( args[++i] );
+         else if( boost::equals(argv[i], "-dv"   ) )   dv      = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-v"    ) )   v       = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-do"   ) )   domega  = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-o"    ) )   omega   = string_to_num<real_t>( args[++i] );
+         else if( boost::equals(argv[i], "-c"    ) )   closed  = string_to_num<int>( args[++i] ) != 0;
+         else if( boost::equals(argv[i], "-r"    ) )   levels += string_to_num<uint_t>( args[++i] );
          else if( boost::equals(argv[i], "--vtk" ) )   useVTK  = true;
          else if( argv[i][0] != '-' ){
             std::cerr << "Usage: -option value" << std::endl; return EXIT_FAILURE;
diff --git a/tests/mesh/MeshAABBIntersectionTest.cpp b/tests/mesh/MeshAABBIntersectionTest.cpp
index fc4711134edbc3068b2011262fa0128330cc4ed6..a516342dee1533d29af1200ca372bf0f05537909 100644
--- a/tests/mesh/MeshAABBIntersectionTest.cpp
+++ b/tests/mesh/MeshAABBIntersectionTest.cpp
@@ -30,7 +30,6 @@
 #include "mesh/TriangleMeshes.h"
 
 #include <random>
-#include <boost/lexical_cast.hpp>
 
 #include <vector>
 #include <string>
@@ -88,7 +87,7 @@ int main( int argc, char * argv[] )
    if( args.size() != 2 )
       WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args[0] << " NUM_AABBS" );
 
-   const uint_t numAABBs = boost::lexical_cast< uint_t >( args[1] );
+   const uint_t numAABBs = string_to_num< uint_t >( args[1] );
 
    runTests< mesh::TriangleMesh >( numAABBs );
    runTests< mesh::FloatTriangleMesh >( numAABBs );
diff --git a/tests/mesh/MeshAABBSelectionTest.cpp b/tests/mesh/MeshAABBSelectionTest.cpp
index 37d421bfe5656e1ec59efe9b8d6e36463f3aa32b..df5b67b5343964f17da3d5c771227029b2e68aca 100644
--- a/tests/mesh/MeshAABBSelectionTest.cpp
+++ b/tests/mesh/MeshAABBSelectionTest.cpp
@@ -32,7 +32,6 @@
 #include "mesh/MeshOperations.h"
 
 #include <random>
-#include <boost/lexical_cast.hpp>
 
 #include <algorithm>
 #include <vector>
@@ -52,8 +51,8 @@ int main( int argc, char * argv[] )
       WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args[0] << " MESH_FILE NUM_BOXES NUM_POINTS_TESTED_PER_BOX" );
 
    const std::string & meshFile = args[1];
-   const uint_t numBoxes = boost::lexical_cast<uint_t>( args[2] );
-   const uint_t numPointsTestedPerBox = boost::lexical_cast<uint_t>( args[3] );
+   const uint_t numBoxes = string_to_num<uint_t>( args[2] );
+   const uint_t numPointsTestedPerBox = string_to_num<uint_t>( args[3] );
 
    auto mesh = make_shared<TriangleMesh>();
    mesh::readAndBroadcast( meshFile, *mesh );
diff --git a/tests/mesh/MeshBlockExclusionTest.cpp b/tests/mesh/MeshBlockExclusionTest.cpp
index 7685ac644dc38bbbbbd408f8ca834aae93214aef..996fdb877a20c6855f50f557426f72df73f39bf2 100644
--- a/tests/mesh/MeshBlockExclusionTest.cpp
+++ b/tests/mesh/MeshBlockExclusionTest.cpp
@@ -39,8 +39,6 @@
 
 #include "mesh/distance_octree/DistanceOctree.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <vector>
 #include <string>
 
@@ -150,7 +148,7 @@ int main( int argc, char * argv[] )
       WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args[0] << " MESH_FILE NUM_BLOCKS" );
 
    const std::string & meshFile       = args[1];
-   const uint_t        numTotalBlocks = boost::lexical_cast< uint_t >( args[2] );
+   const uint_t        numTotalBlocks = string_to_num< uint_t >( args[2] );
 
    run< mesh::TriangleMesh >( meshFile, numTotalBlocks );
    run< mesh::FloatTriangleMesh >( meshFile, numTotalBlocks );
diff --git a/tests/mesh/MeshDistancePlausibilityTest.cpp b/tests/mesh/MeshDistancePlausibilityTest.cpp
index 5ca33efc11439636d2df20583682947407f18580..e0fd1a276206c92c007119ea12c8a098bfe89502 100644
--- a/tests/mesh/MeshDistancePlausibilityTest.cpp
+++ b/tests/mesh/MeshDistancePlausibilityTest.cpp
@@ -44,8 +44,6 @@
 
 #include "stencil/D3Q27.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <vector>
 #include <string>
 
@@ -74,7 +72,7 @@ int main( int argc, char * argv[] )
       args.erase( vtkArgIt );
    }
    const std::string & meshFile = args[1];
-   real_t dx = boost::lexical_cast<real_t>( args[2] );
+   real_t dx = string_to_num<real_t>( args[2] );
 
    auto mesh = make_shared<mesh::TriangleMesh>();
    mesh::readAndBroadcast( meshFile, *mesh);
diff --git a/tests/mesh/MeshInitilizationTest.cpp b/tests/mesh/MeshInitilizationTest.cpp
index 494f212fb68abf7c99e626196168fa156f542092..7deeee2530e65abefb69440fe0c80e79d347e1c2 100644
--- a/tests/mesh/MeshInitilizationTest.cpp
+++ b/tests/mesh/MeshInitilizationTest.cpp
@@ -43,8 +43,6 @@
 
 #include "mesh/distance_octree/DistanceOctree.h"
 
-#include <boost/lexical_cast.hpp>
-
 #include <cmath>
 #include <vector>
 #include <string>
@@ -179,8 +177,8 @@ int main( int argc, char * argv[] )
       WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << args[0] << " MESH_FILE NUM_PROCESSES NUM_BLOCKS" );
 
    const std::string & meshFile       = args[1];
-   const uint_t        numProcesses   = boost::lexical_cast< uint_t >( args[2] );
-   const uint_t        numTotalBlocks = boost::lexical_cast< uint_t >( args[3] );
+   const uint_t        numProcesses   = string_to_num< uint_t >( args[2] );
+   const uint_t        numTotalBlocks = string_to_num< uint_t >( args[3] );
 
    test< mesh::TriangleMesh >( meshFile, numProcesses, numTotalBlocks );
    //test< mesh::FloatTriangleMesh >( meshFile, numProcesses, numTotalBlocks );
diff --git a/tests/mesh/MeshPeRaytracing.cpp b/tests/mesh/MeshPeRaytracing.cpp
index 7fecb85a888036fcf3b23d7e3bdfe39240443a1a..ce8f27992f7d8e202e797f8e716784d852214b36 100644
--- a/tests/mesh/MeshPeRaytracing.cpp
+++ b/tests/mesh/MeshPeRaytracing.cpp
@@ -72,7 +72,7 @@ int CpRayIntersectionTest(const int resolution = 10)
       for (int y = 0; y < resolution; ++y)
       {
          const real_t rand2 = real_c(y) * dx;
-         real_t theta = real_t(2) * M_PI * rand1;
+         real_t theta = real_t(2) * real_t(M_PI) * rand1;
          real_t phi = std::acos(real_t(1) - real_t(2) * rand2);
          Vec3 dir(std::sin(phi) * std::cos(theta), std::sin(phi) * std::sin(theta), std::cos(phi));
 
diff --git a/tests/mesh/PeVTKMeshWriterTest.cpp b/tests/mesh/PeVTKMeshWriterTest.cpp
index 30045036c3046de60e940d1242aed01b912c8c0a..e0a9803879ba80fedd63af3fd414eabaeb7d0105 100644
--- a/tests/mesh/PeVTKMeshWriterTest.cpp
+++ b/tests/mesh/PeVTKMeshWriterTest.cpp
@@ -95,7 +95,7 @@ std::vector<Vector3<real_t>> generatPointCloudOnSphere( const real_t radius, con
    std::vector<Vector3<real_t>> pointCloud( numPoints );
    for( auto & p : pointCloud )
    {
-      real_t theta = 2 * math::PI * distribution(rng);
+      real_t theta = 2 * real_t(M_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;
diff --git a/tests/mesh/QHullTest.cpp b/tests/mesh/QHullTest.cpp
index ddb4d46e51523e8f5608cf8b05f6c58d5d10ba25..78bd8d12001c9d8a3bbe2b133366c371cd99d917 100644
--- a/tests/mesh/QHullTest.cpp
+++ b/tests/mesh/QHullTest.cpp
@@ -209,7 +209,7 @@ std::vector<Vector3<real_t>> generatPointCloudOnSphere( const real_t radius, con
    std::vector<Vector3<real_t>> pointCloud( numPoints );
    for( auto & p : pointCloud )
    {
-      real_t theta = 2 * math::PI * distribution(rng);
+      real_t theta = 2 * real_t(M_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;
diff --git a/tests/pde/BoundaryTest.cpp b/tests/pde/BoundaryTest.cpp
index 17e6b05a3815621c94360ecc2700bd9c5be7a512..226d597b84a9266d27e879255ef6d6397e587ad1 100644
--- a/tests/pde/BoundaryTest.cpp
+++ b/tests/pde/BoundaryTest.cpp
@@ -26,6 +26,7 @@
 
 #include "core/Abort.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 
@@ -134,7 +135,7 @@ void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDa
       for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
       {
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         rhs->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+         rhs->get( *cell ) = real_t(4) * math::M_PI * math::M_PI * std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -180,7 +181,7 @@ void setBoundaryConditionsDirichl( shared_ptr< StructuredBlockForest > & blocks,
       {
 
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         real_t val = std::sin( real_t( 2 ) * math::PI * p[0] ) * std::sinh( real_t( 2 ) * math::PI * p[1] );
+         real_t val = std::sin( real_t( 2 ) * math::M_PI * p[0] ) * std::sinh( real_t( 2 ) * math::M_PI * p[1] );
 
          boundaryHandling->forceBoundary( Dirichlet_Flag, cell->x(), cell->y(), cell->z(), pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( val ) );
 
@@ -310,7 +311,7 @@ int main( int argc, char** argv )
    synchronizeD.addPackInfo( make_shared< field::communication::PackInfo< Field_T > >( dId ) );
 
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
diff --git a/tests/pde/CGTest.cpp b/tests/pde/CGTest.cpp
index eecd1c6a2646291d7509bc1f4db3a35eab9013b3..81dd0e18cb41005b5430634fb55aea5b84cfc8af 100644
--- a/tests/pde/CGTest.cpp
+++ b/tests/pde/CGTest.cpp
@@ -24,6 +24,7 @@
 
 #include "core/Abort.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 
@@ -65,7 +66,7 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
          {
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-            u->get( *cell ) = std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
          }
       }
    }
@@ -82,7 +83,7 @@ void initF( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
       for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
       {
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         f->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+         f->get( *cell ) = real_t(4) * math::M_PI * math::M_PI * std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -164,7 +165,7 @@ int main( int argc, char** argv )
    synchronizeD.addPackInfo( make_shared< field::communication::PackInfo< PdeField_T > >( dId ) );
 
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
@@ -208,4 +209,4 @@ int main( int argc, char** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/pde/JacobiTest.cpp b/tests/pde/JacobiTest.cpp
index 275cbf8dfe50506986f9d3de6d3896ad82844b76..525e0e063b8e5b8a47a08336c4ba6ecae83756d8 100644
--- a/tests/pde/JacobiTest.cpp
+++ b/tests/pde/JacobiTest.cpp
@@ -24,6 +24,7 @@
 
 #include "core/Abort.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 
@@ -69,8 +70,8 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
          {
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-            src->get( *cell ) = std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
-            dst->get( *cell ) = std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+            src->get( *cell ) = std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
+            dst->get( *cell ) = std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
          }
       }
    }
@@ -87,7 +88,7 @@ void initF( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
       for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
       {
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         f->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+         f->get( *cell ) = real_t(4) * math::M_PI * math::M_PI * std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -167,7 +168,7 @@ int main( int argc, char** argv )
    communication.addPackInfo( make_shared< field::communication::PackInfo< PdeField_T > >( srcId ) );
 
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
@@ -215,4 +216,4 @@ int main( int argc, char** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/pde/RBGSTest.cpp b/tests/pde/RBGSTest.cpp
index c3a922ea15176f9fc6ba0674601a683a393807c8..8609ace0af68ecbfeca6b75c22ca8aed39e3ded8 100644
--- a/tests/pde/RBGSTest.cpp
+++ b/tests/pde/RBGSTest.cpp
@@ -24,6 +24,7 @@
 
 #include "core/Abort.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 
@@ -68,7 +69,7 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
          {
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-            u->get( *cell ) = std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
          }
       }
    }
@@ -85,7 +86,7 @@ void initF( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
       for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
       {
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         f->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+         f->get( *cell ) = real_t(4) * math::M_PI * math::M_PI * std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -164,7 +165,7 @@ int main( int argc, char** argv )
    communication.addPackInfo( make_shared< field::communication::PackInfo< PdeField_T > >( uId ) );
 
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
@@ -215,4 +216,4 @@ int main( int argc, char** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/pde/SORTest.cpp b/tests/pde/SORTest.cpp
index 520df1cfa8fd8eb1796f20710d640dae962d0e25..0a5799163292f662469395294bde526e82706dde 100644
--- a/tests/pde/SORTest.cpp
+++ b/tests/pde/SORTest.cpp
@@ -24,6 +24,7 @@
 
 #include "core/Abort.h"
 #include "core/debug/TestSubsystem.h"
+#include "core/math/Constants.h"
 #include "core/mpi/Environment.h"
 #include "core/mpi/MPIManager.h"
 
@@ -68,7 +69,7 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
          for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
          {
             const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-            u->get( *cell ) = std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
          }
       }
    }
@@ -85,7 +86,7 @@ void initF( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
       for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
       {
          const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
-         f->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+         f->get( *cell ) = real_t(4) * math::M_PI * math::M_PI * std::sin( real_t(2) * math::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
       }
    }
 }
@@ -166,7 +167,7 @@ int main( int argc, char** argv )
    communication.addPackInfo( make_shared< field::communication::PackInfo< PdeField_T > >( uId ) );
 
    std::vector< real_t > weights( Stencil_T::Size );
-   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::M_PI * math::M_PI;
    weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
    weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
@@ -217,4 +218,4 @@ int main( int argc, char** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/pe/Collision.cpp b/tests/pe/Collision.cpp
index 958b4de97bf52072a8f70c723b14e694f72e512c..e73d9ac3522cd3cd88e9cfa9f8d08e372d7c598f 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -119,8 +119,8 @@ void BoxTest()
    b4.rotate( Vec3(1,1,0), real_t(atan(sqrt(2))) );
 
    Box b5(123, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
-   b5.rotate( Vec3(0,0,1), real_t(math::PI * 0.25) );
-   b5.rotate( Vec3(1,0,0), real_t(math::PI * 0.25) );
+   b5.rotate( Vec3(0,0,1), real_t(math::M_PI * 0.25) );
+   b5.rotate( Vec3(1,0,0), real_t(math::M_PI * 0.25) );
 
    std::vector<Contact> contacts;
    fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
diff --git a/tests/pe/GJK_EPA.cpp b/tests/pe/GJK_EPA.cpp
index 931d14fd901d97a94059ee8397880b4d5e05ff5a..9293c3c866c440f89d06dde17fae3b8dcb14a525 100644
--- a/tests/pe/GJK_EPA.cpp
+++ b/tests/pe/GJK_EPA.cpp
@@ -128,11 +128,11 @@ int main( int argc, char ** argv )
 
 //      WALBERLA_LOG_INFO("**** (" << dir.cx() << ", " << dir.cy() << ", " << dir.cz() << ") ****");
 //      WALBERLA_LOG_INFO("deltaPoint : |" << contactPoint - contactPointGJK << "| = " << (contactPoint - contactPointGJK).length() );
-//      WALBERLA_LOG_INFO("deltaNormal: |" << normal - normalGJK << "| = " << (normal - normalGJK).length() << " (" << acos(normal*normalGJK) / math::PI * 180 << "°)" );
+//      WALBERLA_LOG_INFO("deltaNormal: |" << normal - normalGJK << "| = " << (normal - normalGJK).length() << " (" << acos(normal*normalGJK) / math::M_PI * 180 << "°)" );
 //      WALBERLA_LOG_INFO("deltaPen   : " << penetrationDepth << " - " << penetrationDepthGJK << " = " << penetrationDepth - penetrationDepthGJK);
 
       point.insert( (contactPoint - contactPointGJK).length() );
-      norm.insert( acos(normal*normalGJK) / math::PI * 180 );
+      norm.insert( acos(normal*normalGJK) / math::M_PI * 180 );
       pen.insert( fabs(penetrationDepth - penetrationDepthGJK) );
    }
 
diff --git a/tests/pe/PeDocumentationSnippets.cpp b/tests/pe/PeDocumentationSnippets.cpp
index d9b3396bbb6125b44133861571b632fb30d642e0..e5067ab1056e12f2c8bcd0daf75b451ef4b2439b 100644
--- a/tests/pe/PeDocumentationSnippets.cpp
+++ b/tests/pe/PeDocumentationSnippets.cpp
@@ -98,14 +98,14 @@ int main( int argc, char ** argv )
    // be used to for instance rotate the box around the global y-axis.
    BoxID box = createBox( *globalBodyStorage, forest->getBlockStorage(), storageID, 1, Vec3(2,3,4), Vec3(2.5,2.5,2.5) );
    if (box != nullptr)
-      box->rotate( 0.0, real_c(math::PI/3.0), 0.0 );
+      box->rotate( 0.0, real_c(math::M_PI/3.0), 0.0 );
    //! [Create a Box]
 
    //! [Create a Capsule]
    // Create a capsule and rotate it after successfull creation.
    CapsuleID capsule = createCapsule( *globalBodyStorage, forest->getBlockStorage(), storageID, 1, Vec3(2,3,4), real_t(1), real_t(1) );
    if (capsule != nullptr)
-      capsule->rotate( 0.0, real_c(math::PI/3.0), 0.0 );
+      capsule->rotate( 0.0, real_c(math::M_PI/3.0), 0.0 );
    //! [Create a Capsule]
 
    //! [Create a Plane]
@@ -118,7 +118,7 @@ int main( int argc, char ** argv )
    // Create a sphere and rotate it after successfull creation.
    SphereID sphere = createSphere( *globalBodyStorage, forest->getBlockStorage(), storageID, 1, Vec3(2,3,4), real_t(1) );
    if (sphere != nullptr)
-      sphere->rotate( 0.0, real_c(math::PI/3.0), 0.0 );
+      sphere->rotate( 0.0, real_c(math::M_PI/3.0), 0.0 );
    //! [Create a Sphere]
 
    //! [Create a Union]
diff --git a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
index 1a4293fec9670eabbe9b0669f24e37b61d674ee3..f3dc7905e01378bffc4e39dea83ceb059758006e 100644
--- a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
@@ -738,7 +738,7 @@ int main( int argc, char **argv )
    const real_t viscosity = ( viscosity_SI/densityFluid_SI ) * dt_SI / ( dx_SI * dx_SI );
    const real_t tau = real_t(1) / lbm::collision_model::omegaFromViscosity( viscosity );
 
-   real_t gravitationalForce = - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::PI / real_t(6);
+   real_t gravitationalForce = - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::M_PI / real_t(6);
 
    // unhindered settling velocity of a single sphere in infinite fluid, would come from experiments or DNS, NOT Stokes settling velocity, from Finn et al, Tab 5
    const real_t velUnhindered_SI = real_t(-0.048); // m/s
@@ -1488,9 +1488,9 @@ int main( int argc, char **argv )
                WALBERLA_LOG_INFO_ON_ROOT("initial simulation ended with relative difference of interaction forces of " << relativeForceDiffLimit << " after " << t << " time steps.");
 
 
-               real_t actingExternalForceOnSpheres = real_c(numSpheres) * ( ( - gravity * densityRatio * diameter * diameter * diameter * math::PI / real_t(6)  ) +
-                                                                            ( gravity * real_t(1) * diameter * diameter * diameter * math::PI / real_t(6) ) +
-                                                                            ( extForce[2] * real_t(1) * diameter * diameter * diameter * math::PI / real_t(6) ) );
+               real_t actingExternalForceOnSpheres = real_c(numSpheres) * ( ( - gravity * densityRatio * diameter * diameter * diameter * math::M_PI / real_t(6)  ) +
+                                                                            ( gravity * real_t(1) * diameter * diameter * diameter * math::M_PI / real_t(6) ) +
+                                                                            ( extForce[2] * real_t(1) * diameter * diameter * diameter * math::M_PI / real_t(6) ) );
                WALBERLA_LOG_INFO_ON_ROOT("f_interaction_z = " << curInteractionForce << ", f_ext_z = " << actingExternalForceOnSpheres );
                if( std::fabs( ( std::fabs( curInteractionForce ) - std::fabs( actingExternalForceOnSpheres ) )/ std::fabs( curInteractionForce ) ) < relativeForceConvergenceLimit )
                {
@@ -1521,7 +1521,7 @@ int main( int argc, char **argv )
    WALBERLA_LOG_INFO_ON_ROOT("===================================================================================" );
    WALBERLA_LOG_INFO_ON_ROOT("Starting simulation with:" );
    WALBERLA_LOG_INFO_ON_ROOT("external forcing on fluid = " << extForce );
-   WALBERLA_LOG_INFO_ON_ROOT("total external forces on all particles = " << real_c(numSpheres) * ( - gravity * ( densityRatio - real_t(1) ) + extForce[2] ) * diameter * diameter * diameter * math::PI / real_t(6) );
+   WALBERLA_LOG_INFO_ON_ROOT("total external forces on all particles = " << real_c(numSpheres) * ( - gravity * ( densityRatio - real_t(1) ) + extForce[2] ) * diameter * diameter * diameter * math::M_PI / real_t(6) );
    WALBERLA_LOG_INFO_ON_ROOT("simulating " << timesteps << " time steps" );
 
 
@@ -1642,9 +1642,9 @@ int main( int argc, char **argv )
 
       // ext forces on bodies
       timeloop.add() << Sweep( DummySweep(), "Dummy Sweep ")
-                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,- gravity * densityRatio * diameter * diameter * diameter * math::PI / real_t(6) )  ), "Gravitational Force Add" )
-                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,gravity * real_t(1) * diameter * diameter * diameter * math::PI / real_t(6) ) ), "Buoyancy Force (due to gravity) Add" )
-                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,extForce[2] * real_t(1) * diameter * diameter * diameter * math::PI / real_t(6) ) ), "Buoyancy Force (due to external fluid force) Add" )
+                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,- gravity * densityRatio * diameter * diameter * diameter * math::M_PI / real_t(6) )  ), "Gravitational Force Add" )
+                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,gravity * real_t(1) * diameter * diameter * diameter * math::M_PI / real_t(6) ) ), "Buoyancy Force (due to gravity) Add" )
+                     << AfterFunction( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, Vector3<real_t>(0,0,extForce[2] * real_t(1) * diameter * diameter * diameter * math::M_PI / real_t(6) ) ), "Buoyancy Force (due to external fluid force) Add" )
                      << AfterFunction( pe_coupling::TimeStep( blocks, bodyStorageID, *cr, syncCall, dtInteractionSubCycle, peSubSteps, lubricationEvaluationFunction ), "Pe Time Step" );
 
       timeloop.add() << Sweep( DummySweep(), "Dummy Sweep ")
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index 2e449157fc7bf338a321fcb794e917d63e227e84..0c717fb156cb1075a37813f778fad4d57d96d352 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -226,9 +226,9 @@ public:
       real_t uBar = computeAverageVel();
 
       // f_total = f_drag + f_buoyancy
-      real_t totalForce = forceX  + real_c(4.0/3.0) * math::PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
 
       // update drag force values
       normalizedDragOld_ = normalizedDragNew_;
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
index 358d165bae76cd446dc8709fd7059b2eb9f6a0ae..1f61e239cc5dd21566990d042aedc9573c14d214 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -297,9 +297,9 @@ class ForceEval
          real_t uBar = getAverageVel();
 
          // f_total = f_drag + f_buoyancy
-         real_t totalForce = forceX  + real_c(4.0/3.0) * math::PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+         real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-         real_t normalizedDragForce = totalForce / real_c( 6.0 * math::PI * setup_->visc * setup_->radius * uBar );
+         real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
 
          // update drag force values
          normalizedDragOld_ = normalizedDragNew_;
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 2ebe21f515c607ccc96698d66d137dc0ce3f2b4e..f1b545ec0a58507829da881819909f000bf37467 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -261,10 +261,10 @@ private:
             WALBERLA_CHECK_FLOAT_EQUAL ( forceSphr2[0], -forceSphr1[0] );
 
             // according to the formula from Ding & Aidun 2003
-            // F = 3/2 * PI * rho_L * nu_L * relative velocity of both spheres * r * r * 1/gap
+            // F = 3/2 * M_PI * rho_L * nu_L * relative velocity of both spheres * r * r * 1/gap
             // the correct analytically calculated value is 339.2920063998
             // in this geometry setup the relative error is 0.1246489711 %
-            real_t analytical = real_c(3.0)/real_c(2.0) * walberla::math::PI * real_c(1.0) * nu_L_ * real_c(2.0) * real_c(vel_[0]) * radius_ * radius_ * real_c(1.0)/gap;
+            real_t analytical = real_c(3.0)/real_c(2.0) * walberla::math::M_PI * real_c(1.0) * nu_L_ * real_c(2.0) * real_c(vel_[0]) * radius_ * radius_ * real_c(1.0)/gap;
             real_t relErr     = std::fabs( analytical - forceSphr2[0] ) / analytical * real_c(100.0);
             WALBERLA_CHECK_LESS( relErr, real_t(1) );
          }
@@ -357,10 +357,10 @@ private:
          if ( timestep == uint_t(26399) )
          {
             // according to the formula from Ding & Aidun 2003
-            // F = 6 * PI * rho_L * nu_L * relative velocity of both bodies=relative velocity of the sphere * r * r * 1/gap
+            // F = 6 * M_PI * rho_L * nu_L * relative velocity of both bodies=relative velocity of the sphere * r * r * 1/gap
             // the correct analytically calculated value is 339.292006996217
             // in this geometry setup the relative error is 0.183515322065561 %
-            real_t analytical = real_c(6.0) * walberla::math::PI * real_c(1.0) * nu_L_ * real_c(-vel_[0]) * radius_ * radius_ * real_c(1.0)/gap;
+            real_t analytical = real_c(6.0) * walberla::math::M_PI * real_c(1.0) * nu_L_ * real_c(-vel_[0]) * radius_ * radius_ * real_c(1.0)/gap;
             real_t relErr     = std::fabs( analytical - forceSphr1[0] ) / analytical * real_c(100.0);
             WALBERLA_CHECK_LESS( relErr, real_t(1) );
          }
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index 68a74cde6dda8d148e7928429a5ee9b4f8dfe168..dd90350c7e9ed3075fafc438896c09a816492902 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -375,7 +375,7 @@ int main( int argc, char **argv )
    setup.checkFrequency = uint_t( 100 );                      // evaluate the torque only every checkFrequency time steps
    setup.radius = real_c(0.5) * chi * real_c( setup.length ); // sphere radius
    setup.visc   = ( setup.tau - real_c(0.5) ) / real_c(3);    // viscosity in lattice units
-   setup.phi    = real_c(4.0/3.0) * math::PI * setup.radius * setup.radius * setup.radius /
+   setup.phi    = real_c(4.0/3.0) * math::M_PI * setup.radius * setup.radius * setup.radius /
                   ( real_c( setup.length * setup.length * setup.length ) ); // solid volume fraction
    const real_t omega      = real_c(1) / setup.tau;           // relaxation rate
    const real_t dx         = real_c(1);                       // lattice dx
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
index b5e53a23fc1856a18548a4929634e94c5b320179..1fadccab92c4c6c88aea97fe1bc902738e968459 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
@@ -178,9 +178,9 @@ public:
       real_t uBar = computeAverageVel();
 
       // f_total = f_drag + f_buoyancy
-      real_t totalForce = forceX  + real_c(4.0/3.0) * math::PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
 
       // update drag force values
       normalizedDragOld_ = normalizedDragNew_;
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
index b0ca02da90f9b912399af654658655dd8f8219b1..bb890e698f560b4d5b6a7c1445c7ed1b45ce5dc2 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
@@ -297,9 +297,9 @@ public:
       real_t uBar = computeAverageVel();
 
       // f_total = f_drag + f_buoyancy
-      real_t totalForce = forceX  + real_c(4.0/3.0) * math::PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
 
       // update drag force values
       normalizedDragOld_ = normalizedDragNew_;
diff --git a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
index 7cd33d47ed694d99b08525e161837bb85b4b8b7a..c63f974e827b1b847716ee05c76a9590c6c89186 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
@@ -309,7 +309,7 @@ int main( int argc, char **argv )
    setup.checkFrequency = uint_t( 100 );                      // evaluate the torque only every checkFrequency time steps
    setup.radius = real_c(0.5) * chi * real_c( setup.length ); // sphere radius
    setup.visc   = ( setup.tau - real_c(0.5) ) / real_c(3);    // viscosity in lattice units
-   setup.phi    = real_c(4.0/3.0) * math::PI * setup.radius * setup.radius * setup.radius /
+   setup.phi    = real_c(4.0/3.0) * math::M_PI * setup.radius * setup.radius * setup.radius /
                   ( real_c( setup.length * setup.length * setup.length ) ); // solid volume fraction
    const real_t omega      = real_c(1) / setup.tau;           // relaxation rate
    const real_t dx         = real_c(1);                       // lattice dx