Commit 73a6bf5d authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

Merge branch '597-compile-with-cuda' into 'master'

Resolve "Compile with CUDA"

Closes #597

See merge request DISC/fargOCA!471
parents d73ee7f6 c3f5e5c6
......@@ -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
......@@ -56,6 +56,7 @@ namespace fargOCA {
HillAreas inner; /**< Force arising from the inner disk */
HillAreas outer; /**< Force arising from the outer disk */
KOKKOS_INLINE_FUNCTION
HillForce() : inner(), outer() {}
HillForce(Disk const& disk, Planet const& planet);
HillForce(HillForce&& other) = default;
......
......@@ -51,13 +51,20 @@ namespace fargOCA {
std::string myODir;
int myNbSteps;
int myStep;
/// \brief The manipulated disk.
Disk& disk() { return *myDisk; }
Disk const& disk() const { return *myDisk; }
SimulationPrivate(shptr<Disk> disk,
shptr<TrackingConfig const> trackingConfig,
std::string odir);
SimulationPrivate(SimulationPrivate&& simu);
void velocityStep(ScalarField const& systemPotential, real dt);
void adiabaticStep(real dt);
void viscosityStep(real dt);
// Cuda want's plain functions.
void ctor();
};
......@@ -65,6 +72,7 @@ namespace fargOCA {
/// \brief Drive the evolution of a planetary disk through time.
class Simulation : private SimulationPrivate {
public:
typedef SimulationPrivate super;
template<class T> using opt = std::optional<T>;
public:
......@@ -81,16 +89,16 @@ namespace fargOCA {
Simulation(Simulation&& simu);
~Simulation();
/// \brief The manipulated disk.
Disk& disk() { return super::disk(); }
Disk const& disk() const { return super::disk(); }
/// \brief Initiate the simulation.
void start();
/// \brief Execute one time step of the simulaton.
void doStep();
/// \brief The manipulated disk.
Disk& disk() { return *myDisk; }
Disk const& disk() const { return *myDisk; }
CompoundTracker& tracker() { return *myTracker; }
CompoundTracker const& tracker() const { return *myTracker; }
......@@ -123,8 +131,6 @@ namespace fargOCA {
std::string odir() const { return myODir; };
private:
void velocityStep(ScalarField const& systemPotential, real dt);
void viscosityStep(real dt);
void adiabaticStep(real dt) { SimulationPrivate::adiabaticStep(dt); }
void energyStep(real dt);
......
......@@ -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);
......
......@@ -98,7 +98,7 @@ namespace fargOCA {
int const nrMax = grid.managedRadiusMax();
kk::View<HillForce*> vf("local Hill contribs", nrMax-nrMin);
for(int i= 0; i < vf.stride(0); ++i) { vf(i) = HillForce(); }
kk::deep_copy(vf, HillForce());
real nan = std::numeric_limits<real>::quiet_NaN();
int const ns = grid.nbSector();
int const ni = grid.nbLayer();
......
......@@ -154,9 +154,10 @@ namespace fargOCA {
arr1d<real const> radii = grid.radii();
arr1d<real const> phiMed = grid.phiMed();
arr3d<real> stellarrad = stellarRadiation;
arr3d<real const> dens = disk.density().data();
arr3d<real const> tau = opticalDepth(disk);
arr3d<real> stellarrad = stellarRadiation;
arr3d<real const> dens = disk.density().data();
arr3d<real const> tau = opticalDepth(disk);
arr3d<real const> simpleOpacity = this->simpleOpacity; // avoid this capture
if (disk.physic().adiabatic->radiative->star) {
real fStar = disk.physic().adiabatic->radiative->star->fStar();
......@@ -275,11 +276,16 @@ namespace fargOCA {
arr3d<real const> temperature = disk.temperature()->data();
arr3d<real const> density = disk.density().data();
DiskPhysic::Adiabatic::Radiative cfg = *disk.physic().adiabatic->radiative;
// avoid this capture:
arr3d<real> simpleOpacity = this->simpleOpacity;
arr3d<real> RosselandOpacity = this->RosselandOpacity;
arr3d<real> PlankOpacity = this->PlankOpacity;
kk::parallel_for("simple_opacity", grid.fullRange(),
LBD(int i, int h, int j) {
simpleOpacity(i,h,j) = cfg.dustToGas*100*3.5*RHO0*R0;
});
kk::parallel_for("Rosse_opacity", grid.fullRange(),
kk::parallel_for("Rosseland_opacity", grid.fullRange(),
LBD(int i, int h, int j) {
// multiplication by a factor 100 of DUSTTOGAS as opacity is already calculated for
// a dust to gas ratio of 0.01
......
......@@ -320,7 +320,7 @@ namespace fargOCA {
}
void
Simulation::velocityStep(ScalarField const& systemPotential, real dt) {
SimulationPrivate::velocityStep(ScalarField const& systemPotential, real dt) {
GridDispatch const& grid = disk().grid();
int const nr = grid.nbRadius();
......@@ -424,7 +424,7 @@ namespace fargOCA {
}
void
Simulation::viscosityStep(real dt)
SimulationPrivate::viscosityStep(real dt)
{
GridDispatch const& grid = disk().grid();
int const nr = grid.nbRadius();
......
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