diff --git a/src/muscl/MHDInitFunctors2D.h b/src/muscl/MHDInitFunctors2D.h
index 7e79b45..b2733e1 100644
--- a/src/muscl/MHDInitFunctors2D.h
+++ b/src/muscl/MHDInitFunctors2D.h
@@ -313,9 +313,8 @@ class InitOrszagTangFunctor2D : public MHDBaseFunctor2D
 
     const real_t gamma0 = params.settings.gamma0;
 
-    const double TwoPi = 4.0 * asin(1.0);
-    const double B0 = 1.0 / sqrt(2.0 * TwoPi);
-    const double p0 = gamma0 / (2.0 * TwoPi);
+    const double B0 = 1.0 / sqrt(4 * PI_v<double>);
+    const double p0 = gamma0 / (4 * PI_v<double>);
     const double d0 = gamma0 * p0;
     const double v0 = 1.0;
 
@@ -332,17 +331,17 @@ class InitOrszagTangFunctor2D : public MHDBaseFunctor2D
       Udata(i, j, ID) = d0;
 
       // rho*vx
-      Udata(i, j, IU) = static_cast<real_t>(-d0 * v0 * sin(yPos * TwoPi));
+      Udata(i, j, IU) = static_cast<real_t>(-d0 * v0 * sin(yPos * TWOPI_F));
 
       // rho*vy
-      Udata(i, j, IV) = static_cast<real_t>(d0 * v0 * sin(xPos * TwoPi));
+      Udata(i, j, IV) = static_cast<real_t>(d0 * v0 * sin(xPos * TWOPI_F));
 
       // rho*vz
       Udata(i, j, IW) = ZERO_F;
 
       // bx, by, bz
-      Udata(i, j, IBX) = -B0 * sin(yPos * TwoPi);
-      Udata(i, j, IBY) = B0 * sin(2.0 * xPos * TwoPi);
+      Udata(i, j, IBX) = -B0 * sin(yPos * TWOPI_F);
+      Udata(i, j, IBY) = B0 * sin(2.0 * xPos * TWOPI_F);
       Udata(i, j, IBZ) = 0.0;
     }
 
@@ -524,7 +523,7 @@ class InitKelvinHelmholtzFunctor2D_MHD : public MHDBaseFunctor2D
 
       const real_t d = rho1 + ramp * (rho2 - rho1);
       const real_t u = v1 + ramp * (v2 - v1);
-      const real_t v = w0 * sin(n * M_PI * x);
+      const real_t v = w0 * sin(n * PI_F * x);
 
       const real_t bx = 0.5;
       const real_t by = 0.0;
@@ -552,7 +551,7 @@ class InitKelvinHelmholtzFunctor2D_MHD : public MHDBaseFunctor2D
 
       const real_t d = (y >= y1 and y <= y2) ? d_in : d_out;
       const real_t u = (y >= y1 and y <= y2) ? vflow_in : vflow_out;
-      const real_t v = w0 * sin(n * M_PI * x);
+      const real_t v = w0 * sin(n * PI_F * x);
 
       const real_t bx = 0.5;
       const real_t by = 0.0;
diff --git a/src/muscl/MHDInitFunctors3D.h b/src/muscl/MHDInitFunctors3D.h
index 560c4e9..c84ba2c 100644
--- a/src/muscl/MHDInitFunctors3D.h
+++ b/src/muscl/MHDInitFunctors3D.h
@@ -347,9 +347,8 @@ class InitOrszagTangFunctor3D : public MHDBaseFunctor3D
 
     const real_t gamma0 = params.settings.gamma0;
 
-    const double TwoPi = 4.0 * asin(1.0);
-    const double B0 = 1.0 / sqrt(2.0 * TwoPi);
-    const double p0 = gamma0 / (2.0 * TwoPi);
+    const double B0 = 1.0 / sqrt(4 * PI_v<double>);
+    const double p0 = gamma0 / (4 * PI_v<double>);
     const double d0 = gamma0 * p0;
     const double v0 = 1.0;
     const double kt = otParams.kt;
@@ -365,19 +364,19 @@ class InitOrszagTangFunctor3D : public MHDBaseFunctor3D
     Udata(i, j, k, ID) = d0;
 
     // rho*vx
-    Udata(i, j, k, IU) = static_cast<real_t>(-d0 * v0 * sin(yPos * TwoPi));
+    Udata(i, j, k, IU) = static_cast<real_t>(-d0 * v0 * sin(yPos * TWOPI_F));
 
     // rho*vy
-    Udata(i, j, k, IV) = static_cast<real_t>(d0 * v0 * sin(xPos * TwoPi));
+    Udata(i, j, k, IV) = static_cast<real_t>(d0 * v0 * sin(xPos * TWOPI_F));
 
     // rho*vz
     Udata(i, j, k, IW) = ZERO_F;
 
     // bx, by, bz
     Udata(i, j, k, IBX) =
-      -B0 * cos(2 * TwoPi * kt * (zPos - zmin) / (zmax - zmin)) * sin(yPos * TwoPi);
+      -B0 * cos(2 * TWOPI_F * kt * (zPos - zmin) / (zmax - zmin)) * sin(yPos * TWOPI_F);
     Udata(i, j, k, IBY) =
-      B0 * cos(2 * TwoPi * kt * (zPos - zmin) / (zmax - zmin)) * sin(2.0 * xPos * TwoPi);
+      B0 * cos(2 * TWOPI_F * kt * (zPos - zmin) / (zmax - zmin)) * sin(2.0 * xPos * TWOPI_F);
     Udata(i, j, k, IBZ) = 0.0;
 
   } // init_all_var_but_energy
@@ -394,9 +393,8 @@ class InitOrszagTangFunctor3D : public MHDBaseFunctor3D
 
     const real_t gamma0 = params.settings.gamma0;
 
-    const double TwoPi = 4.0 * asin(1.0);
-    // const double B0    = 1.0/sqrt(2.0*TwoPi);
-    const double p0 = gamma0 / (2.0 * TwoPi);
+    // const double B0    = 1.0/sqrt(2.0*TWOPI_F);
+    const double p0 = gamma0 / (4 * PI_v<double>);
     // const double d0    = gamma0*p0;
     // const double v0    = 1.0;
 
diff --git a/src/shared/real_type.h b/src/shared/real_type.h
index 3e659df..a1aa965 100644
--- a/src/shared/real_type.h
+++ b/src/shared/real_type.h
@@ -63,6 +63,7 @@ KOKKOS_IMPL_MATH_CONSTANT(HALF,       0.500000000000000000000000000000000000L);
 KOKKOS_IMPL_MATH_CONSTANT(ONE,        1.000000000000000000000000000000000000L);
 KOKKOS_IMPL_MATH_CONSTANT(TWO,        2.000000000000000000000000000000000000L);
 KOKKOS_IMPL_MATH_CONSTANT(ONE_FOURTH, 0.250000000000000000000000000000000000L);
+KOKKOS_IMPL_MATH_CONSTANT(PI, 3.141592653589793238462643383279502884L);
 
 #undef KOKKOS_IMPL_MATH_CONSTANT
 
@@ -71,6 +72,8 @@ constexpr auto HALF_F = HALF_v<real_t>;
 constexpr auto ONE_F  = ONE_v<real_t>;
 constexpr auto TWO_F  = TWO_v<real_t>;
 constexpr auto ONE_FOURTH_F  = ONE_FOURTH_v<real_t>;
+constexpr auto PI_F = PI_v<real_t>;
+constexpr auto TWOPI_F = 2 * PI_v<real_t>;
 
 // math function
 #if defined(USE_DOUBLE) ||  defined(USE_MIXED_PRECISION)