Merge pull request #6055 from gassmoeller/avoid_pow

Replace std::pow where possible

benchmarks/burstedde/burstedde.cc

@@ -360,20 +360,20 @@ namespace aspect
 
       Tensor<1,dim> g;
 
-      g[0]=((y*z+3.*std::pow(x,2)*std::pow(y,3)*z)- mu*(2.+6.*x*y))
-           -dmudx*(2.+4.*x+2.*y+6.*std::pow(x,2)*y)
-           -dmudy*(x+std::pow(x,3)+y+2.*x*std::pow(y,2))
+      g[0]=((y*z+3.*Utilities::fixed_power<2>(x)*Utilities::fixed_power<3>(y)*z)- mu*(2.+6.*x*y))
+           -dmudx*(2.+4.*x+2.*y+6.*Utilities::fixed_power<2>(x)*y)
+           -dmudy*(x+Utilities::fixed_power<3>(x)+y+2.*x*Utilities::fixed_power<2>(y))
            -dmudz*(-3.*z-10.*x*y*z);
 
-      g[1]=((x*z+3.*std::pow(x,3)*std::pow(y,2)*z)- mu*(2.+2.*std::pow(x,2)+2.*std::pow(y,2)))
-           -dmudx*(x+std::pow(x,3)+y+2.*x*std::pow(y,2))
-           -dmudy*(2.+2.*x+4.*y+4.*std::pow(x,2)*y)
-           -dmudz*(-3.*z-5.*std::pow(x,2)*z);
+      g[1]=((x*z+3.*Utilities::fixed_power<3>(x)*Utilities::fixed_power<2>(y)*z)- mu*(2.+2.*Utilities::fixed_power<2>(x)+2.*Utilities::fixed_power<2>(y)))
+           -dmudx*(x+Utilities::fixed_power<3>(x)+y+2.*x*Utilities::fixed_power<2>(y))
+           -dmudy*(2.+2.*x+4.*y+4.*Utilities::fixed_power<2>(x)*y)
+           -dmudz*(-3.*z-5.*Utilities::fixed_power<2>(x)*z);
 
-      g[2]=((x*y+std::pow(x,3)*std::pow(y,3)) - mu*(-10.*y*z))
+      g[2]=((x*y+Utilities::fixed_power<3>(x)*Utilities::fixed_power<3>(y)) - mu*(-10.*y*z))
            -dmudx*(-3.*z-10.*x*y*z)
-           -dmudy*(-3.*z-5.*std::pow(x,2)*z)
-           -dmudz*(-4.-6.*x-6.*y-10.*std::pow(x,2)*y);
+           -dmudy*(-3.*z-5.*Utilities::fixed_power<2>(x)*z)
+           -dmudz*(-4.-6.*x-6.*y-10.*Utilities::fixed_power<2>(x)*y);
 
       return g;
     }

benchmarks/compressibility_formulations/plugins/compressibility_formulations.cc

@@ -155,7 +155,7 @@ namespace aspect
           if (this->get_parameters().formulation_mass_conservation ==
               Parameters<dim>::Formulation::MassConservation::isentropic_compression)
             {
-              out.compressibilities[i] = reference_compressibility - std::pow(thermal_expansivity, 2) * in.temperature[i]
+              out.compressibilities[i] = reference_compressibility - Utilities::fixed_power<2>(thermal_expansivity) * in.temperature[i]
                                          / (out.densities[i] * out.specific_heat[i]);
             }
         }

benchmarks/layeredflow/layeredflow.cc

@@ -54,10 +54,10 @@ namespace aspect
                             - beta*std::log(beta*beta+pow(1.-y0,2.0))+2.*(1.-y0)*std::atan(yminus)
                             + 2.*numbers::PI*(1.+2.*epsilon) );
 
-        const double C2 = ( beta*std::log( beta*beta+std::pow(1+y0,2.) )- 2.*(1.+y0)*std::atan(yplus)
+        const double C2 = ( beta*std::log( beta*beta+Utilities::fixed_power<2>(1+y0) )- 2.*(1.+y0)*std::atan(yplus)
                             + numbers::PI*(1.+2.*epsilon) )*C1 ;
         const double y = pos[1]-1.0;
-        const double v_x = (-beta*C1*std::log(beta*beta+std::pow(y-y0,2.0))+2.*(y-y0)*C1*std::atan((y-y0)/beta)
+        const double v_x = (-beta*C1*std::log(beta*beta+Utilities::fixed_power<2>(y-y0))+2.*(y-y0)*C1*std::atan((y-y0)/beta)
                             + numbers::PI*(1.+2.*epsilon)*y*C1+C2)/(2.*numbers::PI);
         const double v_y = 0;
         return Point<2> (v_x,v_y);

benchmarks/nsinker/nsinker.cc

@@ -292,7 +292,7 @@ namespace aspect
         {
           double dist = p.distance(Point<2>(centers[s](0), centers[s](1)));
           double temp = 1-std::exp(-delta*
-                                   std::pow(std::max(0.0,dist-omega/2.0),2));
+                                   Utilities::fixed_power<2>(std::max(0.0,dist-omega/2.0)));
           chi *= temp;
         }
       return chi;
@@ -310,7 +310,7 @@ namespace aspect
         {
           double dist = p.distance(centers[s]);
           double temp = 1-std::exp(-delta*
-                                   std::pow(std::max(0.0,dist-omega/2.0),2));
+                                   Utilities::fixed_power<2>(std::max(0.0,dist-omega/2.0)));
           chi *= temp;
         }
       return chi;

benchmarks/nsinker_spherical_shell/nsinker.cc

@@ -311,7 +311,7 @@ namespace aspect
         {
           double dist = p.distance(Point<2>(centers[s](0), centers[s](1)));
           double temp = 1-std::exp(-delta*
-                                   std::pow(std::max(0.0,dist-omega/2.0),2));
+                                   Utilities::fixed_power<2>(std::max(0.0,dist-omega/2.0)));
           chi *= temp;
         }
       return chi;
@@ -329,7 +329,7 @@ namespace aspect
         {
           double dist = p.distance(centers[s]);
           double temp = 1-std::exp(-delta*
-                                   std::pow(std::max(0.0,dist-omega/2.0),2));
+                                   Utilities::fixed_power<2>(std::max(0.0,dist-omega/2.0)));
           chi *= temp;
         }
       return chi;

benchmarks/solitary_wave/solitary_wave.cc

@@ -399,13 +399,13 @@ namespace aspect
 
         double length_scaling (const double porosity) const
         {
-          return std::sqrt(reference_permeability * std::pow(porosity,3) * (xi_0 + 4.0/3.0 * eta_0) / eta_f);
+          return std::sqrt(reference_permeability * Utilities::fixed_power<3>(porosity) * (xi_0 + 4.0/3.0 * eta_0) / eta_f);
         }
 
         double velocity_scaling (const double porosity) const
         {
           const Point<dim> surface_point = this->get_geometry_model().representative_point(0.0);
-          return reference_permeability * std::pow(porosity,2) * (reference_rho_s - reference_rho_f)
+          return reference_permeability * Utilities::fixed_power<2>(porosity) * (reference_rho_s - reference_rho_f)
                  * this->get_gravity_model().gravity_vector(surface_point).norm() / eta_f;
         }
 
@@ -451,7 +451,7 @@ namespace aspect
 
                 melt_out->compaction_viscosities[i] = xi_0 * (1.0 - porosity);
                 melt_out->fluid_viscosities[i]= eta_f;
-                melt_out->permeabilities[i]= reference_permeability * std::pow(porosity,3);
+                melt_out->permeabilities[i]= reference_permeability * Utilities::fixed_power<3>(porosity);
                 melt_out->fluid_densities[i]= reference_rho_f;
                 melt_out->fluid_density_gradients[i] = 0.0;
               }

benchmarks/solubility/plugin/solubility.cc

@@ -172,7 +172,7 @@ namespace aspect
               double porosity = std::max(in.composition[q][porosity_idx],0.0);
 
               fluid_out->fluid_viscosities[q] = eta_f;
-              fluid_out->permeabilities[q] = reference_permeability * std::pow(porosity,3) * std::pow(1.0-porosity,2);
+              fluid_out->permeabilities[q] = reference_permeability * Utilities::fixed_power<3>(porosity) * Utilities::fixed_power<2>(1.0-porosity);
 
               fluid_out->fluid_densities[q] = reference_rho_f * std::exp(fluid_compressibility * (in.pressure[q] - this->get_surface_pressure()));
 
@@ -387,7 +387,7 @@ namespace aspect
     reference_darcy_coefficient () const
     {
       // 0.01 = 1% melt
-      return reference_permeability * std::pow(0.01,3.0) / eta_f;
+      return reference_permeability * Utilities::fixed_power<3>(0.01) / eta_f;
     }
 
 

benchmarks/yamauchi_takei_2016_anelasticity/anelasticity_temperature.cc

@@ -290,7 +290,7 @@ namespace aspect
           isothermal_volume_change=r2.first;
           compressibility=isothermal_volume_change*std::exp((gruneisen_parameter+1)*(std::pow(isothermal_volume_change,-1)-1));
           pressure_dependent_density=reference_density*isothermal_volume_change;
-          integrated_thermal_expansivity=(2.832e-5*(temperature-273))+((0.758e-8/2)*(std::pow(temperature,2)-std::pow(273,2)));
+          integrated_thermal_expansivity=(2.832e-5*(temperature-273))+((0.758e-8/2)*(Utilities::fixed_power<2>(temperature)-Utilities::fixed_power<2>(273)));
           density=pressure_dependent_density*(1-(compressibility*integrated_thermal_expansivity));
         }
       // determine J1 term (real part of complex compliance)
@@ -299,8 +299,8 @@ namespace aspect
                                                                       std::erf((std::log(peak_period/normalized_period))/(std::sqrt(2)*peak_width)))));
       // determine J2 term (imaginary part of complex compliance)
       //double loss_compliance=unrelaxed_compliance*(M_PI/2)*(background_amplitude*(std::pow(normalized_period,background_slope))+
-      // (peak_amplitude*std::exp(-1*(std::pow(std::log(peak_period/normalized_period),2)/
-      // (2*std::pow(peak_width,2))))))+(unrelaxed_compliance*normalized_period);
+      // (peak_amplitude*std::exp(-1*(Utilities::fixed_power<2>(std::log(peak_period/normalized_period))/
+      // (2*Utilities::fixed_power<2>(peak_width))))))+(unrelaxed_compliance*normalized_period);
       // calculate anharmonic Vs
       // anharmonic_Vs=1/(std::sqrt(density*unrelaxed_compliance)*1e3);
       // calculate Vs

cookbooks/inner_core_convection/inner_core_convection.cc

@@ -175,7 +175,7 @@ namespace aspect
 
           // The gravity is zero in the center of the Earth, and we assume the density to be constant and equal to 1.
           // We fix the surface pressure to 0 (and we are only interested in pressure differences anyway).
-          return gravity_magnitude * 0.5 * (1.0 - std::pow(radius/max_radius,2));
+          return gravity_magnitude * 0.5 * (1.0 - Utilities::fixed_power<2>(radius/max_radius));
         }
       else
         return hydrostatic_pressure_profile.value(Point<1>(radius));

cookbooks/morency_doin_2004/morency_doin.cc

@@ -70,7 +70,7 @@ namespace aspect
              */
             const double e2inv = std::sqrt((strain_rate[0][0]*strain_rate[0][0] + strain_rate[0][1]*strain_rate[0][1])/2.0);
 
-            const double edot = std::sqrt(std::pow(e2inv,2) + std::pow(min_strain_rate,2));
+            const double edot = std::sqrt(Utilities::fixed_power<2>(e2inv) + Utilities::fixed_power<2>(min_strain_rate));
 
             // Calculate 1/(v_eff^v) and 1/(v_eff^p)
             double one_over_veffv = 0.0;

source/boundary_traction/initial_lithostatic_pressure.cc

@@ -392,7 +392,7 @@ namespace aspect
       prm.leave_subsection();
 
       // Check that we have enough integration points for this mesh.
-      AssertThrow(std::pow(2.0,refinement) <= n_points, ExcMessage("Not enough integration points for this resolution."));
+      AssertThrow(Utilities::pow(2u,refinement) <= n_points, ExcMessage("Not enough integration points for this resolution."));
 
       pressure.resize(n_points,-1);
     }

source/geometry_model/ellipsoidal_chunk.cc

@@ -124,8 +124,8 @@ namespace aspect
         const double p      = std::sqrt(x(0) * x(0) + x(1) * x(1)); // distance from origin projected onto x-y plane
         const double th     = std::atan2(R * x(2), b * p); // starting guess for theta
         const double phi    = std::atan2(x(1), x(0)); // azimuth (geodetic longitude)
-        const double theta  = std::atan2(x(2) + (R * R - b * b) / b * std::pow(std::sin(th),3),
-                                         (p - (eccentricity * eccentricity * R  * std::pow(std::cos(th),3)))); // first iterate for theta
+        const double theta  = std::atan2(x(2) + (R * R - b * b) / b * Utilities::fixed_power<3>(std::sin(th)),
+                                         (p - (eccentricity * eccentricity * R  * Utilities::fixed_power<3>(std::cos(th))))); // first iterate for theta
         const double R_bar  = R / (std::sqrt(1 - eccentricity * eccentricity * std::sin(theta) * std::sin(theta))); // first iterate for R_bar
 
         Point<3> phi_theta_d;

source/initial_temperature/adiabatic.cc

@@ -224,7 +224,7 @@ namespace aspect
             }
           else
             {
-              const double exponential = -kappa * std::pow(numbers::PI, 2) * age_top / std::pow(lithosphere_thickness, 2);
+              const double exponential = -kappa * Utilities::fixed_power<2>(numbers::PI) * age_top / Utilities::fixed_power<2>(lithosphere_thickness);
               double sum_terms = 0.;
               for (unsigned int n=1; n<11; ++n)
                 {

source/initial_temperature/continental_geotherm.cc

@@ -80,10 +80,10 @@ namespace aspect
 
       // Temperature in layer 1
       if (depth <= thicknesses[0])
-        return -0.5*densities[0]*heat_productivities[0]/conductivities[0]*std::pow(depth,2) + (0.5*densities[0]*heat_productivities[0]*thicknesses[0]/conductivities[0] + (T1-T0)/thicknesses[0])*depth + T0;
+        return -0.5*densities[0]*heat_productivities[0]/conductivities[0]*Utilities::fixed_power<2>(depth) + (0.5*densities[0]*heat_productivities[0]*thicknesses[0]/conductivities[0] + (T1-T0)/thicknesses[0])*depth + T0;
       // Temperature in layer 2
       else if (depth <= thicknesses[0]+thicknesses[1])
-        return -0.5*densities[1]*heat_productivities[1]/conductivities[1]*std::pow(depth-thicknesses[0],2.) + (0.5*densities[1]*heat_productivities[1]*thicknesses[1]/conductivities[1] + (T2-T1)/thicknesses[1])*(depth-thicknesses[0]) + T1;
+        return -0.5*densities[1]*heat_productivities[1]/conductivities[1]*Utilities::fixed_power<2>(depth-thicknesses[0]) + (0.5*densities[1]*heat_productivities[1]*thicknesses[1]/conductivities[1] + (T2-T1)/thicknesses[1])*(depth-thicknesses[0]) + T1;
       // Temperature in layer 3
       else if (depth <= thicknesses[0]+thicknesses[1]+thicknesses[2])
         return (LAB_isotherm-T2)/thicknesses[2] *(depth-thicknesses[0]-thicknesses[1]) + T2;

source/initial_temperature/spherical_shell.cc

@@ -230,9 +230,9 @@ namespace aspect
       const double x = (scale - this->depth)*std::cos(angle);
       const double y = (scale - this->depth)*std::sin(angle);
       const double Perturbation = (sign * amplitude *
-                                   std::exp( -( std::pow((position(0)*scale/R1-x),2)
+                                   std::exp( -( Utilities::fixed_power<2>((position(0)*scale/R1-x))
                                                 +
-                                                std::pow((position(1)*scale/R1-y),2) ) / sigma));
+                                                Utilities::fixed_power<2>((position(1)*scale/R1-y)) ) / sigma));
 
       if (r > R1 - 1e-6*R1 || InterpolVal + Perturbation < T1)
         return T1*dT;

source/material_model/equation_of_state/thermodynamic_table_lookup.cc

@@ -196,8 +196,8 @@ namespace aspect
 
                 for (unsigned int j = 0; j < material_lookup.size(); ++j)
                   {
-                    const double mu = material_lookup[j]->density(in.temperature[i],in.pressure[i])*std::pow(material_lookup[j]->seismic_Vs(in.temperature[i],in.pressure[i]), 2.);
-                    const double k =  material_lookup[j]->density(in.temperature[i],in.pressure[i])*std::pow(material_lookup[j]->seismic_Vp(in.temperature[i],in.pressure[i]), 2.) - 4./3.*mu;
+                    const double mu = material_lookup[j]->density(in.temperature[i],in.pressure[i])*Utilities::fixed_power<2>(material_lookup[j]->seismic_Vs(in.temperature[i],in.pressure[i]));
+                    const double k =  material_lookup[j]->density(in.temperature[i],in.pressure[i])*Utilities::fixed_power<2>(material_lookup[j]->seismic_Vp(in.temperature[i],in.pressure[i])) - 4./3.*mu;
 
                     k_voigt += volume_fractions[i][j] * k;
                     mu_voigt += volume_fractions[i][j] * mu;

source/material_model/latent_heat_melt.cc

@@ -194,7 +194,7 @@ namespace aspect
                 = dF_max_dp
                   - dF_max_dp * std::pow((temperature - T_max)/(T_liquidus - T_max),beta)
                   + (1.0 - F_max) * beta * std::pow((temperature - T_max)/(T_liquidus - T_max),beta-1)
-                  * (dT_max_dp * (T_max - T_liquidus) - (dT_liquidus_dp - dT_max_dp) * (temperature - T_max)) / std::pow(T_liquidus - T_max, 2);
+                  * (dT_max_dp * (T_max - T_liquidus) - (dT_liquidus_dp - dT_max_dp) * (temperature - T_max)) / Utilities::fixed_power<2>(T_liquidus - T_max);
             }
 
           double melt_fraction_derivative = 0;