diff --git a/benchmarks/hollow_sphere/hollow_sphere.cc b/benchmarks/hollow_sphere/hollow_sphere.cc
index 6fb3e70dac7..1aca2e0483b 100644
--- a/benchmarks/hollow_sphere/hollow_sphere.cc
+++ b/benchmarks/hollow_sphere/hollow_sphere.cc
@@ -73,17 +73,17 @@ namespace aspect
 
         if (mmm == -1)
           {
-            alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
-            beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
+            alpha=-gammma*(Utilities::fixed_power<3>(R2)-Utilities::fixed_power<3>(R1))/(Utilities::fixed_power<3>(R2)*std::log(R1)-Utilities::fixed_power<3>(R1)*std::log(R2));
+            beta=-3*gammma*(std::log(R2)-std::log(R1))/(Utilities::fixed_power<3>(R1)*std::log(R2)-Utilities::fixed_power<3>(R2)*std::log(R1)) ;
             fr=alpha/(r*r)+beta*r;
-            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
+            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*Utilities::fixed_power<3>(r)+gammma);
           }
         else
           {
-            alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+            alpha=gammma*(mmm+1)*(Utilities::fixed_power<-3>(R1)-Utilities::fixed_power<-3>(R2))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
             beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
             fr=alpha/std::pow(r,mmm+3)+beta*r;
-            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
+            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*Utilities::fixed_power<3>(r)+gammma);
           }
 
         const double v_r    =gr*std::cos(theta);
@@ -115,17 +115,17 @@ namespace aspect
         if (mmm == -1)
           {
             mur=mu0;
-            alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
-            beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
-            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
+            alpha=-gammma*(Utilities::fixed_power<3>(R2)-Utilities::fixed_power<3>(R1))/(Utilities::fixed_power<3>(R2)*std::log(R1)-Utilities::fixed_power<3>(R1)*std::log(R2));
+            beta=-3*gammma*(std::log(R2)-std::log(R1))/(Utilities::fixed_power<3>(R1)*std::log(R2)-Utilities::fixed_power<3>(R2)*std::log(R1)) ;
+            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*Utilities::fixed_power<3>(r)+gammma);
             hr=2/r*gr*mur;
           }
         else
           {
             mur=mu0*std::pow(r,mmm+1);
-            alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+            alpha=gammma*(mmm+1)*(Utilities::fixed_power<-3>(R1)-Utilities::fixed_power<-3>(R2))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
             beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
-            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
+            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*Utilities::fixed_power<3>(r)+gammma);
             hr=(mmm+3)/r*gr*mur;
           }
 
@@ -147,17 +147,17 @@ namespace aspect
 
         if (mmm == -1)
           {
-            alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
-            beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
+            alpha=-gammma*(Utilities::fixed_power<3>(R2)-Utilities::fixed_power<3>(R1))/(Utilities::fixed_power<3>(R2)*std::log(R1)-Utilities::fixed_power<3>(R1)*std::log(R2));
+            beta=-3*gammma*(std::log(R2)-std::log(R1))/(Utilities::fixed_power<3>(R1)*std::log(R2)-Utilities::fixed_power<3>(R2)*std::log(R1)) ;
             fr=alpha/(r*r)+beta*r;
-            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
+            gr=-2/(r*r)*(alpha*std::log(r)+beta/3*Utilities::fixed_power<3>(r)+gammma);
           }
         else
           {
-            alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+            alpha=gammma*(mmm+1)*(Utilities::fixed_power<-3>(R1)-Utilities::fixed_power<-3>(R2))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
             beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
             fr=alpha/std::pow(r,mmm+3)+beta*r;
-            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
+            gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*Utilities::fixed_power<3>(r)+gammma);
           }
 
         return -(6.*gr + 4.*fr) * std::cos(theta) * mu0 / r;
@@ -380,15 +380,15 @@ namespace aspect
 
           if (mmm == -1)
             {
-              alpha = -gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
-              beta  = -3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
-              rho = -(alpha/std::pow(r,4)*(8*std::log(r)-6) + 8./3.*beta/r+8*gammma/std::pow(r,4))*std::cos(theta) + rho_0;
+              alpha = -gammma*(Utilities::fixed_power<3>(R2)-Utilities::fixed_power<3>(R1))/(Utilities::fixed_power<3>(R2)*std::log(R1)-Utilities::fixed_power<3>(R1)*std::log(R2));
+              beta  = -3*gammma*(std::log(R2)-std::log(R1))/(Utilities::fixed_power<3>(R1)*std::log(R2)-Utilities::fixed_power<3>(R2)*std::log(R1)) ;
+              rho = -(alpha/Utilities::fixed_power<4>(r)*(8*std::log(r)-6) + 8./3.*beta/r+8*gammma/Utilities::fixed_power<4>(r))*std::cos(theta) + rho_0;
             }
           else
             {
-              alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+              alpha=gammma*(mmm+1)*(Utilities::fixed_power<-3>(R1)-Utilities::fixed_power<-3>(R2))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
               beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
-              rho= -(2*alpha*std::pow(r,-4)*(mmm+3)/(mmm+1)*(mmm-1)-2*beta/3*(mmm-1)*(mmm+3)*std::pow(r,mmm)-mmm*(mmm+5)*2*gammma*std::pow(r,mmm-3) )*std::cos(theta) + rho_0;
+              rho= -(2*alpha*Utilities::fixed_power<-4>(r)*(mmm+3)/(mmm+1)*(mmm-1)-2*beta/3*(mmm-1)*(mmm+3)*std::pow(r,mmm)-mmm*(mmm+5)*2*gammma*std::pow(r,mmm-3) )*std::cos(theta) + rho_0;
             }
 
           out.densities[i] = rho;
diff --git a/benchmarks/layeredflow/layeredflow.cc b/benchmarks/layeredflow/layeredflow.cc
index e27a592945e..135972f1d08 100644
--- a/benchmarks/layeredflow/layeredflow.cc
+++ b/benchmarks/layeredflow/layeredflow.cc
@@ -50,8 +50,8 @@ namespace aspect
         const double yplus=(1.+y0)/beta;
         const double yminus=(1.-y0)/beta;
         const double C1 = 2.*numbers::PI /
-                          ( beta*std::log(beta*beta+std::pow(1.+y0,2.0))-2.*(1.+y0)*std::atan(yplus)
-                            - beta*std::log(beta*beta+std::pow(1.-y0,2.0))+2.*(1.-y0)*std::atan(yminus)
+                          ( beta*std::log(beta*beta+Utilities::fixed_power<2>(1.+y0))-2.*(1.+y0)*std::atan(yplus)
+                            - beta*std::log(beta*beta+Utilities::fixed_power<2>(1.-y0))+2.*(1.-y0)*std::atan(yminus)
                             + 2.*numbers::PI*(1.+2.*epsilon) );
 
         const double C2 = ( beta*std::log( beta*beta+Utilities::fixed_power<2>(1+y0) )- 2.*(1.+y0)*std::atan(yplus)
diff --git a/benchmarks/viscosity_grooves/viscosity_grooves.cc b/benchmarks/viscosity_grooves/viscosity_grooves.cc
index dd8aa5bf288..a6d007740ea 100644
--- a/benchmarks/viscosity_grooves/viscosity_grooves.cc
+++ b/benchmarks/viscosity_grooves/viscosity_grooves.cc
@@ -62,7 +62,7 @@ namespace aspect
         const GeometryModel::Box<2> *geometry
           = dynamic_cast<const GeometryModel::Box<2>*> (&geometry_model);
         const double L=geometry->get_extents()[0];
-        return x*x*y*y+x*y+5. - std::pow(L,4.)/9.-std::pow(L,2.)/4.-5.;
+        return x*x*y*y+x*y+5. - Utilities::fixed_power<4>(L)/9.-Utilities::fixed_power<2>(L)/4.-5.;
       }
 
       double
diff --git a/source/boundary_temperature/dynamic_core.cc b/source/boundary_temperature/dynamic_core.cc
index 593a9f41803..cb1e888daf5 100644
--- a/source/boundary_temperature/dynamic_core.cc
+++ b/source/boundary_temperature/dynamic_core.cc
@@ -562,7 +562,7 @@ namespace aspect
     DynamicCore<dim>::
     get_mass(const double r) const
     {
-      return 4.*numbers::PI*Rho_cen*(-std::pow(L,2)/2.*r*std::exp(-std::pow(r/L,2))+std::pow(L,3)/4.*std::sqrt(numbers::PI)*std::erf(r/L));
+      return 4.*numbers::PI*Rho_cen*(-Utilities::fixed_power<2>(L)/2.*r*std::exp(-Utilities::fixed_power<2>(r/L))+Utilities::fixed_power<3>(L)/4.*std::sqrt(numbers::PI)*std::erf(r/L));
     }
 
 
@@ -591,9 +591,9 @@ namespace aspect
     DynamicCore<dim>::
     get_pressure(const double r) const
     {
-      return P_CMB-(4*numbers::PI*constants::big_g*std::pow(Rho_cen,2))/3
-             *((3*std::pow(r,2)/10.-std::pow(L,2)/5)*std::exp(-std::pow(r/L,2))
-               -(3*std::pow(Rc,2)/10-std::pow(L,2)/5)*std::exp(-std::pow(Rc/L,2)));
+      return P_CMB-(4*numbers::PI*constants::big_g*Utilities::fixed_power<2>(Rho_cen))/3
+             *((3*Utilities::fixed_power<2>(r)/10.-Utilities::fixed_power<2>(L)/5)*std::exp(-Utilities::fixed_power<2>(r/L))
+               -(3*Utilities::fixed_power<2>(Rc)/10-Utilities::fixed_power<2>(L)/5)*std::exp(-Utilities::fixed_power<2>(Rc/L)));
     }
 
 
@@ -603,7 +603,7 @@ namespace aspect
     DynamicCore<dim>::
     get_rho(const double r) const
     {
-      return Rho_cen*std::exp(-std::pow(r/L,2));
+      return Rho_cen*std::exp(-Utilities::fixed_power<2>(r/L));
     }
 
 
@@ -613,7 +613,7 @@ namespace aspect
     DynamicCore<dim>::
     get_g(const double r) const
     {
-      return (4*numbers::PI/3)*constants::big_g*Rho_cen*r*(1-3*std::pow(r,2)/(5*std::pow(L,2)));
+      return (4*numbers::PI/3)*constants::big_g*Rho_cen*r*(1-3*Utilities::fixed_power<2>(r)/(5*Utilities::fixed_power<2>(L)));
     }
 
 
@@ -623,7 +623,7 @@ namespace aspect
     DynamicCore<dim>::
     get_T(const double Tc, const double r) const
     {
-      return Tc*std::exp((std::pow(Rc,2)-std::pow(r,2))/std::pow(D,2));
+      return Tc*std::exp((Utilities::fixed_power<2>(Rc)-Utilities::fixed_power<2>(r))/Utilities::fixed_power<2>(D));
     }
 
 
@@ -633,8 +633,8 @@ namespace aspect
     DynamicCore<dim>::
     get_gravity_potential(const double r) const
     {
-      return 2./3.*numbers::PI*constants::big_g*Rho_cen*(std::pow(r,2)*(1.-3.*std::pow(r,2)
-                                                                        /(10.*std::pow(L,2)))-std::pow(Rc,2)*(1.-3.*std::pow(Rc,2)/(10.*std::pow(L,2))));
+      return 2./3.*numbers::PI*constants::big_g*Rho_cen*(Utilities::fixed_power<2>(r)*(1.-3.*Utilities::fixed_power<2>(r)
+                                                         /(10.*Utilities::fixed_power<2>(L)))-Utilities::fixed_power<2>(Rc)*(1.-3.*Utilities::fixed_power<2>(Rc)/(10.*Utilities::fixed_power<2>(L))));
     }
 
 
@@ -644,8 +644,8 @@ namespace aspect
     DynamicCore<dim>::
     get_specific_heating(const double Tc, double &Qs,double &Es) const
     {
-      const double A = std::sqrt(1./(std::pow(L,-2)+std::pow(D,-2)));
-      const double Is = 4.*numbers::PI*get_T(Tc,0.)*Rho_cen*(-std::pow(A,2)*Rc/2.*std::exp(-std::pow(Rc/A,2))+std::pow(A,3)*std::sqrt(numbers::PI)/4.*std::erf(Rc/A));
+      const double A = std::sqrt(1./(Utilities::fixed_power<-2>(L)+Utilities::fixed_power<-2>(D)));
+      const double Is = 4.*numbers::PI*get_T(Tc,0.)*Rho_cen*(-Utilities::fixed_power<2>(A)*Rc/2.*std::exp(-Utilities::fixed_power<2>(Rc/A))+Utilities::fixed_power<3>(A)*std::sqrt(numbers::PI)/4.*std::erf(Rc/A));
 
       Qs = -Cp/Tc*Is;
       Es = Cp/Tc*(Mc-Is/Tc);
@@ -661,13 +661,13 @@ namespace aspect
       double It = numbers::signaling_nan<double>();
       if (D>L)
         {
-          const double B = std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,2)));
-          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
+          const double B = std::sqrt(1/(1/Utilities::fixed_power<2>(L)-1/Utilities::fixed_power<2>(D)));
+          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-Utilities::fixed_power<2>(B)*Rc/2*std::exp(-Utilities::fixed_power<2>(Rc/B))+Utilities::fixed_power<3>(B)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
         }
       else
         {
-          const double B = std::sqrt(1/(std::pow(D,-2)-std::pow(L,-2)));
-          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2);
+          const double B = std::sqrt(1/(Utilities::fixed_power<-2>(D)-Utilities::fixed_power<-2>(L)));
+          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*Rc/2*std::exp(Utilities::fixed_power<2>(Rc/B))-Utilities::fixed_power<2>(B)*fun_Sn(B,Rc,100)/2);
         }
 
       Qr = Mc*core_data.H;
@@ -684,17 +684,17 @@ namespace aspect
       double It = numbers::signaling_nan<double>();
       if (D>L)
         {
-          const double B = std::sqrt(1./(1./std::pow(L,2)-1./std::pow(D,2)));
-          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
-          It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*r/2*std::exp(-std::pow(r/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(r/B));
+          const double B = std::sqrt(1./(1./Utilities::fixed_power<2>(L)-1./Utilities::fixed_power<2>(D)));
+          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-Utilities::fixed_power<2>(B)*Rc/2*std::exp(-Utilities::fixed_power<2>(Rc/B))+Utilities::fixed_power<3>(B)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
+          It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-Utilities::fixed_power<2>(B)*r/2*std::exp(-Utilities::fixed_power<2>(r/B))+Utilities::fixed_power<3>(B)/std::sqrt(numbers::PI)/4*std::erf(r/B));
         }
       else
         {
-          const double B = std::sqrt(1./(std::pow(D,-2)-std::pow(L,-2)));
-          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2);
-          It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*r/2*std::exp(std::pow(r/B,2))-std::pow(B,2)*fun_Sn(B,r,100)/2);
+          const double B = std::sqrt(1./(Utilities::fixed_power<-2>(D)-Utilities::fixed_power<-2>(L)));
+          It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*Rc/2*std::exp(Utilities::fixed_power<2>(Rc/B))-Utilities::fixed_power<2>(B)*fun_Sn(B,Rc,100)/2);
+          It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*r/2*std::exp(Utilities::fixed_power<2>(r/B))-Utilities::fixed_power<2>(B)*fun_Sn(B,r,100)/2);
         }
-      const double Cc = 4*numbers::PI*std::pow(r,2)*get_rho(r)*X/(Mc-get_mass(r));
+      const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
       Eh = Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc;
     }
 
@@ -705,17 +705,17 @@ namespace aspect
     DynamicCore<dim>::
     get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const
     {
-      const double Cc = 4*numbers::PI*std::pow(r,2)*get_rho(r)*X/(Mc-get_mass(r));
-      const double C_2 = 3./16.*std::pow(L,2) - 0.5*std::pow(Rc,2)*(1.-3./10.*std::pow(Rc/L,2));
+      const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
+      const double C_2 = 3./16.*Utilities::fixed_power<2>(L) - 0.5*Utilities::fixed_power<2>(Rc)*(1.-3./10.*Utilities::fixed_power<2>(Rc/L));
       if (r==Rc)
         Qg = 0.;
       else
         {
-          Qg = (8./3.*std::pow(numbers::PI*Rho_cen,2)*constants::big_g*(
-                  ((3./20.*std::pow(Rc,5)-std::pow(L,2)*std::pow(Rc,3)/8.-C_2*std::pow(L,2)*Rc)*std::exp(-std::pow(Rc/L,2))
-                   +C_2/2.*std::pow(L,3)*std::sqrt(numbers::PI)*std::erf(Rc/L))
-                  -((3./20.*std::pow(r,5)-std::pow(L,2)*std::pow(r,3)/8.-C_2*std::pow(L,2)*r)*std::exp(-std::pow(r/L,2))
-                    +C_2/2.*std::pow(L,3)*std::sqrt(numbers::PI)*std::erf(r/L)))
+          Qg = (8./3.*Utilities::fixed_power<2>(numbers::PI*Rho_cen)*constants::big_g*(
+                  ((3./20.*Utilities::fixed_power<5>(Rc)-Utilities::fixed_power<2>(L)*Utilities::fixed_power<3>(Rc)/8.-C_2*Utilities::fixed_power<2>(L)*Rc)*std::exp(-Utilities::fixed_power<2>(Rc/L))
+                   +C_2/2.*Utilities::fixed_power<3>(L)*std::sqrt(numbers::PI)*std::erf(Rc/L))
+                  -((3./20.*Utilities::fixed_power<5>(r)-Utilities::fixed_power<2>(L)*Utilities::fixed_power<3>(r)/8.-C_2*Utilities::fixed_power<2>(L)*r)*std::exp(-Utilities::fixed_power<2>(r/L))
+                    +C_2/2.*Utilities::fixed_power<3>(L)*std::sqrt(numbers::PI)*std::erf(r/L)))
                 -(Mc-get_mass(r))*get_gravity_potential(r))*Beta_c*Cc;
         }
 
@@ -729,8 +729,8 @@ namespace aspect
     DynamicCore<dim>::
     get_adiabatic_heating(const double Tc, double &Ek, double &Qk) const
     {
-      Ek = 16*numbers::PI*k_c*std::pow(Rc,5)/5/std::pow(D,4);
-      Qk = 8*numbers::PI*std::pow(Rc,3)*k_c*Tc/std::pow(D,2);
+      Ek = 16*numbers::PI*k_c*Utilities::fixed_power<5>(Rc)/5/Utilities::fixed_power<4>(D);
+      Qk = 8*numbers::PI*Utilities::fixed_power<3>(Rc)*k_c*Tc/Utilities::fixed_power<2>(D);
     }
 
 
@@ -740,7 +740,7 @@ namespace aspect
     DynamicCore<dim>::
     get_latent_heating(const double Tc, const double r, double &El, double &Ql) const
     {
-      Ql = 4.*numbers::PI*std::pow(r,2)*Lh*get_rho(r);
+      Ql = 4.*numbers::PI*Utilities::fixed_power<2>(r)*Lh*get_rho(r);
       El = Ql*(get_T(Tc,r)-Tc)/(Tc*get_T(Tc,r));
     }
 
diff --git a/source/geometry_model/ellipsoidal_chunk.cc b/source/geometry_model/ellipsoidal_chunk.cc
index 2cc2c0afc0f..288b061acd2 100644
--- a/source/geometry_model/ellipsoidal_chunk.cc
+++ b/source/geometry_model/ellipsoidal_chunk.cc
@@ -495,7 +495,7 @@ namespace aspect
           bottom_depth = prm.get_double("Depth");
           semi_major_axis_a = prm.get_double("Semi-major axis");
           eccentricity = prm.get_double("Eccentricity");
-          semi_minor_axis_b = std::sqrt((1 - std::pow(eccentricity,2.)) * std::pow(semi_major_axis_a,2.));
+          semi_minor_axis_b = std::sqrt((1 - Utilities::fixed_power<2>(eccentricity)) * Utilities::fixed_power<2>(semi_major_axis_a));
           EW_subdiv = prm.get_integer("East-West subdivisions");
           NS_subdiv = prm.get_integer("North-South subdivisions");
           depth_subdiv = prm.get_integer("Depth subdivisions");
diff --git a/source/material_model/latent_heat_melt.cc b/source/material_model/latent_heat_melt.cc
index 8ef32f7352f..f51aba907d6 100644
--- a/source/material_model/latent_heat_melt.cc
+++ b/source/material_model/latent_heat_melt.cc
@@ -172,7 +172,7 @@ namespace aspect
             = beta * std::pow((temperature - T_solidus)/(T_lherz_liquidus - T_solidus),beta-1)
               * (dT_solidus_dp * (temperature - T_lherz_liquidus)
                  + dT_lherz_liquidus_dp * (T_solidus - temperature))
-              / std::pow(T_lherz_liquidus - T_solidus,2);
+              / Utilities::fixed_power<2>(T_lherz_liquidus - T_solidus);
 
           // melt fraction after melting of all clinopyroxene
           const double R_cpx = r1 + r2 * pressure;
@@ -181,7 +181,7 @@ namespace aspect
           if (peridotite_melt_fraction(temperature, pressure, compositional_fields, position) > F_max)
             {
               const double T_max = std::pow(F_max,1.0/beta) * (T_lherz_liquidus - T_solidus) + T_solidus;
-              const double dF_max_dp = - M_cpx * std::pow(r1 + r2 * pressure,-2) * r2;
+              const double dF_max_dp = - M_cpx * Utilities::fixed_power<-2>(r1 + r2 * pressure) * r2;
               const double dT_max_dp = dT_solidus_dp
                                        + 1.0/beta * std::pow(F_max,1.0/beta - 1.0) * dF_max_dp * (T_lherz_liquidus - T_solidus)
                                        + std::pow(F_max,1.0/beta) * (dT_lherz_liquidus_dp - dT_solidus_dp);
diff --git a/source/material_model/melt_boukare.cc b/source/material_model/melt_boukare.cc
index 2b9b335a569..e7db6da2276 100644
--- a/source/material_model/melt_boukare.cc
+++ b/source/material_model/melt_boukare.cc
@@ -190,7 +190,7 @@ namespace aspect
                                                       (a + (1. - a) * std::pow(1 + b * (pressure - Pth), c)));
 
           const double Cp_ref = reference_specific_heats[i] + specific_heat_linear_coefficients[i] * in.temperature[q]
-                                + specific_heat_second_coefficients[i] * std::pow(in.temperature[q], -2.)
+                                + specific_heat_second_coefficients[i] * Utilities::fixed_power<-2>(in.temperature[q])
                                 + specific_heat_third_coefficients[i] * std::pow(in.temperature[q], -0.5);
 
           const long double dSdT0 = reference_volumes[i] * reference_bulk_moduli[i] * Utilities::fixed_power<2>(heat_capacity_ratio * reference_thermal_expansivities[i])
diff --git a/source/material_model/reaction_model/katz2003_mantle_melting.cc b/source/material_model/reaction_model/katz2003_mantle_melting.cc
index d7a61eb8b7b..2b80c681c11 100644
--- a/source/material_model/reaction_model/katz2003_mantle_melting.cc
+++ b/source/material_model/reaction_model/katz2003_mantle_melting.cc
@@ -129,7 +129,7 @@ namespace aspect
               = beta * std::pow((temperature - T_solidus)/(T_lherz_liquidus - T_solidus),beta-1)
                 * (dT_solidus_dp * (temperature - T_lherz_liquidus)
                    + dT_lherz_liquidus_dp * (T_solidus - temperature))
-                / std::pow(T_lherz_liquidus - T_solidus,2);
+                / Utilities::fixed_power<2>(T_lherz_liquidus - T_solidus);
 
             // melt fraction after melting of all clinopyroxene
             const double R_cpx = r1 + r2 * std::max(0.0, pressure);
@@ -138,7 +138,7 @@ namespace aspect
             if (melt_fraction(temperature, pressure) > F_max)
               {
                 const double T_max = std::pow(F_max,1.0/beta) * (T_lherz_liquidus - T_solidus) + T_solidus;
-                const double dF_max_dp = - M_cpx * std::pow(r1 + r2 * pressure,-2) * r2;
+                const double dF_max_dp = - M_cpx * Utilities::fixed_power<-2>(r1 + r2 * pressure) * r2;
                 const double dT_max_dp = dT_solidus_dp
                                          + 1.0/beta * std::pow(F_max,1.0/beta - 1.0) * dF_max_dp * (T_lherz_liquidus - T_solidus)
                                          + std::pow(F_max,1.0/beta) * (dT_lherz_liquidus_dp - dT_solidus_dp);
diff --git a/source/volume_of_fluid/utilities.cc b/source/volume_of_fluid/utilities.cc
index a95a2b5e8fa..05630f02a2e 100644
--- a/source/volume_of_fluid/utilities.cc
+++ b/source/volume_of_fluid/utilities.cc
@@ -404,7 +404,7 @@ namespace aspect
             //Trapezoid Y
             coeffs[0]=(d + 0.5*n_xp)/n_xp; // 1
             coeffs[1]=-0.5*n_yp/n_xp*sign_n_y; // 2*y - 1
-            coeffs[2]=0.25*(12.0*std::pow(d, 2) - 3.0*n_xp*n_xp + n_yp*n_yp)/(n_xp*n_xp)*sign_n_x; // 2*x - 1
+            coeffs[2]=0.25*(12.0*dealii::Utilities::fixed_power<2>(d) - 3.0*n_xp*n_xp + n_yp*n_yp)/(n_xp*n_xp)*sign_n_x; // 2*x - 1
             coeffs[3]=-3.0*d*n_yp/(n_xp*n_xp)*sign_n_x*sign_n_y; // (2*x - 1)*(2*y - 1)
           }
         else
diff --git a/tests/ellipsoidal_chunk_geometry.cc b/tests/ellipsoidal_chunk_geometry.cc
index cd14a0627d2..2e44fe3d28b 100644
--- a/tests/ellipsoidal_chunk_geometry.cc
+++ b/tests/ellipsoidal_chunk_geometry.cc
@@ -84,7 +84,7 @@ int f()
 
   const double semi_major_axis_a = 6378137.0;
   const double eccentricity = 8.1819190842622e-2;
-  const double semi_minor = std::sqrt((1 - std::pow(eccentricity,2)) * std::pow(semi_major_axis_a,2));
+  const double semi_minor = std::sqrt((1 - dealii::Utilities::fixed_power<2>(eccentricity)) * dealii::Utilities::fixed_power<2>(semi_major_axis_a));
   {
     std::cout << "WGS84 test" << std::endl;