Skip to content

Commit

Permalink
Replace std::pow where possible
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Dec 12, 2024
1 parent 7d5bac3 commit cc1452a
Show file tree
Hide file tree
Showing 10 changed files with 63 additions and 63 deletions.
40 changes: 20 additions & 20 deletions benchmarks/hollow_sphere/hollow_sphere.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
}

Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 2 additions & 2 deletions benchmarks/layeredflow/layeredflow.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion benchmarks/viscosity_grooves/viscosity_grooves.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
64 changes: 32 additions & 32 deletions source/boundary_temperature/dynamic_core.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}


Expand Down Expand Up @@ -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)));
}


Expand All @@ -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));
}


Expand All @@ -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)));
}


Expand All @@ -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));
}


Expand All @@ -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))));
}


Expand All @@ -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);
Expand All @@ -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;
Expand All @@ -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;
}

Expand All @@ -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;
}

Expand All @@ -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);
}


Expand All @@ -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));
}

Expand Down
2 changes: 1 addition & 1 deletion source/geometry_model/ellipsoidal_chunk.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down
4 changes: 2 additions & 2 deletions source/material_model/latent_heat_melt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion source/material_model/melt_boukare.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down
Loading

0 comments on commit cc1452a

Please sign in to comment.