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;