Commit a395aa6c authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

Make LinOpacity shorter as it will be inline (and yes, it;s just aesthetic).

refs #597
parent a08b08a4
......@@ -447,6 +447,82 @@ namespace fargOCA {
/// \brief When not provided by user
real computeDensityStartWithAlphaViscosity(real userAccretionRate, real viscosity, real aspectRatio);
// That stuff need to be inline (sort of) to please CUDA (sort of)
namespace details {
KOKKOS_INLINE_FUNCTION
real
LinOpacity(real temp, real rho)
{
using std::pow;
using std::cbrt;
real const t234 = 1.6e3;
real const t456 = 5.7e3;
real const t678 = 2.28e6;
real ts4 = 1.e-4*temp; // to avoid overflow
if (temp < t234*pow(rho,4.44444444e-2)){
// coefficients for opacity laws in cgs units.
real const ak[] = { 2.e-4, 2.e16, 5.e-3};
// disjoint opacity laws
real o1 = ak[0]*pow(temp, 2);
real o2 = ak[1]*temp/pow(temp, 8);
real o3 = ak[2]*temp;
// parameters used for smoothing
real o1an = pow(o1, 2);
real o2an = pow(o2, 2);
// smoothed and continuous opacity law for regions 1, 2, and 3.
return pow((pow(o1an*o2an/(o1an+o2an),2)
+pow(o3/(1+1.e22/pow(temp,10)),4)),
0.25);
} else if (temp < t456*pow(rho,2.381e-2)) {
// coefficients for opacity laws T_4 units.
real const bk[] = { 50, 2.e-2, 2.e4 };
// disjoint opacity laws for 3, 4, and 5.
real o4 = bk[1]*pow(cbrt(rho),2)/pow(ts4,9);
real o5 = bk[2]*pow(cbrt(rho),2)*pow(ts4,3);
// parameters used for smoothing
real o4an = pow(o4,4);
real o3an = pow(bk[0]*ts4,4);
// smoothed and continuous opacity law for regions 3, 4, and 5.
return pow(((o4an*o3an/(o4an+o3an))
+pow(o5/(1+6.561e-5/pow(ts4,8)),4)),
0.25);
} else if (temp < t678*pow(rho,2.267e-1) || rho <= 1.0e-10){
// coefficients for opacity laws 3, 4, 5, 6, 7, and 8 in T_4 units.
real const bk[] { 2.e4, 1.e4, 1.5e10 };
// disjoint opacity laws for 5, 6, and 7.
real o5 = bk[0]*pow(cbrt(rho),2)*pow(ts4,3);
real o6 = bk[1]*cbrt(rho)*pow(ts4,10);
real o7 = bk[2]*rho/(pow(ts4, 2)*sqrt(ts4));
// parameters used for smoothing
real o6an = pow(o6,2);
real o7an = pow(o7,2);
return pow((pow(o6an*o7an/(o6an+o7an), 2)
+pow(o5/(1+pow(ts4/(1.1*pow(rho,0.04762)),10)),4)),
0.25);
} else {
// no scattering! (So undefined behavior is ok ?)
return std::numeric_limits<real>::infinity();
}
}
}
KOKKOS_INLINE_FUNCTION
real
DiskPhysic::Adiabatic::Radiative::kappa(real temp, real rho) const {
real kappa;
switch (opacityLaw) {
case OL_BELL_LIN:
kappa = details::LinOpacity(temp,rho);
break;
case OL_CONSTANT:
kappa = constKappa;
break;
}
return kappa*RHO0*R0+1.0e-300;
}
}
#endif
......@@ -754,98 +754,6 @@ namespace fargOCA {
h5::writeAttribute(grp, "shadow_angle", shadowAngle);
}
namespace details {
KOKKOS_INLINE_FUNCTION
real
oplin(real temp, real rho)
{
using std::pow;
using std::cbrt;
// data t234,t456,t678/1.6e3, 5.7e3, 2.28e6/
real t234 = 1.6e3;
real t456 = 5.7e3;
real t678 = 2.28e6;
// coefficients for opacity laws 1, 2, and 3 in cgs units.
// data ak1,ak2,ak3/2.e-4, 2.e16, 5.e-3/
real ak1 = 2.e-4;
real ak2 = 2.e16;
real ak3 = 5.e-3;
/* coefficients for opacity laws 3, 4, 5, 6, 7, and 8 in T_4 units.
data bk3,bk4,bk5,bk6,bk7 50., 2.e-2, 2.e4, 1.d4, 1.5d10 */
real bk3 = 50.;
real bk4 = 2.e-2;
real bk5 = 2.e4;
real bk6 = 1.e4;
real bk7 = 1.5e10;
/* test T against (T_23 * T_34 * T_34)**0.333333333 */
real opacity = std::numeric_limits<real>::quiet_NaN();
if( temp < t234*pow(rho,4.44444444e-2)){
// disjoint opacity laws
real o1 = ak1*pow(temp, 2);
real o2 = ak2*temp/pow(temp, 8);
real o3 = ak3*temp;
// parameters used for smoothing
real o1an = pow(o1, 2);
real o2an = pow(o2, 2);
// smoothed and continuous opacity law for regions 1, 2, and 3.
opacity = pow((pow(o1an*o2an/(o1an+o2an),2)
+pow(o3/(1+1.e22/pow(temp,10)),4)),
0.25);
} else if ( temp < t456*pow(rho,2.381e-2)){
// to avoid overflow
real ts4 = 1.e-4*temp;
// disjoint opacity laws for 3, 4, and 5.
real o4 = bk4*pow(cbrt(rho),2)/pow(ts4,9);
real o5 = bk5*pow(cbrt(rho),2)*pow(ts4,3);
// parameters used for smoothing
real o4an = pow(o4,4);
real o3an = pow(bk3*ts4,4);
// smoothed and continuous opacity law for regions 3, 4, and 5.
opacity = pow(((o4an*o3an/(o4an+o3an))
+pow(o5/(1+6.561e-5/pow(ts4,8)),4)),
0.25);
} else if (temp < t678*pow(rho,2.267e-1) || rho <= 1.0e-10){
// to avoid overflow
real ts4 = 1.e-4*temp;
// disjoint opacity laws for 5, 6, and 7.
real o5 = bk5*pow(cbrt(rho),2)*pow(ts4,3);
real o6 = bk6*cbrt(rho)*pow(ts4,10);
real o7 = bk7*rho/(pow(ts4, 2)*sqrt(ts4));
// parameters used for smoothing
real o6an = pow(o6,2);
real o7an = pow(o7,2);
// smoothed and continuous opacity law for regions 5, 6, and 7.
/* FortranForm
opacty=((o6an*o7an/(o6an+o7an))**2+(o5/(1+(ts4/(1.1*rho**0.04762))**10)**4)**0.25 */
opacity = pow(pow(o6an*o7an/(o6an+o7an), 2)+pow(o5/(1+pow(ts4/(1.1*pow(rho,0.04762)),10)),4.0),0.25);
} else {
// no scattering! (So undefined behavior is ok ?)
opacity = std::numeric_limits<real>::infinity();
}
return opacity;
}
}
KOKKOS_FUNCTION
real
DiskPhysic::Adiabatic::Radiative::kappa(real temp, real rho) const {
real kappa;
switch (opacityLaw) {
case OL_BELL_LIN:
kappa = details::oplin(temp,rho);
break;
case OL_CONSTANT:
kappa = constKappa;
break;
}
return kappa*RHO0*R0+1.0e-300;
}
void
DiskPhysic::Adiabatic::Cooling::readH5(H5::Group const& grp) {
h5::readAttribute(grp, "time_scale", timeScale);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment