diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
index e7a4a0fead15db23016d18cccd0254b788d1594c..935d226ebca87cf1e064e203c12a7456366f7db6 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
@@ -1229,7 +1229,7 @@ int main( int argc, char **argv )
 
    const Vector3<uint_t> domainSize( XBlocks * blockSize, YBlocks * blockSize, ZBlocks * blockSize );
    const auto domainVolume = real_t(domainSize[0] * domainSize[1] * domainSize[2]);
-   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+   const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
    const uint_t numberOfSediments = uint_c(std::ceil(solidVolumeFraction * domainVolume / sphereVolume));
 
    real_t expectedSedimentVolumeFraction = (useBox||useHopper) ? real_t(0.45) : real_t(0.52);
diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
index 7b87d670a294f757cea8750597a8a34bd52b36ff..9e03ef1a9089495c0287af8343df40f9ee9bb5c2 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
@@ -463,8 +463,8 @@ int main( int argc, char **argv )
    real_t domainHeight = real_c(ZCells) - topWallOffset;
    real_t fluidVolume =  real_c( XCells * YCells ) * domainHeight;
    real_t solidVolume = solidVolumeFraction * fluidVolume;
-   uint_t numberOfParticles = uint_c(std::ceil(solidVolume / ( math::M_PI / real_t(6) * diameter * diameter * diameter )));
-   diameter = std::cbrt( solidVolume / ( real_c(numberOfParticles) * math::M_PI / real_t(6) ) );
+   uint_t numberOfParticles = uint_c(std::ceil(solidVolume / ( math::pi / real_t(6) * diameter * diameter * diameter )));
+   diameter = std::cbrt( solidVolume / ( real_c(numberOfParticles) * math::pi / real_t(6) ) );
 
    auto densityRatio = real_t(2.5);
 
@@ -821,7 +821,7 @@ int main( int argc, char **argv )
 
    }
 
-   real_t sphereVolume = diameter * diameter * diameter * math::M_PI / real_t(6);
+   real_t sphereVolume = diameter * diameter * diameter * math::pi / real_t(6);
    Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densityRatio - real_t(1)) * gravitationalAcceleration * sphereVolume );
    timeloop.addFuncAfterTimeStep(pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, gravitationalForce ), "Gravitational force" );
 
diff --git a/apps/benchmarks/DEM/DEM.cpp b/apps/benchmarks/DEM/DEM.cpp
index d3a0decf71375e78b2fa109d888f87a68e6e20c7..a76f4993f8514030d92c52e6d9cb09573f7ffa8c 100644
--- a/apps/benchmarks/DEM/DEM.cpp
+++ b/apps/benchmarks/DEM/DEM.cpp
@@ -33,13 +33,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::M_PI / std::sqrt(k / meff - a*a));
+   return std::exp(-a * math::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::M_PI / std::sqrt( k/meff - a*a);
+   return math::pi / std::sqrt( k/meff - a*a);
 }
 }
 
diff --git a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
index b78fd4395482d64c340b7acc27011a432fc898aa..7e5c37e2e6d0f56a765b987f774bb721df67925f 100644
--- a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
+++ b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
@@ -675,7 +675,7 @@ int main( int argc, char **argv )
    auto refinementTimestep = lbm::refinement::makeTimeStep< LatticeModel_T, BoundaryHandling_T >( blocks, sweep, pdfFieldID, boundaryHandlingID );
 
    // add force evaluation and logging
-   real_t normalizationFactor = ( zeroShearTest ) ? real_t(1) : ( math::M_PI / real_t(8) * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter );
+   real_t normalizationFactor = ( zeroShearTest ) ? real_t(1) : ( math::pi / real_t(8) * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter );
    std::string loggingFileName( baseFolderLogging + "/LoggingForcesNearPlane");
    loggingFileName += "_lvl" + std::to_string(numberOfLevels);
    loggingFileName += "_D" + std::to_string(uint_c(diameter));
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index 725ff77b5e3894255c570f8d1f1465e7b1562351..9c7b016f0d059a1f471df256c884b04055fe9e16 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -1080,7 +1080,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::M_PI / real_t(6) );
+      gravity = forceEval->getForce() / ( (densityRatio - real_t(1) ) * diameter * diameter * diameter * math::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 );
@@ -1234,7 +1234,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::M_PI / real_t(6));
+   Vector3<real_t> extForcesOnSphere( real_t(0), real_t(0), - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::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 ff9d0f46beaaeb6c567244fa0e738ca01d4c72ad..1a0482f5781a272ae6159653f4e39435556f2aa7 100644
--- a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
+++ b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
@@ -762,7 +762,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::M_PI * setup.radius_L * setup.radius_L;
+      setup.flowRate_L = setup.meanVelocity_L * math::pi * setup.radius_L * setup.radius_L;
    }
    else
    {
diff --git a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
index 2f7a5f6b5858cf51664abc1d0938abf590aab18d..fcd8ba66932c1433b8c1ed4f7e11933982f277c3 100644
--- a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
+++ b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
@@ -744,7 +744,7 @@ public:
 
    void operator()( const real_t t )
    {
-      tConstTerm_ = ( sinPeriod_ > real_t(0) ) ? ( std::abs( std::sin( math::M_PI * t / sinPeriod_ ) ) ) : real_t(1);
+      tConstTerm_ = ( sinPeriod_ > real_t(0) ) ? ( std::abs( std::sin( math::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/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp b/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
index b91418447e3a8a5b0e7bcf503aa1bfbe98c4f580..d3e6fe63a866c8df9d30a9c5ea23d8522d868a9f 100644
--- a/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
+++ b/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
@@ -371,7 +371,7 @@ uint_t createSpheresRandomly( StructuredBlockForest & forest, pe::BodyStorage &
 
       pe::createSphere( globalBodyStorage, forest.getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( xParticle, yParticle, zParticle ), creationDiameter * real_t(0.5), material );
 
-      currentSphereVolume += math::M_PI / real_t(6) * creationDiameter * creationDiameter * creationDiameter;
+      currentSphereVolume += math::pi / real_t(6) * creationDiameter * creationDiameter * creationDiameter;
 
       ++numberOfSpheres;
    }
@@ -851,16 +851,16 @@ int main( int argc, char **argv ) {
    const real_t restitutionCoeff = real_t(0.88);
    const real_t frictionCoeff = real_t(0.25);
 
-   real_t sphereVolume = diameterAvg * diameterAvg * diameterAvg * math::M_PI / real_t(6); // based on avg. diameter
+   real_t sphereVolume = diameterAvg * diameterAvg * diameterAvg * math::pi / real_t(6); // based on avg. diameter
    const real_t particleMass = densityRatio * sphereVolume;
    const real_t Mij = particleMass * particleMass / (real_t(2) * particleMass);
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
    const real_t collisionTime = real_t(0.5);
-   const real_t stiffnessCoeff = math::M_PI * math::M_PI * Mij / (collisionTime * collisionTime *
-                                 (real_t(1) - lnDryResCoeff * lnDryResCoeff / (math::M_PI * math::M_PI + lnDryResCoeff * lnDryResCoeff)));
+   const real_t stiffnessCoeff = math::pi * math::pi * Mij / (collisionTime * collisionTime *
+                                 (real_t(1) - lnDryResCoeff * lnDryResCoeff / (math::pi * math::pi + lnDryResCoeff * lnDryResCoeff)));
    const real_t dampingCoeff = -real_t(2) * std::sqrt(Mij * stiffnessCoeff) *
-                               (std::log(restitutionCoeff) / std::sqrt(math::M_PI * math::M_PI + (std::log(restitutionCoeff) * std::log(restitutionCoeff))));
-   const real_t contactDuration = real_t(2) * math::M_PI * Mij / (std::sqrt(real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff)); //formula from Uhlman
+                               (std::log(restitutionCoeff) / std::sqrt(math::pi * math::pi + (std::log(restitutionCoeff) * std::log(restitutionCoeff))));
+   const real_t contactDuration = real_t(2) * math::pi * Mij / (std::sqrt(real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff)); //formula from Uhlman
 
    WALBERLA_LOG_INFO_ON_ROOT("Created particle material with:\n"
                                    << " - coefficient of restitution = " << restitutionCoeff << "\n"
diff --git a/apps/tutorials/pde/01_SolvingPDE.cpp b/apps/tutorials/pde/01_SolvingPDE.cpp
index 11b47c0f0c989122816afa4060f7d69aab306856..4b1e7d92f5edfaf572085862668ea085e6c22d96 100644
--- a/apps/tutorials/pde/01_SolvingPDE.cpp
+++ b/apps/tutorials/pde/01_SolvingPDE.cpp
@@ -70,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*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 );
+            // 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 );
          }
       }
    }
@@ -98,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*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] );
+         // 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] );
       }
    }
 }
@@ -147,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::M_PI * math::M_PI );
+      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::pi * math::pi );
    )
 
    // swap source and destination fields
@@ -282,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::M_PI * math::M_PI;
+   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::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 8a7e33dce7da53943f4fbdffd4a5f387bfd1165f..a978dbf0cf294c20b92903e495a5de833b4ab735 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::M_PI * math::M_PI );
+      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::pi * math::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::M_PI * math::M_PI;
+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::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*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] );
+         // 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] );
       }
    }
 }
@@ -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*M_PI*x)*sinh(2*M_PI)
+            // 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::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 );
+            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 );
          }
       }
    }
diff --git a/apps/tutorials/pde/02_HeatEquation.cpp b/apps/tutorials/pde/02_HeatEquation.cpp
index c4ec151ea024d15505e13429e64c82e642111700..dea4433808530f6d64fb25e5d881b214ca09ce6d 100644
--- a/apps/tutorials/pde/02_HeatEquation.cpp
+++ b/apps/tutorials/pde/02_HeatEquation.cpp
@@ -66,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(M_PI*x)*sin(M_PI*y)
-         u->get( *cell ) = std::sin( math::M_PI * p[0] ) * std::sin( math::M_PI * p[1] );
+         // 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] );
       }
    }
 }
diff --git a/apps/tutorials/pde/03_HeatEquation_Extensions.cpp b/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
index 92f3b79acf19310c2baa6e5c192ad1409dfc3248..9f65e2b125fddd3b50ddfb468cbeb05c07f342b0 100644
--- a/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
+++ b/apps/tutorials/pde/03_HeatEquation_Extensions.cpp
@@ -67,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(M_PI*x)*sin(M_PI*y)
-         u->get( *cell ) = std::sin( math::M_PI * p[0] ) * std::sin( math::M_PI * p[1] );
+         // 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] );
       }
    }
 }
diff --git a/python/mesa_pd/templates/kernel/SpringDashpot.templ.h b/python/mesa_pd/templates/kernel/SpringDashpot.templ.h
index 1a6f39ca1239a949f9712017668386d4192b3aac..93df9f865d5ca019e0c62c189bc498994f1e981e 100644
--- a/python/mesa_pd/templates/kernel/SpringDashpot.templ.h
+++ b/python/mesa_pd/templates/kernel/SpringDashpot.templ.h
@@ -90,7 +90,7 @@ public:
                                        const real_t meff)
    {
       auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
-      return std::exp(-a * math::M_PI / std::sqrt(getStiffness(type1, type2) / meff - a*a));
+      return std::exp(-a * math::pi / std::sqrt(getStiffness(type1, type2) / meff - a*a));
    }
 
    inline
@@ -99,7 +99,7 @@ public:
                             const real_t meff)
    {
       auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
-      return math::M_PI / std::sqrt( getStiffness(type1, type2)/meff - a*a);
+      return math::pi / std::sqrt( getStiffness(type1, type2)/meff - a*a);
    }
 
    inline
@@ -110,8 +110,8 @@ public:
                              const real_t meff)
    {
       const real_t lnDryResCoeff = std::log(cor);
-      setStiffness(type1, type2, math::M_PI * math::M_PI * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  ));
-      setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffness(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) ));
+      setStiffness(type1, type2, math::pi * math::pi * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  ));
+      setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffness(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) ));
    }
 private:
    uint_t numParticleTypes_;
diff --git a/src/core/math/Constants.h b/src/core/math/Constants.h
index b02700a39f59928d6ac87ec051d671bed44973b5..79e97e97ed716524712878b36566379a01d36c24 100644
--- a/src/core/math/Constants.h
+++ b/src/core/math/Constants.h
@@ -16,57 +16,43 @@
 //! \file Constants.h
 //! \author Klaus Iglberger
 //! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
 //! \brief Header file for mathematical constants
 //
 //======================================================================================================================
 
 #pragma once
 
-
 //*************************************************************************************************
 // Includes
 //*************************************************************************************************
 
-#include <core/DataTypes.h>
-
 #include <cmath>
+#include <core/DataTypes.h>
 
 // Disable false warnings in GCC 5
-#if ( defined __GNUC__ ) && ( __GNUC__ == 5 ) && ( __GNUC_MINOR__ == 1 )
-#  pragma GCC diagnostic push
-#  pragma GCC diagnostic ignored "-Wunused-variable"
+#if (defined __GNUC__) && (__GNUC__ == 5) && (__GNUC_MINOR__ == 1)
+#   pragma GCC diagnostic push
+#   pragma GCC diagnostic ignored "-Wunused-variable"
 #endif
 
-namespace walberla {
-namespace math {
-
-# undef M_E
-# undef M_LOG2E
-# undef M_LOG10E
-# undef M_LN2
-# undef M_LN10
-# undef M_PI
-# undef M_PI_2
-# undef M_PI_4
-# undef M_1_PI
-# undef M_2_PI
-# undef M_2_SQRTPI
-# undef M_SQRT2
-# undef M_SQRT1_2
-
-constexpr real_t M_E          = real_t( 2.7182818284590452354  );  /* e */
-constexpr real_t M_LOG2E      = real_t( 1.4426950408889634074  );  /* log_2 e */
-constexpr real_t M_LOG10E     = real_t( 0.43429448190325182765 );  /* log_10 e */
-constexpr real_t M_LN2        = real_t( 0.69314718055994530942 );  /* log_e 2 */
-constexpr real_t M_LN10       = real_t( 2.30258509299404568402 );  /* log_e 10 */
-constexpr real_t M_PI         = real_t( 3.14159265358979323846 );  /* pi */
-constexpr real_t M_PI_2       = real_t( 1.57079632679489661923 );  /* pi/2 */
-constexpr real_t M_PI_4       = real_t( 0.78539816339744830962 );  /* pi/4 */
-constexpr real_t M_1_PI       = real_t( 0.31830988618379067154 );  /* 1/pi */
-constexpr real_t M_2_PI       = real_t( 0.63661977236758134308 );  /* 2/pi */
-constexpr real_t M_2_SQRTPI   = real_t( 1.12837916709551257390 );  /* 2/sqrt(pi) */
-constexpr real_t M_SQRT2      = real_t( 1.41421356237309504880 );  /* sqrt(2) */
-constexpr real_t M_SQRT1_2    = real_t( 0.70710678118654752440 );  /* 1/sqrt(2) */
+namespace walberla
+{
+namespace math
+{
+constexpr real_t e                = real_t(2.7182818284590452354);  /* e */
+constexpr real_t log2_e           = real_t(1.4426950408889634074);  /* log_2 e */
+constexpr real_t log10_e          = real_t(0.43429448190325182765); /* log_10 e */
+constexpr real_t ln_two           = real_t(0.69314718055994530942); /* log_e 2 */
+constexpr real_t ln_ten           = real_t(2.30258509299404568402); /* log_e 10 */
+constexpr real_t pi               = real_t(3.14159265358979323846); /* pi */
+constexpr real_t half_pi          = real_t(1.57079632679489661923); /* pi/2 */
+constexpr real_t fourth_pi        = real_t(0.78539816339744830962); /* pi/4 */
+constexpr real_t one_div_pi       = real_t(0.31830988618379067154); /* 1/pi */
+constexpr real_t two_div_pi       = real_t(0.63661977236758134308); /* 2/pi */
+constexpr real_t two_div_root_pi  = real_t(1.12837916709551257390); /* 2/sqrt(pi) */
+constexpr real_t root_two         = real_t(1.41421356237309504880); /* sqrt(2) */
+constexpr real_t one_div_root_two = real_t(0.70710678118654752440); /* 1/sqrt(2) */
 
 } // namespace math
-}
+} // namespace walberla
\ No newline at end of file
diff --git a/src/core/math/GenericAABB.h b/src/core/math/GenericAABB.h
index 9b321e738f53548dba55094722425937967bef12..a6e18882af3e2df37e7f9fb00983d360ed3d6d14 100644
--- a/src/core/math/GenericAABB.h
+++ b/src/core/math/GenericAABB.h
@@ -111,8 +111,8 @@ public:
    inline bool containsClosedInterval( const vector_type & point ) const;
    inline bool containsClosedInterval( const vector_type & point, const value_type dx ) const;
 
-   inline GenericAABB getExtended( const value_type e ) const;
-   inline GenericAABB getExtended( const vector_type & e ) const;
+   inline GenericAABB getExtended( const value_type eps ) const;
+   inline GenericAABB getExtended( const vector_type & eps ) const;
 
    inline GenericAABB getTranslated( const vector_type & translation ) const;
 
@@ -162,8 +162,8 @@ public:
 
    inline void setAxisBounds( const uint_t index, const value_type value1, const value_type value2 );
 
-   inline void extend( const value_type e );
-   inline void extend( const vector_type & e );
+   inline void extend( const value_type eps );
+   inline void extend( const vector_type & eps );
 
    inline void setCenter( const vector_type & center );
    inline void translate( const vector_type & translation );
diff --git a/src/core/math/GenericAABB.impl.h b/src/core/math/GenericAABB.impl.h
index 1fae800a88c6f7bcfb1d7a2b24cc9a30411c8066..fd1d39996bbb110060cb5b57a48fdbb241f50aa9 100644
--- a/src/core/math/GenericAABB.impl.h
+++ b/src/core/math/GenericAABB.impl.h
@@ -674,19 +674,19 @@ bool GenericAABB< T >::containsClosedInterval( const vector_type & point, const
 /**
  * \brief Creates a new GenericAABB by extending this one
  *
- * \param e epsilon by which the bounding box is extended in each direction
+ * \param eps epsilon by which the bounding box is extended in each direction
  *
  * \returns The extended GenericAABB
  */
 template< typename T >
-GenericAABB< T > GenericAABB< T >::getExtended( const value_type e ) const
+GenericAABB< T > GenericAABB< T >::getExtended( const value_type eps ) const
 {
-   vector_type newMinCorner( minCorner_[0] - e, minCorner_[1] - e, minCorner_[2] - e );
+   vector_type newMinCorner( minCorner_[0] - eps, minCorner_[1] - eps, minCorner_[2] - eps );
 
    return createFromMinMaxCorner( newMinCorner[0], newMinCorner[1], newMinCorner[2],
-                                  std::max( newMinCorner[0], maxCorner_[0] + e ),
-                                  std::max( newMinCorner[1], maxCorner_[1] + e ),
-                                  std::max( newMinCorner[2], maxCorner_[2] + e ) );
+                                  std::max( newMinCorner[0], maxCorner_[0] + eps ),
+                                  std::max( newMinCorner[1], maxCorner_[1] + eps ),
+                                  std::max( newMinCorner[2], maxCorner_[2] + eps ) );
 }
 
 
@@ -694,19 +694,19 @@ GenericAABB< T > GenericAABB< T >::getExtended( const value_type e ) const
 /**
  * \brief Creates a new GenericAABB by extending this one
  *
- * \param e epsilon vector by which the bounding box is extended. The box is extended in each direction by
+ * \param eps epsilon vector by which the bounding box is extended. The box is extended in each direction by
  *          the corresponding vector component.
  *
  * \returns The extended GenericAABB
  */
 template< typename T >
-GenericAABB< T > GenericAABB< T >::getExtended( const vector_type & e ) const
+GenericAABB< T > GenericAABB< T >::getExtended( const vector_type & eps ) const
 {
-   vector_type newMinCorner( minCorner_ - e );
+   vector_type newMinCorner( minCorner_ - eps );
    return createFromMinMaxCorner( newMinCorner[0], newMinCorner[1], newMinCorner[2],
-                                  std::max( newMinCorner[0], maxCorner_[0] + e[0] ),
-                                  std::max( newMinCorner[1], maxCorner_[1] + e[0] ),
-                                  std::max( newMinCorner[2], maxCorner_[2] + e[0] ) );
+                                  std::max( newMinCorner[0], maxCorner_[0] + eps[0] ),
+                                  std::max( newMinCorner[1], maxCorner_[1] + eps[0] ),
+                                  std::max( newMinCorner[2], maxCorner_[2] + eps[0] ) );
 }
 
 
@@ -1519,18 +1519,18 @@ void GenericAABB< T >::setAxisBounds( const uint_t index, const value_type value
 /**
  * \brief Extends this GenericAABB
  *
- * \param e epsilon by which the bounding box is extended in each direction
+ * \param eps epsilon by which the bounding box is extended in each direction
  */
 template< typename T >
-void GenericAABB< T >::extend( const value_type e )
+void GenericAABB< T >::extend( const value_type eps )
 {
-   minCorner_[0] -= e;
-   minCorner_[1] -= e;
-   minCorner_[2] -= e;
+   minCorner_[0] -= eps;
+   minCorner_[1] -= eps;
+   minCorner_[2] -= eps;
 
-   maxCorner_[0] += e;
-   maxCorner_[1] += e;
-   maxCorner_[2] += e;
+   maxCorner_[0] += eps;
+   maxCorner_[1] += eps;
+   maxCorner_[2] += eps;
 
    for( uint_t i = 0; i < 3; ++i )
       if( minCorner_[i] > maxCorner_[i] )
@@ -1544,14 +1544,14 @@ void GenericAABB< T >::extend( const value_type e )
 /**
  * \brief Extends this GenericAABB
  *
- * \param e epsilon vector by which the bounding box is extended. The box is extended in each direction by
+ * \param eps epsilon vector by which the bounding box is extended. The box is extended in each direction by
  *          the corresponding vector component.
  */
 template< typename T >
-void GenericAABB< T >::extend( const vector_type & e )
+void GenericAABB< T >::extend( const vector_type & eps )
 {
-   minCorner_ -= e;
-   maxCorner_ += e;
+   minCorner_ -= eps;
+   maxCorner_ += eps;
 
    for( uint_t i = 0; i < 3; ++i )
       if( minCorner_[i] > maxCorner_[i] )
diff --git a/src/core/math/equation_system/EquationParser.cpp b/src/core/math/equation_system/EquationParser.cpp
index 4a3894268b19b4f4b88cd4248caffd42b32d4456..aaf1b146d72c91f4bc429ba8d628e7010b897fa7 100644
--- a/src/core/math/equation_system/EquationParser.cpp
+++ b/src/core/math/equation_system/EquationParser.cpp
@@ -189,12 +189,12 @@ NodePtr EquationParser::parseFunction( const std::string& str, size_t& index ) c
    {
    case OP_FUNC_EXP:
       funcPtr = std::make_shared<Node>( OP_PROD );
-      funcPtr->left()  = std::make_shared<Node>( M_E  );
+      funcPtr->left()  = std::make_shared<Node>( math::e  );
       funcPtr->right() = nodePtr;
       return funcPtr;
    case OP_FUNC_LN:
       funcPtr = std::make_shared<Node>( OP_LOG );
-      funcPtr->right() = std::make_shared<Node>( M_E  );
+      funcPtr->right() = std::make_shared<Node>( math::e  );
       funcPtr->left()  = nodePtr;
       return funcPtr;
    case OP_FUNC_SQRT:
diff --git a/src/mesa_pd/data/shape/Sphere.h b/src/mesa_pd/data/shape/Sphere.h
index 50106bbbc7ecc5918b574a921152ce326be98df7..79137ea18b38369439263f616bd7d0f439c417f9 100644
--- a/src/mesa_pd/data/shape/Sphere.h
+++ b/src/mesa_pd/data/shape/Sphere.h
@@ -37,7 +37,7 @@ public:
 
    void updateMassAndInertia(const real_t density) override;
 
-   real_t getVolume() const override { return (real_t(4) / real_t(3)) * math::M_PI * getRadius() * getRadius() * getRadius(); }
+   real_t getVolume() const override { return (real_t(4) / real_t(3)) * math::pi * getRadius() * getRadius() * getRadius(); }
 
    static const int SHAPE_TYPE = 1; ///< Unique shape type identifier for spheres.\ingroup mesa_pd_shape
 
@@ -48,7 +48,7 @@ private:
 inline
 void Sphere::updateMassAndInertia(const real_t density)
 {
-   const real_t m = (real_c(4.0)/real_c(3.0) * math::M_PI) * getRadius() * getRadius() * getRadius() * density;
+   const real_t m = (real_c(4.0)/real_c(3.0) * math::pi) * getRadius() * getRadius() * getRadius() * density;
    const Mat3   I = Mat3::makeDiagonalMatrix( real_c(0.4) * m * getRadius() * getRadius() );
 
    mass_         = m;
diff --git a/src/mesa_pd/kernel/SpringDashpot.h b/src/mesa_pd/kernel/SpringDashpot.h
index fdc9fc6a7702298fa0e409199267580afaa384e9..60957e0fa374d883232a65531434c519183df168 100644
--- a/src/mesa_pd/kernel/SpringDashpot.h
+++ b/src/mesa_pd/kernel/SpringDashpot.h
@@ -100,7 +100,7 @@ public:
                                        const real_t meff)
    {
       auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
-      return std::exp(-a * math::M_PI / std::sqrt(getStiffness(type1, type2) / meff - a*a));
+      return std::exp(-a * math::pi / std::sqrt(getStiffness(type1, type2) / meff - a*a));
    }
 
    inline
@@ -109,7 +109,7 @@ public:
                             const real_t meff)
    {
       auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
-      return math::M_PI / std::sqrt( getStiffness(type1, type2)/meff - a*a);
+      return math::pi / std::sqrt( getStiffness(type1, type2)/meff - a*a);
    }
 
    inline
@@ -120,8 +120,8 @@ public:
                              const real_t meff)
    {
       const real_t lnDryResCoeff = std::log(cor);
-      setStiffness(type1, type2, math::M_PI * math::M_PI * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  ));
-      setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffness(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) ));
+      setStiffness(type1, type2, math::pi * math::pi * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  ));
+      setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffness(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) ));
    }
 private:
    uint_t numParticleTypes_;
diff --git a/src/pe/collision/EPA.cpp b/src/pe/collision/EPA.cpp
index 8e93e89a5521d27bf6f2e2afcec574ad692632b8..85775e5a89893aa2638ac7ae0cb58d64c750d7a6 100644
--- a/src/pe/collision/EPA.cpp
+++ b/src/pe/collision/EPA.cpp
@@ -536,7 +536,7 @@ inline void EPA::createInitialSimplex( size_t numPoints, GeomPrimitive &geom1, G
       }
 
       Vec3 direction1 = (d % axis).getNormalized();
-      Quat q(d, (real_t(2.0)/real_t(3.0)) * real_t(walberla::math::M_PI));
+      Quat q(d, (real_t(2.0)/real_t(3.0)) * real_t(walberla::math::pi));
       Mat3 rot = q.toRotationMatrix();
       Vec3 direction2 = (rot*direction1).getNormalized();
       Vec3 direction3 = (rot*direction2).getNormalized();
diff --git a/src/pe/cr/HCSITS.impl.h b/src/pe/cr/HCSITS.impl.h
index a03c3107fdc7e6162a21af9fd353ec19358e6d3e..ba8fd6aa29d0ec028448b838838406c9d49c50f2 100644
--- a/src/pe/cr/HCSITS.impl.h
+++ b/src/pe/cr/HCSITS.impl.h
@@ -1203,8 +1203,8 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticGenerali
                                                       alpha_left = -angleI - shiftI;
                                                       alpha_right = +angleI - shiftI;
                                                       if( alpha_left < 0 ) {
-                                                         alpha_left += 2 * math::M_PI;
-                                                         alpha_right += 2 * math::M_PI;
+                                                         alpha_left += 2 * math::pi;
+                                                         alpha_right += 2 * math::pi;
                                                       }
                                                    }
                                                    else if( contactCache.diag_nto_[i](0, 0) > contactCache.mu_[i] * a3 ) {
@@ -1213,8 +1213,8 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticGenerali
                alpha_left = -angleJ - shiftJ;
                alpha_right = +angleJ - shiftJ;
                if( alpha_left < 0 ) {
-                  alpha_left += 2 * math::M_PI;
-                  alpha_right += 2 * math::M_PI;
+                  alpha_left += 2 * math::pi;
+                  alpha_right += 2 * math::pi;
                }
             }
             else {
@@ -1223,15 +1223,15 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticGenerali
                real_t alpha1_left( -angleJ - shiftJ );
                real_t alpha1_right( +angleJ - shiftJ );
                if( alpha1_left < 0 ) {
-                  alpha1_left += 2 * math::M_PI;
-                  alpha1_right += 2 * math::M_PI;
+                  alpha1_left += 2 * math::pi;
+                  alpha1_right += 2 * math::pi;
                }
                const real_t angleI( std::acos( fractionI ) );
                real_t alpha2_left( -angleI - shiftI );
                real_t alpha2_right( +angleI - shiftI );
                if( alpha2_left < 0 ) {
-                  alpha2_left += 2 * math::M_PI;
-                  alpha2_right += 2 * math::M_PI;
+                  alpha2_left += 2 * math::pi;
+                  alpha2_right += 2 * math::pi;
                }
 
                // Swap intervals if second interval does not start right of the first interval.
@@ -1241,7 +1241,7 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticGenerali
                }
 
                if( alpha2_left > alpha1_right ) {
-                  alpha2_right -= 2*math::M_PI;
+                  alpha2_right -= 2*math::pi;
                   if( alpha2_right > alpha1_right ) {
                      // [alpha1_left; alpha1_right] \subset [alpha2_left; alpha2_right]
                   }
@@ -1807,7 +1807,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::integratePositions( Body
             v = v * (edge * getSpeedLimitFactor() / dt / speed );
          }
 
-         const real_t maxPhi = real_t(2) * math::M_PI * getSpeedLimitFactor();
+         const real_t maxPhi = real_t(2) * math::pi * getSpeedLimitFactor();
          const real_t phi    = w.length() * dt;
          if (phi > maxPhi)
          {
diff --git a/src/pe/raytracing/Raytracer.cpp b/src/pe/raytracing/Raytracer.cpp
index 6c23292f9ef86459a836ca2de7833f44bead1361..833089fcd874c135859b2f30119b03c7960d166d 100644
--- a/src/pe/raytracing/Raytracer.cpp
+++ b/src/pe/raytracing/Raytracer.cpp
@@ -197,7 +197,7 @@ void Raytracer::setupView_() {
    // viewing plane setup
    d_ = (cameraPosition_ - lookAtPoint_).length();
    aspectRatio_ = real_t(pixelsHorizontal_) / real_t(pixelsVertical_);
-   real_t fov_vertical_rad = fov_vertical_ * math::M_PI / real_t(180.0);
+   real_t fov_vertical_rad = fov_vertical_ * math::pi / real_t(180.0);
    viewingPlaneHeight_ = real_c(tan(fov_vertical_rad/real_t(2.))) * real_t(2.) * d_;
    viewingPlaneWidth_ = viewingPlaneHeight_ * aspectRatio_;
    viewingPlaneOrigin_ = lookAtPoint_ - u_*viewingPlaneWidth_/real_t(2.) - v_*viewingPlaneHeight_/real_t(2.);
diff --git a/src/pe/rigidbody/Capsule.cpp b/src/pe/rigidbody/Capsule.cpp
index f1d9ad9ee146d699bb7f3b062303454e394ed815..b544f08518b9dfe6ca6f66064e63dcc20b553445 100644
--- a/src/pe/rigidbody/Capsule.cpp
+++ b/src/pe/rigidbody/Capsule.cpp
@@ -194,8 +194,8 @@ void Capsule::calcBoundingBox()
  */
 Mat3 Capsule::calcInertia( const real_t radius, const real_t length, const real_t density)
 {
-   const real_t  sphereMass( real_t (4)/real_t (3) * math::M_PI * radius*radius*radius * density );
-   const real_t  cylinderMass( math::M_PI * radius*radius * length * density );
+   const real_t  sphereMass( real_t (4)/real_t (3) * math::pi * radius*radius*radius * density );
+   const real_t  cylinderMass( math::pi * radius*radius * length * density );
 
    // 'Ia' represent the moment of inertia along the x-axis. 'Ia' contains the following two parts:
    //  - cylinder :  I = (1/2)*mass*radius^2
diff --git a/src/pe/rigidbody/Capsule.h b/src/pe/rigidbody/Capsule.h
index 8718f4c90768026a8a90999c8ed569635eb333ed..937aadcda652d09631d4ecb0702eda7a27acafcb 100644
--- a/src/pe/rigidbody/Capsule.h
+++ b/src/pe/rigidbody/Capsule.h
@@ -223,7 +223,7 @@ inline real_t  Capsule::getVolume() const
  */
 inline real_t  Capsule::calcVolume( real_t  radius, real_t  length )
 {
-   return math::M_PI*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length );
+   return math::pi*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length );
 }
 //*************************************************************************************************
 
@@ -238,7 +238,7 @@ inline real_t  Capsule::calcVolume( real_t  radius, real_t  length )
  */
 inline real_t  Capsule::calcMass( real_t  radius, real_t  length, real_t  density )
 {
-   return math::M_PI*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length ) * density;
+   return math::pi*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length ) * density;
 }
 //*************************************************************************************************
 
@@ -253,7 +253,7 @@ inline real_t  Capsule::calcMass( real_t  radius, real_t  length, real_t  densit
  */
 inline real_t  Capsule::calcDensity( real_t  radius, real_t  length, real_t  mass )
 {
-   return mass / ( math::M_PI*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length ) );
+   return mass / ( math::pi*radius*radius*( ( static_cast<real_t >( 4 ) / static_cast<real_t >( 3 ) )*radius + length ) );
 }
 //*************************************************************************************************
 
diff --git a/src/pe/rigidbody/Ellipsoid.h b/src/pe/rigidbody/Ellipsoid.h
index 1cda4738157ba98c7aabe5439a4ff9547ceddbe3..958c78ab31f474f99c5bdc310023fa0ca79edc31 100644
--- a/src/pe/rigidbody/Ellipsoid.h
+++ b/src/pe/rigidbody/Ellipsoid.h
@@ -203,7 +203,7 @@ inline real_t Ellipsoid::getVolume() const
  */
 inline real_t Ellipsoid::calcVolume(const Vec3& semiAxes ) 
 {
-   return real_c(4.0)/real_c(3.0) * math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2];
+   return real_c(4.0)/real_c(3.0) * math::pi * semiAxes[0] * semiAxes[1] * semiAxes[2];
 }
 //*************************************************************************************************
 
@@ -217,7 +217,7 @@ inline real_t Ellipsoid::calcVolume(const Vec3& semiAxes )
  */
 inline real_t Ellipsoid::calcMass(const Vec3& semiAxes, real_t density )
 {
-   return real_c(4.0)/real_c(3.0) * math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2] * density;
+   return real_c(4.0)/real_c(3.0) * math::pi * semiAxes[0] * semiAxes[1] * semiAxes[2] * density;
 }
 //*************************************************************************************************
 
@@ -231,7 +231,7 @@ inline real_t Ellipsoid::calcMass(const Vec3& semiAxes, real_t density )
  */
 inline real_t Ellipsoid::calcDensity(const Vec3& semiAxes, real_t mass )
 {
-   return real_c(0.75) * mass / ( math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2] );
+   return real_c(0.75) * mass / ( math::pi * semiAxes[0] * semiAxes[1] * semiAxes[2] );
 }
 //*************************************************************************************************
 
diff --git a/src/pe/rigidbody/Sphere.h b/src/pe/rigidbody/Sphere.h
index c4a1c95fc710dbf10010730df7cb7ef6dbbaab53..669cf07bcb15805f80c45f7992a4c31c7233678d 100644
--- a/src/pe/rigidbody/Sphere.h
+++ b/src/pe/rigidbody/Sphere.h
@@ -211,7 +211,7 @@ inline real_t Sphere::getVolume() const
  */
 inline real_t Sphere::calcVolume( real_t radius )
 {
-   return real_c(4.0)/real_c(3.0) * math::M_PI * radius * radius * radius;
+   return real_c(4.0)/real_c(3.0) * math::pi * radius * radius * radius;
 }
 //*************************************************************************************************
 
@@ -225,7 +225,7 @@ inline real_t Sphere::calcVolume( real_t radius )
  */
 inline real_t Sphere::calcMass( real_t radius, real_t density )
 {
-   return real_c(4.0)/real_c(3.0) * math::M_PI * radius * radius * radius * density;
+   return real_c(4.0)/real_c(3.0) * math::pi * radius * radius * radius * density;
 }
 //*************************************************************************************************
 
@@ -239,7 +239,7 @@ inline real_t Sphere::calcMass( real_t radius, real_t density )
  */
 inline real_t Sphere::calcDensity( real_t radius, real_t mass )
 {
-   return real_c(0.75) * mass / ( math::M_PI * radius * radius * radius );
+   return real_c(0.75) * mass / ( math::pi * radius * radius * radius );
 }
 //*************************************************************************************************
 
diff --git a/src/pe/utility/Overlap.cpp b/src/pe/utility/Overlap.cpp
index 2ed2fe1c943f3349b71d86e1917baed710194e77..6c7658dce7cd4c30d99ce67d02be38932b752e50 100644
--- a/src/pe/utility/Overlap.cpp
+++ b/src/pe/utility/Overlap.cpp
@@ -34,7 +34,7 @@ real_t getSphereSphereOverlap(const real_t d, const real_t r1, const real_t r2)
    //http://math.stackexchange.com/questions/297751/overlapping-spheres
    if (d > r1+r2)
       return 0; else
-      return math::M_PI / (real_c(12.0)  * d) * (r1 + r2 - d) * (r1 + r2 - d) * (d*d + 2*d*(r1+r2) - 3*(r1-r2)*(r1-r2));
+      return math::pi / (real_c(12.0)  * d) * (r1 + r2 - d) * (r1 + r2 - d) * (d*d + 2*d*(r1+r2) - 3*(r1-r2)*(r1-r2));
 }
 
 }  // namespace pe
diff --git a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
index 8d7fc55c0fa014f6eddd1189618b5a985de1dc67..1db90c9457a398846ba1f478ee33ebc3a5da2333 100644
--- a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
+++ b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
@@ -54,10 +54,10 @@ real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
 }
 
 // Coefficient from Stokes' law for drag, only valid for Stokes regime (low Reynolds numbers)
-// = 3 * M_PI * mu * D * fluidVolumeFraction
+// = 3 * math::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;
+   return real_t(3) * math::pi * diameter * fluidDynamicViscosity * fluidVolumeFraction;
 }
 
 // threshold value for absolute relative velocity
@@ -172,7 +172,7 @@ Vector3<real_t> dragForceFelice( const Vector3<real_t> & fluidVel, const Vector3
    real_t temp2 = real_t(1.5) - std::log10( reynoldsNumber );
    real_t chi = real_t(3.7) - std::pow( real_t(0.65), (- real_t(0.5) * temp2 * temp2 ) );
 
-   return real_t(0.125) * dragCoeff * fluidDensity * math::M_PI * diameter * diameter * absVelDiff *
+   return real_t(0.125) * dragCoeff * fluidDensity * math::pi * diameter * diameter * absVelDiff *
           std::pow( fluidVolumeFraction, real_t(2) - chi) * velDiff;
 
 }
@@ -206,7 +206,7 @@ Vector3<real_t> dragForceTenneti( const Vector3<real_t> & fluidVel, const Vector
 
    const real_t CdRe = fluidVolumeFraction * ( CdRe0Sphere / fvfCubed + A + B );
 
-   return real_t(3) * math::M_PI * diameter * fluidDynamicViscosity * fluidVolumeFraction * CdRe * velDiff;
+   return real_t(3) * math::pi * diameter * fluidDynamicViscosity * fluidVolumeFraction * CdRe * velDiff;
 
 }
 
diff --git a/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h b/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
index cd036bfc064357692fa645c7ad99df529f5910f3..2de86148379f3ed7391a09477be7712b44f264e7 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::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 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(9)/real_t(84) ) * ( h / d ) * std::log( d/( real_t(2)*h ) ) );
-   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);
+   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);
    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::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 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(9)/real_t(84) ) * ( h / d ) * std::log( d/( real_t(2)*h ) ) );
-   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);
+   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);
    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/geometry/SphereEquivalentDiameter.h b/src/pe_coupling/geometry/SphereEquivalentDiameter.h
index d3750379167bc8506ec25c6e2d103a19c811b986..7e40a843dac7f0d5aca20cb996a5c28eb2914e08 100644
--- a/src/pe_coupling/geometry/SphereEquivalentDiameter.h
+++ b/src/pe_coupling/geometry/SphereEquivalentDiameter.h
@@ -39,7 +39,7 @@ real_t getSphereEquivalentDiameter( pe::RigidBody & body )
       real_t radius = sphere.getRadius();
       return real_t(2) * radius;
    } else {
-      const real_t preFac = real_t(6) / math::M_PI;
+      const real_t preFac = real_t(6) / math::pi;
       return std::cbrt( body.getVolume() * preFac );
    }
 }
diff --git a/src/pe_coupling/utility/LubricationCorrection.cpp b/src/pe_coupling/utility/LubricationCorrection.cpp
index 6cee7eb439df47af3c8051e9a8fb5d2c882e175d..a20a0ca4be3c0951278695e6588ef1c2fa051360 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::M_PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+   pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::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::M_PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+   pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::pi * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
 
    WALBERLA_LOG_DETAIL_SECTION() {
       std::stringstream ss;
diff --git a/tests/core/math/PlaneTest.cpp b/tests/core/math/PlaneTest.cpp
index 2f895119184e14b303536a3a84a44cdd707a8d68..6c0a49a95886e62d830f782f2c91ee0fbb49c2a8 100644
--- a/tests/core/math/PlaneTest.cpp
+++ b/tests/core/math/PlaneTest.cpp
@@ -138,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::M_PI / real_t(180) )
+      if( (p0 - p1).sqrLength() < 1e-6 || (p0 - p2).sqrLength() < 1e-6 || (p2 - p1).sqrLength() < 1e-6 || angle < math::pi / real_t(180) )
       {
          --i;
          continue;
diff --git a/tests/core/timing/TimerTest.cpp b/tests/core/timing/TimerTest.cpp
index a6ef83f9ee0fe8680089813d429c575ebf997acf..b9ee91ec15be2caa8e2c3637702a4f058747c7ba 100644
--- a/tests/core/timing/TimerTest.cpp
+++ b/tests/core/timing/TimerTest.cpp
@@ -31,7 +31,7 @@ namespace walberla {
 static double burnTime()
 {
    double sum = 0.0;
-   for( double d = 0.0; d < math::M_PI; d += 0.000001 )
+   for( double d = 0.0; d < math::pi; d += 0.000001 )
    {
       sum += std::atan( std::tan( d ) );
       sum += std::asin( std::sin( d ) );
diff --git a/tests/core/timing/TimingPoolTest.cpp b/tests/core/timing/TimingPoolTest.cpp
index e3ff6f3912453cabde8a975b74ca4c519bb16197..db13dc3bd7e2a6a3fa140de4dada2f61f8ac07b4 100644
--- a/tests/core/timing/TimingPoolTest.cpp
+++ b/tests/core/timing/TimingPoolTest.cpp
@@ -62,7 +62,7 @@ void scopedTimer()
       auto scopedTimer = pool.getScopeTimer( "scope timer" );
 
       double sum = 0.0;
-      for( double d = 0.0; d < math::M_PI; d += 0.00001 )
+      for( double d = 0.0; d < math::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 aaa669c2ac0f79048e379e09919afdba22412ec6..cb8797f187c217f1ac787288a4737ee6a8f84456 100644
--- a/tests/fft/GreensTest.cpp
+++ b/tests/fft/GreensTest.cpp
@@ -56,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(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])) -
+      return real_c(0.5) / ( std::cos( real_c(2) * real_c(math::pi) * real_c(x) / real_c(dim[0])) +
+                             std::cos( real_c(2) * real_c(math::pi) * real_c(y) / real_c(dim[1])) +
+                             std::cos( real_c(2) * real_c(math::pi) * real_c(z) / real_c(dim[2])) -
                              real_c(3) ) / real_c(dim[0]*dim[1]*dim[2]);
    };
    
diff --git a/tests/geometry/ScalarFieldFromBodyTest.cpp b/tests/geometry/ScalarFieldFromBodyTest.cpp
index 0aa382715a8d664570635d958a86fd3f23bc58cc..47d8d8274d32db43343c7984646643e7aef0fe44 100644
--- a/tests/geometry/ScalarFieldFromBodyTest.cpp
+++ b/tests/geometry/ScalarFieldFromBodyTest.cpp
@@ -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] * math::M_PI;
+   const real_t expectedVolume = real_t(4) / real_t(3) * radii[0] * radii[1] * radii[2] * math::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 * math::M_PI;
+   const real_t expectedVolume = real_t(4) / real_t(3) * radius * radius * radius * math::pi;
 
    resetField( storage, fieldID );
 
diff --git a/tests/lbm/DiffusionTest.cpp b/tests/lbm/DiffusionTest.cpp
index cdb6c71a057d209f0c3f4a9c38797ec4a486c34d..a01187c455108ab7d76597685d2e835183df4b67 100644
--- a/tests/lbm/DiffusionTest.cpp
+++ b/tests/lbm/DiffusionTest.cpp
@@ -236,8 +236,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::M_PI / real_c( timesteps  );
-   const real_t cperiod = real_t(2) * math::M_PI / real_c( cells[dim] );
+   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] );
 
    // --- create blockstorage --- //
    auto blockStorage = blockforest::createUniformBlockGrid(
diff --git a/tests/lbm/boundary/BoundaryForce.cpp b/tests/lbm/boundary/BoundaryForce.cpp
index cc1f9e63acf10569164968c64b55d24eb683c03c..2182030d08dfc5d370949e680a4105cf782e9efe 100644
--- a/tests/lbm/boundary/BoundaryForce.cpp
+++ b/tests/lbm/boundary/BoundaryForce.cpp
@@ -211,7 +211,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::M_PI * visc * R * velocity;
+   Vector3<real_t> stokes = 6 * math::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/DiffusionDirichlet.cpp b/tests/lbm/boundary/DiffusionDirichlet.cpp
index ec0f1c51071aaa14108a97313b833df34e027e18..b9215d4321f5445f0c04360ede0c8b8dd674aaa3 100644
--- a/tests/lbm/boundary/DiffusionDirichlet.cpp
+++ b/tests/lbm/boundary/DiffusionDirichlet.cpp
@@ -109,7 +109,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::M_PI/L ),
+      period_( real_t(2)*math::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) ),
@@ -254,7 +254,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::M_PI/real_c(length) );
+   auto cbc = make_shared<CosBoundaryConfiguration>( real_t(2)*math::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 699d2fe1bff055b51fa00c2f0adbdd499976bbff..ab7efa8e037dc4adb0a60c9bd57731a6c53cb799 100644
--- a/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp
+++ b/tests/lbm/boundary/SimpleDiffusionDirichlet.cpp
@@ -237,8 +237,8 @@ public:
          delta_( maxValue - minValue),
          length_(real_c(length)),
          lengthInv_(real_t(1)/real_c(length)),
-         pi_(math::M_PI),
-         piInv_(real_t(1)/math::M_PI),
+         pi_(math::pi),
+         piInv_(real_t(1)/math::pi),
          valid_(uint_c(std::ceil(omega*omega*omega*real_t(10)))),
          time_( time ),
          expArray(),
diff --git a/tests/lbm/evaluations/PermeabilityTest.cpp b/tests/lbm/evaluations/PermeabilityTest.cpp
index 7826539357a236d6cc3976631e9cf6163fcd65fc..cb224e30fe8f837110de353ca10a0051c56881d1 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::M_PI * r * real_t(2) * drag );
+   return ( L * L * L ) / ( real_t(6) * math::pi * r * real_t(2) * drag );
 }
 
 
diff --git a/tests/mesa_pd/common/IntersectionRatio.cpp b/tests/mesa_pd/common/IntersectionRatio.cpp
index 0451362c529aee088ef83b3e8158b8196fa0aac5..324f91558382411bef521c53f652df5f427d43d5 100644
--- a/tests/mesa_pd/common/IntersectionRatio.cpp
+++ b/tests/mesa_pd/common/IntersectionRatio.cpp
@@ -140,7 +140,7 @@ int main( int argc, char **argv )
       auto planeShape = shapeStorage->create<mesa_pd::data::HalfSpace>( normal.getNormalized() );
 
       // rotate to same position as half space before
-      Vector3<real_t> rotationAngles( -math::M_PI / real_t(4), real_t(0), real_t(0));
+      Vector3<real_t> rotationAngles( -math::pi / real_t(4), real_t(0), real_t(0));
       Quaternion<real_t> quat( rotationAngles );
 
       mesa_pd::data::Particle&& p = *ps->create(true);
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
index 53cff74ea690c4e55abb18894f80145063cae2ff..894cc47ee0722e9d836ed71cf4013329108ce57e 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
@@ -108,8 +108,8 @@ int main( int argc, char** argv )
    const real_t particleMass = real_t(1) / ss->shapes[sphereShape]->getInvMass();
    const real_t Mij = particleMass; // Mij = M for sphere-wall collision
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
-   const real_t stiffnessN = math::M_PI * math::M_PI * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  );
-   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) );
+   const real_t stiffnessN = math::pi * math::pi * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  );
+   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) );
 
    WALBERLA_LOG_INFO("dt = " << dt << ", Tc = " << collisionTime << ", coefficient of restitution = " << restitutionCoeff);
    WALBERLA_LOG_INFO(" -> mass " << particleMass << ", collision duration = " << collisionTime / dt << ", stiffness = " << stiffnessN << ", damping = " << dampingN);
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
index 5b25815de2fd7714b3eb96032cf16eeabca28b0d..0fa4a40cd8ec2edc9a55cbd59f14ab2544f473b5 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
@@ -109,8 +109,8 @@ int main( int argc, char** argv )
    const real_t particleMass = real_t(1) / ss->shapes[sphereShape]->getInvMass();
    const real_t Mij = particleMass; // Mij = M for sphere-wall collision
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
-   const real_t stiffnessN = math::M_PI * math::M_PI * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  );
-   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) );
+   const real_t stiffnessN = math::pi * math::pi * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  );
+   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) );
 
    WALBERLA_LOG_INFO("dt = " << dt << ", Tc = " << collisionTime << ", coefficient of restitution = " << restitutionCoeff);
    WALBERLA_LOG_INFO(" -> mass " << particleMass << ", collision duration = " << collisionTime / dt << ", stiffness = " << stiffnessN << ", damping = " << dampingN);
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
index 6c45e39cff3d462206fc51100ad0a304e8ba0124..ab3e92e67c71ba7c35b25744fae1a6522296833a 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
@@ -108,8 +108,8 @@ int main( int argc, char** argv )
    const real_t particleMass = real_t(1) / ss->shapes[sphereShape]->getInvMass();
    const real_t Mij = particleMass; // Mij = M for sphere-wall collision
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
-   const real_t stiffnessN = math::M_PI * math::M_PI * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  );
-   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) );
+   const real_t stiffnessN = math::pi * math::pi * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  );
+   const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) );
 
    WALBERLA_LOG_INFO("dt = " << dt << ", Tc = " << collisionTime << ", coefficient of restitution = " << restitutionCoeff);
    WALBERLA_LOG_INFO(" -> mass " << particleMass << ", collision duration = " << collisionTime / dt << ", stiffness = " << stiffnessN << ", damping = " << dampingN);
diff --git a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
index 48f62346600f57a3404f880b9b6684bc1de3be75..0992675235cb5cc4eba1f9eed90f4140130ede8b 100644
--- a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
+++ b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
@@ -90,12 +90,12 @@ int main( int argc, char ** argv )
    }
 
 
-   real_t phase = phaseFraction * math::M_PI;
+   real_t phase = phaseFraction * math::pi;
    real_t omega = std::sqrt(k / mass);
-   real_t durationOnePeriod = real_t(2) * math::M_PI / omega;
+   real_t durationOnePeriod = real_t(2) * math::pi / omega;
    uint_t timeSteps = uint_c(periods * durationOnePeriod / dt);
 
-   WALBERLA_LOG_INFO("omega = " << omega << ", T = " << real_t(2) * math::M_PI / omega << ", time steps = " << timeSteps << ", phase = " << phase << ", periods = " << periods);
+   WALBERLA_LOG_INFO("omega = " << omega << ", T = " << real_t(2) * math::pi / omega << ", time steps = " << timeSteps << ", phase = " << phase << ", periods = " << periods);
 
    //initialize particle
    const auto pos = Vec3(0,0,analyticalTrajectory(amplitude, real_t(0), omega, phase));
diff --git a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
index 6ff5bede797db82bd9bf02aa970f082c5b7391ce..0513e8086b36d10519291470960b59d2f4036eed 100644
--- a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
@@ -117,16 +117,16 @@ int main( int argc, char ** argv )
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
 
    // normal material aprameters
-   const real_t stiffnessN = math::M_PI * math::M_PI * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ))  );
+   const real_t stiffnessN = math::pi * math::pi * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ))  );
    const real_t dampingN = - real_t(2) * std::sqrt( Mij * stiffnessN ) *
-   ( lnDryResCoeff / std::sqrt( math::M_PI * math::M_PI + ( lnDryResCoeff * lnDryResCoeff ) ) );
+   ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) );
 
    WALBERLA_LOG_INFO_ON_ROOT("normal: stiffness = " << stiffnessN << ", damping = " << dampingN);
 
    const real_t lnDryResCoeffTangential = lnDryResCoeff; // std::log(0.31); //TODO: was same as in normal direction
    const real_t kappa = real_t(2) * ( real_t(1) - nu ) / ( real_t(2) - nu ) ;
-   const real_t stiffnessT = kappa * Mij * math::M_PI * math::M_PI / ( collisionTime *  collisionTime );
-   const real_t dampingT = real_t(2) * std::sqrt(Mij * stiffnessT) * ( - lnDryResCoeffTangential ) / ( std::sqrt( math::M_PI * math::M_PI + lnDryResCoeffTangential * lnDryResCoeffTangential ));
+   const real_t stiffnessT = kappa * Mij * math::pi * math::pi / ( collisionTime *  collisionTime );
+   const real_t dampingT = real_t(2) * std::sqrt(Mij * stiffnessT) * ( - lnDryResCoeffTangential ) / ( std::sqrt( math::pi * math::pi + lnDryResCoeffTangential * lnDryResCoeffTangential ));
 
    WALBERLA_LOG_INFO_ON_ROOT("tangential: kappa = " << kappa << ", stiffness T = " << stiffnessT << ", damping T = " << dampingT);
 
diff --git a/tests/mesh/MeshPeRaytracing.cpp b/tests/mesh/MeshPeRaytracing.cpp
index a31dd985a9747eec26394629fe0ed58540ecd584..038c78a34cbfb0f9b6d405851ecc020c18dc8563 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) * real_t(M_PI) * rand1;
+         real_t theta = real_t(2) * real_t(math::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));
 
@@ -94,7 +94,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) * math::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 04d44473b570d53fed525cbfa57783ecb79aeeff..589f46642b6ca902f52cc4e379b2f32200dbc245 100644
--- a/tests/mesh/PeVTKMeshWriterTest.cpp
+++ b/tests/mesh/PeVTKMeshWriterTest.cpp
@@ -96,7 +96,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 * real_t(M_PI) * distribution(rng);
+      real_t theta = 2 * real_t(math::pi) * distribution(rng);
       real_t phi = std::acos( real_t(1.0) - real_t(2.0) * distribution(rng) );
       p[0] = std::sin(phi) * std::cos(theta) * radius;
       p[1] = std::sin(phi) * std::sin(theta) * radius;
diff --git a/tests/mesh/QHullTest.cpp b/tests/mesh/QHullTest.cpp
index 78bd8d12001c9d8a3bbe2b133366c371cd99d917..5e17b97e844d3b0523c52c410e6649134ac6b7c8 100644
--- a/tests/mesh/QHullTest.cpp
+++ b/tests/mesh/QHullTest.cpp
@@ -23,6 +23,7 @@
 #include "core/logging/Logging.h"
 #include "core/mpi/Environment.h"
 #include "core/timing/Timer.h"
+#include "core/math/Constants.h"
 #include "core/math/Utility.h"
 
 #include "mesh/TriangleMeshes.h"
@@ -209,7 +210,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 * real_t(M_PI) * distribution(rng);
+      real_t theta = 2 * real_t(math::pi) * distribution(rng);
       real_t phi = std::acos( real_t(1.0) - real_t(2.0) * distribution(rng) );
       p[0] = std::sin(phi) * std::cos(theta) * radius;
       p[1] = std::sin(phi) * std::sin(theta) * radius;
diff --git a/tests/pde/BoundaryTest.cpp b/tests/pde/BoundaryTest.cpp
index c232591ada427488f37cf8679e5949703396f5f1..df548787a45c42fe7ade1d91a878601c7772361b 100644
--- a/tests/pde/BoundaryTest.cpp
+++ b/tests/pde/BoundaryTest.cpp
@@ -132,7 +132,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::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] );
+         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] );
       }
    }
 }
@@ -178,7 +178,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::M_PI * p[0] ) * std::sinh( real_t( 2 ) * math::M_PI * p[1] );
+         real_t val = std::sin( real_t( 2 ) * math::pi * p[0] ) * std::sinh( real_t( 2 ) * math::pi * p[1] );
 
          boundaryHandling->forceBoundary( Dirichlet_Flag, cell->x(), cell->y(), cell->z(), pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( val ) );
 
@@ -308,7 +308,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::M_PI * math::M_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::pi * math::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 81dd0e18cb41005b5430634fb55aea5b84cfc8af..a14832467fa0615e4720c330372cec29aea9b7c2 100644
--- a/tests/pde/CGTest.cpp
+++ b/tests/pde/CGTest.cpp
@@ -66,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::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::pi * p[0] ) * std::sinh( real_t(2) * math::pi * p[1] );
          }
       }
    }
@@ -83,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::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] );
+         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] );
       }
    }
 }
@@ -165,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::M_PI * math::M_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::pi * math::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/JacobiTest.cpp b/tests/pde/JacobiTest.cpp
index 525e0e063b8e5b8a47a08336c4ba6ecae83756d8..2e65dc6ef3c69c9c19469fff47e15f2b91ffad2d 100644
--- a/tests/pde/JacobiTest.cpp
+++ b/tests/pde/JacobiTest.cpp
@@ -70,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::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] );
+            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] );
          }
       }
    }
@@ -88,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::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] );
+         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] );
       }
    }
 }
@@ -168,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::M_PI * math::M_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::pi * math::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/RBGSTest.cpp b/tests/pde/RBGSTest.cpp
index 8609ace0af68ecbfeca6b75c22ca8aed39e3ded8..6417395dbd1d10a12bb9e243212e1b88c4a8911f 100644
--- a/tests/pde/RBGSTest.cpp
+++ b/tests/pde/RBGSTest.cpp
@@ -69,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::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::pi * p[0] ) * std::sinh( real_t(2) * math::pi * p[1] );
          }
       }
    }
@@ -86,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::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] );
+         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] );
       }
    }
 }
@@ -165,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::M_PI * math::M_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::pi * math::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/SORTest.cpp b/tests/pde/SORTest.cpp
index 0a5799163292f662469395294bde526e82706dde..374aa0b7c332e07da760a686681d446783082719 100644
--- a/tests/pde/SORTest.cpp
+++ b/tests/pde/SORTest.cpp
@@ -69,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::M_PI * p[0] ) * std::sinh( real_t(2) * math::M_PI * p[1] );
+            u->get( *cell ) = std::sin( real_t(2) * math::pi * p[0] ) * std::sinh( real_t(2) * math::pi * p[1] );
          }
       }
    }
@@ -86,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::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] );
+         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] );
       }
    }
 }
@@ -167,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::M_PI * math::M_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::pi * math::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/pe/Collision.cpp b/tests/pe/Collision.cpp
index a69762d7fb928ef9e96b2b93e2712ab259c9e0c2..6cd6f2377f89860a1600766406bbb979f40c4f62 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -96,7 +96,7 @@ void SphereTest()
    WALBERLA_CHECK(  collideFunc(&sp4, &cb1) );
    checkContact( contacts.at(0),
                  Contact( &sp4, &cb1, Vec3(0,real_t(2),real_t(0)), Vec3(0, -1, 0).getNormalized(), real_t(-0.5)) );
-   cb1.rotateAroundOrigin( Vec3( 0,0,1), math::M_PI * real_t(0.25) );
+   cb1.rotateAroundOrigin( Vec3( 0,0,1), math::pi * real_t(0.25) );
    WALBERLA_CHECK(  !collideFunc(&sp1, &cb1) );
    WALBERLA_CHECK(  collideFunc(&sp2, &cb1) );
    WALBERLA_CHECK(  collideFunc(&sp4, &cb1) );
@@ -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), Quat(), Vec3(2,2,2), iron, false, true, false);
-   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) );
+   b5.rotate( Vec3(0,0,1), real_t(math::pi * 0.25) );
+   b5.rotate( Vec3(1,0,0), real_t(math::pi * 0.25) );
 
    std::vector<Contact> contacts;
    fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
@@ -211,7 +211,7 @@ void CapsuleTest2()
    MaterialID     material = createMaterial( "granular", real_t( 1.0 ), 0, static_cof, dynamic_cof, real_t( 0.5 ), 1, 1, 0, 0 );
    //create obstacle
    Capsule c1(100, 100, Vec3(10,10,0), Quat(), 3, 40, material, false, true, false);
-   c1.rotate( Vec3(0,1,0), math::M_PI * real_t(0.5) );
+   c1.rotate( Vec3(0,1,0), math::pi * real_t(0.5) );
    Sphere sp1(123, 123, Vec3(real_t(6.5316496854295262864), real_t(10.099999999999999645), real_t(0.46999999991564372914) ), Quat(), real_t(0.47), material, false, true, false);
 
    std::vector<Contact> contacts;
diff --git a/tests/pe/CollisionTobiasGJK.cpp b/tests/pe/CollisionTobiasGJK.cpp
index 41bdf1f691d44ab36dfbe8b4435cdfdbc6531909..a2de970709a85388bef0015e45d3be666f08d844 100644
--- a/tests/pe/CollisionTobiasGJK.cpp
+++ b/tests/pe/CollisionTobiasGJK.cpp
@@ -207,7 +207,7 @@ void MainTest()
    //Testcase 04 Cube with turned Cube
    WALBERLA_LOG_INFO("Test 04: CUBE <-> TURNED CUBE");
    //compute rotation.
-   real_t angle = walberla::math::M_PI/real_t(4.0);
+   real_t angle = walberla::math::pi/real_t(4.0);
    Vec3 zaxis(0, 0, 1);
    Quat q4(zaxis, angle);
 
@@ -263,7 +263,7 @@ void MainTest()
    //Testcase 08:
    // CAPSULE <-> CAPSULE
    WALBERLA_LOG_INFO("Test 08: CAPSULE <-> CAPSULE");
-   Quat q8(Vec3(0,1,0), walberla::math::M_PI/real_t(2.0)); //creates a y-axis aligned capsule
+   Quat q8(Vec3(0,1,0), walberla::math::pi/real_t(2.0)); //creates a y-axis aligned capsule
    Capsule cap8_1(139, 18, Vec3(0,0,0), Quat(), real_t(4.0), real_t(10.0), iron, false, true, false);
    Capsule cap8_2(140, 19, Vec3(0,0, real_t(13.0)),  q8, real_t(4.0), real_t(10.0), iron, false, true, false);
    Vec3 wpm8(0, 0, real_t(-0.5));
diff --git a/tests/pe/GJK_EPA.cpp b/tests/pe/GJK_EPA.cpp
index 9293c3c866c440f89d06dde17fae3b8dcb14a525..918cf594221a37235fba4c01f8b6763a2681df45 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::M_PI * 180 << "°)" );
+//      WALBERLA_LOG_INFO("deltaNormal: |" << normal - normalGJK << "| = " << (normal - normalGJK).length() << " (" << acos(normal*normalGJK) / math::pi * 180 << "°)" );
 //      WALBERLA_LOG_INFO("deltaPen   : " << penetrationDepth << " - " << penetrationDepthGJK << " = " << penetrationDepth - penetrationDepthGJK);
 
       point.insert( (contactPoint - contactPointGJK).length() );
-      norm.insert( acos(normal*normalGJK) / math::M_PI * 180 );
+      norm.insert( acos(normal*normalGJK) / math::pi * 180 );
       pen.insert( fabs(penetrationDepth - penetrationDepthGJK) );
    }
 
diff --git a/tests/pe/PeDocumentationSnippets.cpp b/tests/pe/PeDocumentationSnippets.cpp
index eb667d526cefedef744e63ecba9cf8f533b55a60..6c4e49c51ef1a900efe6b0f5c8bcc7631ee519f4 100644
--- a/tests/pe/PeDocumentationSnippets.cpp
+++ b/tests/pe/PeDocumentationSnippets.cpp
@@ -97,14 +97,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::M_PI/3.0), 0.0 );
+      box->rotate( 0.0, real_c(math::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::M_PI/3.0), 0.0 );
+      capsule->rotate( 0.0, real_c(math::pi/3.0), 0.0 );
    //! [Create a Capsule]
 
    //! [Create a Plane]
@@ -117,7 +117,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::M_PI/3.0), 0.0 );
+      sphere->rotate( 0.0, real_c(math::pi/3.0), 0.0 );
    //! [Create a Sphere]
 
    //! [Create a Union]
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
index 15dd6ae60f0e7c8ea8e2d70310672c107cca90c4..f7aa2ca20cb45da786005f489e848f16f2bd6166 100644
--- a/tests/pe/Raytracing.cpp
+++ b/tests/pe/Raytracing.cpp
@@ -168,7 +168,7 @@ void BoxIntersectsTest() {
    WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(9.7068), real_t(1e-4));
    
    Box box5(128, 5, Vec3(4, 0, 0), Quat(), Vec3(4, 4, 4), iron, false, true, false);
-   box5.rotate(0,0,math::M_PI/4);
+   box5.rotate(0,0,math::pi/4);
    Ray ray5(Vec3(0,1.5,0), Vec3(1,0,0));
    WALBERLA_CHECK(intersects(&box5, ray5, t, n));
    WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(2.67157), real_t(1e-4));
@@ -310,7 +310,7 @@ void RaytracerTest(Raytracer::Algorithm raytracingAlgorithm = Raytracer::RAYTRAC
    createSphere(*globalBodyStorage, *forest, storageID, 3, Vec3(4,real_t(5.5),5), real_t(1));
    createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,real_t(8.5),5), real_t(1));
    BoxID box = createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(5,real_t(6.5),5), Vec3(2,4,3));
-   if (box != nullptr) box->rotate(0,math::M_PI/4,math::M_PI/4);
+   if (box != nullptr) box->rotate(0,math::pi/4,math::pi/4);
    createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,1,8), Vec3(2,2,2));
    // Test scene v1 end
    
@@ -318,12 +318,12 @@ void RaytracerTest(Raytracer::Algorithm raytracingAlgorithm = Raytracer::RAYTRAC
    createBox(*globalBodyStorage, *forest, storageID, 9, Vec3(9,9,5), Vec3(1,1,10));
    createCapsule(*globalBodyStorage, *forest, storageID, 10, Vec3(3, 9, 1), real_t(0.5), real_t(7), iron);
    CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, 11, Vec3(7, real_t(3.5), real_t(7.5)), real_t(1), real_t(2), iron);
-   if (capsule != nullptr) capsule->rotate(0,math::M_PI/3,math::M_PI/4-math::M_PI/8);
+   if (capsule != nullptr) capsule->rotate(0,math::pi/3,math::pi/4-math::pi/8);
    // Test scene v2 end
    
    // Test scene v3 additions start
    EllipsoidID ellipsoid = createEllipsoid(*globalBodyStorage, *forest, storageID, 12, Vec3(6,2,real_t(2.5)), Vec3(3,2,real_t(1.2)));
-   ellipsoid->rotate(0, math::M_PI/real_t(6), 0);
+   ellipsoid->rotate(0, math::pi/real_t(6), 0);
    // Test scene v3 end
    
    //raytracer.setTBufferOutputDirectory("tbuffer");
@@ -460,7 +460,7 @@ void HashGridsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t a
       BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), Vec3(len, len, len));
       WALBERLA_CHECK(box_ != nullptr);
       if (boxRotation) {
-         box_->rotate(0, math::realRandom(real_t(0), real_t(1))*math::M_PI, math::realRandom(real_t(0), real_t(1))*math::M_PI);
+         box_->rotate(0, math::realRandom(real_t(0), real_t(1))*math::pi, math::realRandom(real_t(0), real_t(1))*math::pi);
       }
       bodies.push_back(box_);
       bodySIDs.push_back(box_->getSystemID());
@@ -475,7 +475,7 @@ void HashGridsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t a
       walberla::id_t id = walberla::id_t(boxes+i);
       CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), radius, len);
       WALBERLA_CHECK(capsule != nullptr);
-      capsule->rotate(0, math::realRandom(real_t(0), real_t(1))*math::M_PI, math::realRandom(real_t(0), real_t(1))*math::M_PI);
+      capsule->rotate(0, math::realRandom(real_t(0), real_t(1))*math::pi, math::realRandom(real_t(0), real_t(1))*math::pi);
       bodies.push_back(capsule);
       bodySIDs.push_back(capsule->getSystemID());
    }
diff --git a/tests/pe/RigidBody.cpp b/tests/pe/RigidBody.cpp
index 083c6ab43dba6b7cf9fa639b57668529e316599e..1492f75d00b1415a71690dbcbaf3c72d3f0b05db 100644
--- a/tests/pe/RigidBody.cpp
+++ b/tests/pe/RigidBody.cpp
@@ -87,34 +87,34 @@ void checkRotationFunctions()
    auto sp3 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
    auto sp4 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
 
-   sp1->rotate( 1, 0, 0, math::M_PI * real_t(0.5));
-   sp1->rotate( 0, 1, 0, math::M_PI * real_t(0.5));
-   sp1->rotate( 0, 0, 1, math::M_PI * real_t(0.5));
+   sp1->rotate( 1, 0, 0, math::pi * real_t(0.5));
+   sp1->rotate( 0, 1, 0, math::pi * real_t(0.5));
+   sp1->rotate( 0, 0, 1, math::pi * real_t(0.5));
 
-   sp2->rotate( 1, 0, 0, math::M_PI * real_t(0.5));
-   sp2->rotate( 0, 1, 0, math::M_PI * real_t(0.5));
-   sp2->rotate( 0, 0, 1, math::M_PI * real_t(0.5));
+   sp2->rotate( 1, 0, 0, math::pi * real_t(0.5));
+   sp2->rotate( 0, 1, 0, math::pi * real_t(0.5));
+   sp2->rotate( 0, 0, 1, math::pi * real_t(0.5));
 
-   sp3->rotate( math::M_PI * real_t(0.5), math::M_PI * real_t(0.5), math::M_PI * real_t(0.5) );
-   sp4->rotate( Vec3(math::M_PI * real_t(0.5), math::M_PI * real_t(0.5), math::M_PI * real_t(0.5)) );
+   sp3->rotate( math::pi * real_t(0.5), math::pi * real_t(0.5), math::pi * real_t(0.5) );
+   sp4->rotate( Vec3(math::pi * real_t(0.5), math::pi * real_t(0.5), math::pi * real_t(0.5)) );
 
-   WALBERLA_CHECK_FLOAT_EQUAL( sp1->getQuaternion(), Quat(math::M_SQRT2 * real_t(0.5), 0, math::M_SQRT2 * real_t(0.5), 0) );
-   WALBERLA_CHECK_FLOAT_EQUAL( sp2->getQuaternion(), Quat(math::M_SQRT2 * real_t(0.5), 0, math::M_SQRT2 * real_t(0.5), 0) );
-   WALBERLA_CHECK_FLOAT_EQUAL( sp3->getQuaternion(), Quat(math::M_SQRT2 * real_t(0.5), 0, math::M_SQRT2 * real_t(0.5), 0) );
-   WALBERLA_CHECK_FLOAT_EQUAL( sp4->getQuaternion(), Quat(math::M_SQRT2 * real_t(0.5), 0, math::M_SQRT2 * real_t(0.5), 0) );
+   WALBERLA_CHECK_FLOAT_EQUAL( sp1->getQuaternion(), Quat(math::root_two * real_t(0.5), 0, math::root_two * real_t(0.5), 0) );
+   WALBERLA_CHECK_FLOAT_EQUAL( sp2->getQuaternion(), Quat(math::root_two * real_t(0.5), 0, math::root_two * real_t(0.5), 0) );
+   WALBERLA_CHECK_FLOAT_EQUAL( sp3->getQuaternion(), Quat(math::root_two * real_t(0.5), 0, math::root_two * real_t(0.5), 0) );
+   WALBERLA_CHECK_FLOAT_EQUAL( sp4->getQuaternion(), Quat(math::root_two * real_t(0.5), 0, math::root_two * real_t(0.5), 0) );
 
    WALBERLA_CHECK_FLOAT_EQUAL( sp1->getPosition(), Vec3(0, 0, 0) );
    WALBERLA_CHECK_FLOAT_EQUAL( sp2->getPosition(), Vec3(0, 0, 0) );
    WALBERLA_CHECK_FLOAT_EQUAL( sp3->getPosition(), Vec3(0, 0, 0) );
    WALBERLA_CHECK_FLOAT_EQUAL( sp4->getPosition(), Vec3(0, 0, 0) );
 
-   sp1->rotateAroundPoint( Vec3(-10, 0, 0), Vec3(0, 1, 0), math::M_PI );
+   sp1->rotateAroundPoint( Vec3(-10, 0, 0), Vec3(0, 1, 0), math::pi );
    WALBERLA_CHECK_FLOAT_EQUAL( sp1->getPosition(), Vec3(-20, 0, 0) );
 
-   sp2->rotateAroundPoint( Vec3(-10, 0, 0), Vec3(0, 0, 1), math::M_PI );
+   sp2->rotateAroundPoint( Vec3(-10, 0, 0), Vec3(0, 0, 1), math::pi );
    WALBERLA_CHECK_FLOAT_EQUAL( sp2->getPosition(), Vec3(-20, 0, 0) );
 
-   sp3->rotateAroundPoint( Vec3(0, 10, 0), Vec3(0, 0, 1), math::M_PI );
+   sp3->rotateAroundPoint( Vec3(0, 10, 0), Vec3(0, 0, 1), math::pi );
    WALBERLA_CHECK_FLOAT_EQUAL( sp3->getPosition(), Vec3(0, 20, 0) );
 }
 
@@ -155,9 +155,9 @@ int main( int argc, char** argv )
    Vec3 v0 = Vec3(-1,-1,1);
 
    for (auto it = storage.begin(); it != storage.end(); ++it){
-      it->rotateAroundPoint(Vec3(1,1,0), Vec3(0,0,1), math::M_PI);
+      it->rotateAroundPoint(Vec3(1,1,0), Vec3(0,0,1), math::pi);
       WALBERLA_CHECK_FLOAT_EQUAL(it->getPosition(), Vec3(2,2,0));
-      it->rotateAroundOrigin(Vec3(0,0,1), real_c(0.5) * math::M_PI);
+      it->rotateAroundOrigin(Vec3(0,0,1), real_c(0.5) * math::pi);
       WALBERLA_CHECK_FLOAT_EQUAL(it->getPosition(), x0);
       it->setLinearVel(v0);
    }
diff --git a/tests/pe/Union.cpp b/tests/pe/Union.cpp
index de959af081196a78e64d369beb32aab9bb1e61dd..4977a4f980f55c3d3e9eeddb8cbcd47b7e5f320f 100644
--- a/tests/pe/Union.cpp
+++ b/tests/pe/Union.cpp
@@ -126,7 +126,7 @@ void UnionConstruction()
    WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)));
 
    // Check mass volume and inertia
-   WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(4./3.)*real_t(math::M_PI)));
+   WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(4./3.)*real_t(math::pi)));
    WALBERLA_CHECK_FLOAT_EQUAL(un->getMass(), un->getVolume()*Material::getDensity(iron));
    real_t scalar_inertia = real_t(0.4)*un->getMass(); // for sphere: I = 2/5*m*r*r
    WALBERLA_CHECK_EQUAL(un->getInertia(), Mat3(scalar_inertia,0,0,0,scalar_inertia,0,0,0,scalar_inertia));
@@ -143,10 +143,10 @@ void UnionConstruction()
    WALBERLA_CHECK_FLOAT_EQUAL(subb2.getPosition(), Vec3(-1,0,0));
 
    // Check mass volume and inertia
-   WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(8./3.)*real_t(math::M_PI)));
+   WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(8./3.)*real_t(math::pi)));
    WALBERLA_CHECK_FLOAT_EQUAL(un->getMass(), un->getVolume()*Material::getDensity(iron));
    // Mass of one sphere
-   real_t masssphere = real_t(4./3.)*real_t(math::M_PI)*Material::getDensity(iron);
+   real_t masssphere = real_t(4./3.)*real_t(math::pi)*Material::getDensity(iron);
    Mat3 bodyinertia(real_t(2.0)*scalar_inertia, 0, 0, 0, real_t(2.0)*(scalar_inertia + masssphere),0, 0, 0, real_t(2.0)*(scalar_inertia + masssphere));
    WALBERLA_CHECK_FLOAT_EQUAL(un->getInertia(), bodyinertia);
 
@@ -156,7 +156,7 @@ void UnionConstruction()
 
    WALBERLA_LOG_INFO("- Performing Rotation.");
    //Check values for rotated union
-   Quat rotz30(Vec3(0,0,1), real_t(math::M_PI/6.0)); // rotate by 30 deg via z axis
+   Quat rotz30(Vec3(0,0,1), real_t(math::pi/6.0)); // rotate by 30 deg via z axis
    real_t sin30 = real_t(0.5);
    real_t cos30 = real_t(sqrt(3.0)/2.0);
    un->setOrientation(rotz30);
@@ -259,7 +259,7 @@ void UnionAABB() {
    aabb = un->getAABB();
    checkAABB(aabb, AABB(real_t(8),real_t(9),real_t(9),real_t(12),real_t(11),real_t(11)));
 
-   Quat rotz30(Vec3(0,0,1), real_t(math::M_PI/6.0)); // rotate by 30 deg via z axis
+   Quat rotz30(Vec3(0,0,1), real_t(math::pi/6.0)); // rotate by 30 deg via z axis
    real_t sin30 = real_t(0.5);
    real_t cos30 = real_t(sqrt(3.0)/2.0);
    un->setOrientation(rotz30);
diff --git a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
index 29b3d0ccab8083d85aba0077b28df6ed7a8c0c22..81350061ffd987857518de3efec925df2fbdd4bc 100644
--- a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
@@ -247,7 +247,7 @@ uint_t createSpheresRandomly( StructuredBlockForest & forest, pe::BodyStorage &
 {
    real_t domainVolume = generationDomain.volume();
    real_t totalSphereVolume = domainVolume * solidVolumeFraction;
-   real_t sphereVolume = diameter * diameter * diameter * math::M_PI / real_t(6);
+   real_t sphereVolume = diameter * diameter * diameter * math::pi / real_t(6);
    uint_t numberOfSpheres = uint_c( totalSphereVolume / sphereVolume );
 
    real_t xParticle = real_t(0);
@@ -286,7 +286,7 @@ uint_t createSphereLattice( StructuredBlockForest & forest, pe::BodyStorage & gl
                             const BlockDataID & bodyStorageID, const AABB & generationDomain,
                             real_t diameter, real_t solidVolumeFraction, pe::MaterialID & material, real_t initialZVelocity )
 {
-   real_t sphereVolume = math::M_PI * diameter * diameter * diameter / real_t(6);
+   real_t sphereVolume = math::pi * diameter * diameter * diameter / real_t(6);
    real_t numSpheresDesired = solidVolumeFraction * generationDomain.volume() / sphereVolume;
    uint_t spheresPerDirection = uint_c(std::cbrt(numSpheresDesired) );
 
@@ -722,7 +722,7 @@ int main( int argc, char **argv )
    if( solidVolumeFraction < real_t(1e-10) )
    {
       // create only a single sphere
-      solidVolumeFraction = math::M_PI * diameter * diameter * diameter / ( real_t(6) * real_c( xlength * ylength * zlength));
+      solidVolumeFraction = math::pi * diameter * diameter * diameter / ( real_t(6) * real_c( xlength * ylength * zlength));
    }
 
    const real_t diameter_SI = real_t(0.00035); // m, Finn et al, Tab 5
@@ -737,7 +737,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::M_PI / real_t(6);
+   real_t gravitationalForce = - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::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
@@ -826,15 +826,15 @@ int main( int argc, char **argv )
    const real_t restitutionCoeff = real_t(0.88);
    const real_t frictionCoeff = real_t(0.25);
 
-   real_t sphereVolume = diameter * diameter * diameter * math::M_PI / real_t(6);
+   real_t sphereVolume = diameter * diameter * diameter * math::pi / real_t(6);
    const real_t particleMass = densityRatio * sphereVolume;
    const real_t Mij = particleMass * particleMass / ( real_t(2) * particleMass );
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
    const real_t collisionTime = real_t(0.5);
-   const real_t stiffnessCoeff = math::M_PI * math::M_PI * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ) ) );
+   const real_t stiffnessCoeff = math::pi * math::pi * Mij / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ) ) );
    const real_t dampingCoeff = - real_t(2) * std::sqrt( Mij * stiffnessCoeff ) *
-                               ( std::log(restitutionCoeff) / std::sqrt( math::M_PI * math::M_PI + (std::log(restitutionCoeff) * std::log(restitutionCoeff) ) ) );
-   const real_t contactDuration = real_t(2) * math::M_PI * Mij / ( std::sqrt( real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff )); //formula from Uhlman
+                               ( std::log(restitutionCoeff) / std::sqrt( math::pi * math::pi + (std::log(restitutionCoeff) * std::log(restitutionCoeff) ) ) );
+   const real_t contactDuration = real_t(2) * math::pi * Mij / ( std::sqrt( real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff )); //formula from Uhlman
 
    if( !funcTest ) {
       WALBERLA_LOG_INFO_ON_ROOT("Created sediment material with:\n"
@@ -901,7 +901,7 @@ int main( int argc, char **argv )
 
 
    const real_t domainVolume = real_c( xlength * ylength * zlength );
-   real_t actualSolidVolumeFraction = real_c( numSpheres ) * diameter * diameter * diameter * math::M_PI / ( real_t(6) * domainVolume );
+   real_t actualSolidVolumeFraction = real_c( numSpheres ) * diameter * diameter * diameter * math::pi / ( real_t(6) * domainVolume );
    real_t ReynoldsNumber = std::fabs(velUnhindered) * diameter / viscosity;
 
    // apply external forcing on fluid to approximately balance the force from the settling particles to avoid too large fluid or particle velocities
@@ -1487,9 +1487,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::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) ) );
+               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) ) );
                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 )
                {
@@ -1520,7 +1520,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::M_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::pi / real_t(6) );
    WALBERLA_LOG_INFO_ON_ROOT("simulating " << timesteps << " time steps" );
 
 
@@ -1641,9 +1641,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::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::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::TimeStep( blocks, bodyStorageID, *cr, syncCall, dtInteractionSubCycle, peSubSteps, lubricationEvaluationFunction ), "Pe Time Step" );
 
       timeloop.add() << Sweep( DummySweep(), "Dummy Sweep ")
diff --git a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
index 88b5950fdd88a6e0b34538507bb4a5b44cb8b2b5..2ae19a3bc5d34ced1119dfc4686d7718b27b87c6 100644
--- a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
@@ -617,7 +617,7 @@ int main( int argc, char **argv )
    const real_t viscosity = ug * diameter / GalileoNumber;
    const real_t tau = real_t(1) / lbm::collision_model::omegaFromViscosity(viscosity);
 
-   const real_t sphereVolume = math::M_PI * diameter * diameter * diameter / real_t(6);
+   const real_t sphereVolume = math::pi * diameter * diameter * diameter / real_t(6);
    Vector3<real_t> gravitationalForce ( real_t(0), real_t(0), ( densityRatio - real_t(1) ) * sphereVolume * gravity );
 
    if( !funcTest )
@@ -710,13 +710,13 @@ int main( int argc, char **argv )
    const real_t particleMass = densityRatio * sphereVolume;
    const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
    const real_t lnDryResCoeff = std::log(restitutionCoeff);
-   const real_t stiffnessCoeff = math::M_PI * math::M_PI * Mij / ( scaledCollisionTime * scaledCollisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::M_PI * math::M_PI + lnDryResCoeff* lnDryResCoeff ) ) );
+   const real_t stiffnessCoeff = math::pi * math::pi * Mij / ( scaledCollisionTime * scaledCollisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff ) ) );
    const real_t normalizedStiffnessCoeff = stiffnessCoeff / ( ( densityRatio - real_t(1) ) * gravity * sphereVolume / diameter );
 
    const real_t dampingCoeff = - real_t(2) * std::sqrt( Mij * stiffnessCoeff ) *
-                               ( std::log(restitutionCoeff) / std::sqrt( math::M_PI * math::M_PI + (std::log(restitutionCoeff) * std::log(restitutionCoeff) ) ) );
-   const real_t contactDuration = real_t(2) * math::M_PI * Mij / ( std::sqrt( real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff )); //formula from Uhlman
-   const real_t contactDuration2 = std::sqrt(( math::M_PI * math::M_PI + std::log(restitutionCoeff) * std::log(restitutionCoeff)) / ( stiffnessCoeff / Mij)); //formula from Finn
+                               ( std::log(restitutionCoeff) / std::sqrt( math::pi * math::pi + (std::log(restitutionCoeff) * std::log(restitutionCoeff) ) ) );
+   const real_t contactDuration = real_t(2) * math::pi * Mij / ( std::sqrt( real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff )); //formula from Uhlman
+   const real_t contactDuration2 = std::sqrt(( math::pi * math::pi + std::log(restitutionCoeff) * std::log(restitutionCoeff)) / ( stiffnessCoeff / Mij)); //formula from Finn
 
    if( !funcTest )
    {
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
index 3a12ec962315f268e05d2ce59e304d8c465eebd4..b86a57b5f2986c47ee9a15546a89f94be5733db2 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
@@ -137,7 +137,7 @@ public:
                   const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID, real_t sphereRadius) :
          blocks_( blocks ), bodyStorageID_( bodyStorageID ), globalBodyStorage_( globalBodyStorage ),
          boundaryHandlingID_( boundaryHandlingID ), bodyFieldID_( bodyFieldID ),
-         sphereVolume_( math::M_PI * real_t(4) / real_t(3) * sphereRadius * sphereRadius * sphereRadius )
+         sphereVolume_( math::pi * real_t(4) / real_t(3) * sphereRadius * sphereRadius * sphereRadius )
    { }
 
    // check the mapping in the inner domain of the block and check mapped volume against real sphere volume
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index 2f6ac710af30d329c4f6f1b1b4790c30b3a7f55b..5999df24dd04fd3cfe88138d641a3af255a6ef6c 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -225,9 +225,9 @@ public:
       real_t uBar = computeAverageVel();
 
       // f_total = f_drag + f_buoyancy
-      real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::pi * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::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 c60e11186ff6f378f71f73bcc0fbda8e64efdf3e..f77f8076fb7be877189dea7c7ac819d704005535 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -296,9 +296,9 @@ class ForceEval
          real_t uBar = getAverageVel();
 
          // f_total = f_drag + f_buoyancy
-         real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+         real_t totalForce = forceX  + real_c(4.0/3.0) * math::pi * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-         real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
+         real_t normalizedDragForce = totalForce / real_c( 6.0 * math::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 f3b768189dc9f5912161855b75f200077b582f72..be0d9c4d369d0105fab22a2f2e2df053ed5d5b87 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -259,10 +259,10 @@ private:
             WALBERLA_CHECK_FLOAT_EQUAL ( forceSphr2[0], -forceSphr1[0] );
 
             // according to the formula from Ding & Aidun 2003
-            // F = 3/2 * M_PI * rho_L * nu_L * relative velocity of both spheres * r * r * 1/gap
+            // F = 3/2 * math::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::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 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 relErr     = std::fabs( analytical - forceSphr2[0] ) / analytical * real_c(100.0);
             WALBERLA_CHECK_LESS( relErr, real_t(1) );
          }
@@ -355,10 +355,10 @@ private:
          if ( timestep == uint_t(26399) )
          {
             // according to the formula from Ding & Aidun 2003
-            // F = 6 * M_PI * rho_L * nu_L * relative velocity of both bodies=relative velocity of the sphere * r * r * 1/gap
+            // F = 6 * math::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::M_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::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/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 0042c361a7d0209d8d4df3f17fd607d6e9b3b34f..0c90a90d7ecd823585d0e628b1b2b82da8cb1028 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -729,7 +729,7 @@ int main( int argc, char **argv )
    // The buoyancy force on the body due to this pressure gradient has to be added 'artificially'
    // F_{buoyancy} = - V_{body} * grad ( p ) = V_{body} * \rho_{fluid} *  a
    // ( V_{body} := volume of body,  a := acceleration driving the flow )
-   Vector3<real_t> buoyancyForce(math::M_PI / real_t(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
+   Vector3<real_t> buoyancyForce(math::pi / real_t(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
                                  real_t(0), real_t(0));
    timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, buoyancyForce ), "Buoyancy force" );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
index f8a4e3c731f8aabb7ad58253e437dc9f63b3b206..6b2d0c5517f533c0be3b47c03b1129f90fe42c4b 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
@@ -407,7 +407,7 @@ int main( int argc, char **argv )
                                      uint_c(floor(domainSize_SI[1] / dx_SI + real_t(0.5)) ),
                                      uint_c(floor(domainSize_SI[2] / dx_SI + real_t(0.5)) ) );
    const real_t diameter = diameter_SI / dx_SI;
-   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+   const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
 
    const real_t expectedSettlingVelocity = real_t(0.01);
    const real_t dt_SI = expectedSettlingVelocity / expectedSettlingVelocity_SI * dx_SI;
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
index 19fc95b569fe61cd8a6442ee35c2f0ce5d62f939..5bcb1a9efe4e9e39156c2546e03682d4e79e89dc 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
@@ -515,7 +515,7 @@ int main( int argc, char **argv )
                                      uint_c(floor(domainSize_SI[1] / dx_SI + real_t(0.5)) ),
                                      uint_c(floor(domainSize_SI[2] / dx_SI + real_t(0.5)) ) );
    const real_t diameter = diameter_SI / dx_SI;
-   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+   const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
 
    const real_t expectedSettlingVelocity = real_t(0.01);
    const real_t dt_SI = expectedSettlingVelocity / expectedSettlingVelocity_SI * dx_SI;
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
index e7bc2538a4e784b5032a1a60b6d7fc9ac06027fc..84fcc93be0798cfd92cd3726161ccd76e73b76fe 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
@@ -476,7 +476,7 @@ int main( int argc, char **argv )
                                      uint_c(floor(domainSize_SI[1] / dx_SI + real_t(0.5)) ),
                                      uint_c(floor(domainSize_SI[2] / dx_SI + real_t(0.5)) ) );
    const real_t diameter = diameter_SI / dx_SI;
-   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+   const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
 
    const real_t expectedSettlingVelocity = real_t(0.01);
    const real_t dt_SI = expectedSettlingVelocity / expectedSettlingVelocity_SI * dx_SI;
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index 0b649234430aa2820acae8544694316445f0590a..ee405e7b449d899909ddac558a413fa6f1e52615 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -374,7 +374,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::M_PI * setup.radius * setup.radius * setup.radius /
+   setup.phi    = real_c(4.0/3.0) * math::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 969b81d6274883096db4df87531957f7b86f9d01..a2042350550ef9d957890c79138076a5c8a629c6 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::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::pi * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::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 b35c41db12a57bcea104e4b6fa81213ca0ca3f14..115bcaffab070ba8c5d7c82290bd107ec2574177 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
@@ -296,9 +296,9 @@ public:
       real_t uBar = computeAverageVel();
 
       // f_total = f_drag + f_buoyancy
-      real_t totalForce = forceX  + real_c(4.0/3.0) * math::M_PI * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
+      real_t totalForce = forceX  + real_c(4.0/3.0) * math::pi * setup_->radius * setup_->radius * setup_->radius * setup_->extForce ;
 
-      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::M_PI * setup_->visc * setup_->radius * uBar );
+      real_t normalizedDragForce = totalForce / real_c( 6.0 * math::pi * setup_->visc * setup_->radius * uBar );
 
       // update drag force values
       normalizedDragOld_ = normalizedDragNew_;
diff --git a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
index fbec373004928d25d7b1d548417e32a2d58f221b..38ae149c9bb0394eee46418a634af575190f3544 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
@@ -623,7 +623,7 @@ int main( int argc, char **argv )
    // The buoyancy force on the body due to this pressure gradient has to be added 'artificially'
    // F_{buoyancy} = - V_{body} * grad ( p ) = V_{body} * \rho_{fluid} *  a
    // ( V_{body} := volume of body,  a := acceleration driving the flow )
-   Vector3<real_t> buoyancyForce(math::M_PI / real_t(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
+   Vector3<real_t> buoyancyForce(math::pi / real_t(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
                                  real_t(0), real_t(0));
    timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, buoyancyForce ), "Buoyancy force" );
 
diff --git a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
index 60a0e97be6fedc0da74921190a3b2f1330f73a9f..1a2e16a8937ef95f9f3037415679e9bb0069c01f 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::M_PI * setup.radius * setup.radius * setup.radius /
+   setup.phi    = real_c(4.0/3.0) * math::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