Skip to content

Commit

Permalink
Merge pull request #80 from MChabanov/fishbone-moncrief-id
Browse files Browse the repository at this point in the history 
Fishbone moncrief
  • Loading branch information
@MChabanov
MChabanov authored Dec 12, 2024
2 parents f81841f + caa627f commit ec368fc
Show file tree
Hide file tree
Showing 6 changed files with 102 additions and 88 deletions.
2 changes: 1 addition & 1 deletion FishboneMoncriefIDX/configuration.ccl
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
REQUIRES Loop
REQUIRES Loop AsterUtils
2 changes: 1 addition & 1 deletion FishboneMoncriefIDX/interface.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,5 +2,5 @@ IMPLEMENTS: FishboneMoncriefIDX

INHERITS: ADMBaseX HydroBaseX AsterX

#USES INCLUDE HEADER: fixmath.hxx
USES INCLUDE HEADER: loop_device.hxx
USES INCLUDE HEADER: aster_utils.hxx
106 changes: 73 additions & 33 deletions FishboneMoncriefIDX/src/FM_disk_ID.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,32 @@
#include <cstdio>
#include <cstdbool>
#include <cmath>
#include <cstdlib> // Needed for rand()

#include "FM_disk_implementation.hxx"
#include "FM_disk_utils.hxx"

namespace FMdisk {

using namespace AsterUtils;

extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {

DECLARE_CCTK_ARGUMENTSX_FishboneMoncrief_ET_GRHD_initial;
DECLARE_CCTK_PARAMETERS;

// Set type of atmosphere
atmosphere_t atm_t;
if (CCTK_EQUALS(atmo_type, "isentropic-graded")) {
atm_t = atmosphere_t::isentropic_graded;
} else if (CCTK_EQUALS(atmo_type, "free-graded")) {
atm_t = atmosphere_t::free_graded;
} else if (CCTK_EQUALS(atmo_type, "constant")) {
atm_t = atmosphere_t::constant;
} else {
CCTK_ERROR("Unknown value for parameter \"atmo_type\"");
}

CCTK_VINFO("Fishbone-Moncrief Disk Initial data");
CCTK_VINFO("Using input parameters of\n a = %e,\n M = %e,\nr_in = "
"%e,\nr_at_max_density = %e\nkappa = %e\ngamma = %e",
Expand All @@ -27,7 +44,7 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {

// First compute maximum pressure and density
const CCTK_REAL hm1_init =
FMdisk::GRHD_hm1(xcoord_init, ycoord_init, zcoord_init);
GRHD_hm1(xcoord_init, ycoord_init, zcoord_init, M, r_in, a, r_at_max_density);
const CCTK_REAL rho_max =
pow(hm1_init * (gamma - 1.0) / (kappa * gamma), 1.0 / (gamma - 1.0));
const CCTK_REAL P_max = kappa * pow(rho_max, gamma);
Expand Down Expand Up @@ -79,10 +96,11 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {
CCTK_REAL kDD12_L{ 0. };
CCTK_REAL kDD22_L{ 0. };

FMdisk::KerrSchild(xcoord, ycoord, zcoord, alp_L, betaU0_L, betaU1_L,
KerrSchild(xcoord, ycoord, zcoord, alp_L, betaU0_L, betaU1_L,
betaU2_L, gDD00_L, gDD01_L, gDD02_L, gDD11_L,
gDD12_L, gDD22_L, kDD00_L, kDD01_L, kDD02_L, kDD11_L,
kDD12_L, kDD22_L);
kDD12_L, kDD22_L,
M, a);

alp(p.I) = alp_L;
betax(p.I) = betaU0_L;
Expand Down Expand Up @@ -126,7 +144,7 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {

if (rr > r_in) {

CCTK_REAL hm1 = FMdisk::GRHD_hm1(xcoord, ycoord, zcoord);
CCTK_REAL hm1 = GRHD_hm1(xcoord, ycoord, zcoord, M, r_in, a, r_at_max_density);

if (hm1 > 0) {

Expand All @@ -141,8 +159,9 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {
CCTK_REAL velU1_L{ 0. };
CCTK_REAL velU2_L{ 0. };

FMdisk::GRHD_velocities(xcoord, ycoord, zcoord, velU0_L, velU1_L,
velU2_L);
GRHD_velocities(xcoord, ycoord, zcoord, velU0_L, velU1_L,
velU2_L,
M, a, r_at_max_density);

velx(p.I) = velU0_L;
vely(p.I) = velU1_L;
Expand All @@ -160,10 +179,11 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {

// Outside the disk? Set to atmosphere all hydrodynamic variables!
if (set_to_atmosphere) {

switch (atm_t) {
case atmosphere_t::isentropic_graded : {

if (CCTK_EQUALS(atmo_type, "isentropic-graded")) {

// Choose an atmosphere such that
// Choose an atmosphere such that
// rho = 1e-5 * r^(-3/2), and
// P = k rho^gamma
// Add 1e-100 or 1e-300 to rr or rho to avoid divisions by zero.
Expand All @@ -173,28 +193,32 @@ extern "C" void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {
velx(p.I) = 0.0;
vely(p.I) = 0.0;
velz(p.I) = 0.0;

} else if (CCTK_EQUALS(atmo_type, "free-graded")) {
break;

};
case atmosphere_t::free_graded : {

rho(p.I) = rho_min * pow(rr + 1e-100, -nrho);
press(p.I) = press_min * pow(rr + 1e-100, -npress);
eps(p.I) = press(p.I) / ((rho(p.I) + 1e-300) * (gamma - 1.0));
velx(p.I) = 0.0;
vely(p.I) = 0.0;
velz(p.I) = 0.0;
break;

} else if (CCTK_EQUALS(atmo_type, "constant")) {
};
case atmosphere_t::constant : {

rho(p.I) = rho_min;
press(p.I) = press_min;
eps(p.I) = press_min / ((rho_min + 1e-300) * (gamma - 1.0));
velx(p.I) = 0.0;
vely(p.I) = 0.0;
velz(p.I) = 0.0;
break;

} else {
CCTK_ERROR("Unknown value for parameter \"atmo_type\"");
}
};
}
}
});
}
Expand All @@ -205,9 +229,10 @@ FishboneMoncrief_ET_GRHD_initial__perturb_pressure(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTSX_FishboneMoncrief_ET_GRHD_initial__perturb_pressure;
DECLARE_CCTK_PARAMETERS;

grid.loop_all_device<1, 1, 1>(
// rand() is a host function, thus we cannot run the following code on device
grid.loop_all<1, 1, 1>(
grid.nghostzones,
[=] CCTK_DEVICE(const Loop::PointDesc & p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
[=] CCTK_HOST(const Loop::PointDesc & p) CCTK_ATTRIBUTE_ALWAYS_INLINE {

CCTK_REAL xcoord = p.x;
CCTK_REAL ycoord = p.y;
Expand All @@ -219,15 +244,24 @@ FishboneMoncrief_ET_GRHD_initial__perturb_pressure(CCTK_ARGUMENTS) {

if (rr > r_in) {

hm1 = FMdisk::GRHD_hm1(xcoord, ycoord, zcoord);
hm1 = GRHD_hm1(xcoord, ycoord, zcoord, M, r_in, a, r_at_max_density);

if (hm1 > 0) {

CCTK_REAL press_L = press(p.I);
CCTK_REAL eps_L = eps(p.I);
CCTK_REAL rho_L = rho(p.I);

FMdisk::GRHD_perturb_pressure(press_L, eps_L, rho_L);
// Generate random number in range [0,1),
// snippet courtesy http://daviddeley.com/random/crandom.htm
const CCTK_REAL random_number_between_0_and_1 =
((CCTK_REAL)rand() / ((CCTK_REAL)(RAND_MAX) + (CCTK_REAL)(1)));

const CCTK_REAL random_number_between_min_and_max =
random_min + (random_max - random_min) * random_number_between_0_and_1;

GRHD_perturb_pressure(press_L, eps_L, rho_L,
random_number_between_min_and_max, gamma);

press(p.I) = press_L;
eps(p.I) = eps_L;
Expand Down Expand Up @@ -259,15 +293,16 @@ extern "C" void FishboneMoncrief_Set_A(CCTK_ARGUMENTS) {
CCTK_REAL ytilde =
wrt_rho_max ? ycoord - r_at_max_density * sinphi : ycoord;

CCTK_REAL pressL_stag = FM_Utils::calc_avg_c2e(press, p, 0);
CCTK_REAL rhoL_stag = FM_Utils::calc_avg_c2e(rho, p, 0);
CCTK_REAL pressL_stag = calc_avg_c2e(press, p, 0);
CCTK_REAL rhoL_stag = calc_avg_c2e(rho, p, 0);

CCTK_REAL AxL = 0.;
CCTK_REAL AyL = 0.;
CCTK_REAL AzL = 0.;

FMdisk::GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL);
GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL,
use_pressure, A_b, A_n, A_c, press_cut, rho_cut);

Avec_x(p.I) = AxL;
});
Expand All @@ -289,15 +324,16 @@ extern "C" void FishboneMoncrief_Set_A(CCTK_ARGUMENTS) {
CCTK_REAL ytilde =
wrt_rho_max ? ycoord - r_at_max_density * sinphi : ycoord;

CCTK_REAL pressL_stag = FM_Utils::calc_avg_c2e(press, p, 1);
CCTK_REAL rhoL_stag = FM_Utils::calc_avg_c2e(rho, p, 0);
CCTK_REAL pressL_stag = calc_avg_c2e(press, p, 1);
CCTK_REAL rhoL_stag = calc_avg_c2e(rho, p, 0);

CCTK_REAL AxL = 0.;
CCTK_REAL AyL = 0.;
CCTK_REAL AzL = 0.;

FMdisk::GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL);
GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL,
use_pressure, A_b, A_n, A_c, press_cut, rho_cut);

Avec_y(p.I) = AyL;
});
Expand All @@ -319,15 +355,16 @@ extern "C" void FishboneMoncrief_Set_A(CCTK_ARGUMENTS) {
CCTK_REAL ytilde =
wrt_rho_max ? ycoord - r_at_max_density * sinphi : ycoord;

CCTK_REAL pressL_stag = FM_Utils::calc_avg_c2e(press, p, 2);
CCTK_REAL rhoL_stag = FM_Utils::calc_avg_c2e(rho, p, 0);
CCTK_REAL pressL_stag = calc_avg_c2e(press, p, 2);
CCTK_REAL rhoL_stag = calc_avg_c2e(rho, p, 0);

CCTK_REAL AxL = 0.;
CCTK_REAL AyL = 0.;
CCTK_REAL AzL = 0.;

FMdisk::GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL);
GRMHD_set_A(pressL_stag, rhoL_stag, xtilde, ytilde, AxL, AyL,
AzL,
use_pressure, A_b, A_n, A_c, press_cut, rho_cut);

Avec_z(p.I) = AzL;
});
Expand Down Expand Up @@ -364,10 +401,11 @@ extern "C" void FishboneMoncrief_Set_Spacetime(CCTK_ARGUMENTS) {
CCTK_REAL kDD12_L{ 0. };
CCTK_REAL kDD22_L{ 0. };

FMdisk::KerrSchild(xcoord, ycoord, zcoord, alp_L, betaU0_L, betaU1_L,
KerrSchild(xcoord, ycoord, zcoord, alp_L, betaU0_L, betaU1_L,
betaU2_L, gDD00_L, gDD01_L, gDD02_L, gDD11_L,
gDD12_L, gDD22_L, kDD00_L, kDD01_L, kDD02_L, kDD11_L,
kDD12_L, kDD22_L);
kDD12_L, kDD22_L,
M, a);

alp(p.I) = alp_L;
betax(p.I) = betaU0_L;
Expand Down Expand Up @@ -396,3 +434,5 @@ extern "C" void FishboneMoncrief_Set_Spacetime(CCTK_ARGUMENTS) {
kyz(p.I) = kDD12_L;
});
}

} // namespace FMdisk
53 changes: 18 additions & 35 deletions FishboneMoncriefIDX/src/FM_disk_implementation.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,15 @@
#define FM_IMPL_HXX

#include <cmath>
#include <cstdlib> // Needed for rand()
#include <cctk_Parameters.h>

namespace FMdisk {

CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
GRHD_hm1(CCTK_REAL const &xcoord, CCTK_REAL const &ycoord,
CCTK_REAL const &zcoord) {

DECLARE_CCTK_PARAMETERS;
CCTK_REAL const &zcoord,
CCTK_REAL const M, CCTK_REAL const r_in, CCTK_REAL const a,
CCTK_REAL const r_at_max_density) {

const CCTK_REAL tmp_2 = 2 * M * r_in;
const CCTK_REAL tmp_3 = ((a) * (a));
Expand Down Expand Up @@ -63,9 +62,8 @@ KerrSchild(CCTK_REAL const &xcoord, CCTK_REAL const &ycoord,
CCTK_REAL &gammaDD01, CCTK_REAL &gammaDD02, CCTK_REAL &gammaDD11,
CCTK_REAL &gammaDD12, CCTK_REAL &gammaDD22, CCTK_REAL &KDD00,
CCTK_REAL &KDD01, CCTK_REAL &KDD02, CCTK_REAL &KDD11,
CCTK_REAL &KDD12, CCTK_REAL &KDD22) {

DECLARE_CCTK_PARAMETERS;
CCTK_REAL &KDD12, CCTK_REAL &KDD22,
CCTK_REAL const M, CCTK_REAL const a) {

const CCTK_REAL FDPart3_0 = ((a) * (a));
const CCTK_REAL FDPart3_1 = ((zcoord) * (zcoord));
Expand Down Expand Up @@ -238,9 +236,9 @@ CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
GRHD_velocities(CCTK_REAL const &xcoord, CCTK_REAL const &ycoord,
CCTK_REAL const &zcoord, CCTK_REAL &Valencia3velocityU0GF,
CCTK_REAL &Valencia3velocityU1GF,
CCTK_REAL &Valencia3velocityU2GF) {

DECLARE_CCTK_PARAMETERS;
CCTK_REAL &Valencia3velocityU2GF,
CCTK_REAL const M, CCTK_REAL const a,
CCTK_REAL const r_at_max_density) {

const CCTK_REAL FDPart3_0 = ((a) * (a));
const CCTK_REAL FDPart3_2 =
Expand Down Expand Up @@ -304,41 +302,26 @@ GRHD_velocities(CCTK_REAL const &xcoord, CCTK_REAL const &ycoord,
}

CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
GRHD_perturb_pressure(CCTK_REAL &press, CCTK_REAL &eps, CCTK_REAL const &rho) {

DECLARE_CCTK_PARAMETERS;

// Generate random number in range [0,1),
// snippet courtesy http://daviddeley.com/random/crandom.htm
const CCTK_REAL random_number_between_0_and_1 =
((CCTK_REAL)rand() / ((CCTK_REAL)(RAND_MAX) + (CCTK_REAL)(1)));

const CCTK_REAL random_number_between_min_and_max =
random_min + (random_max - random_min) * random_number_between_0_and_1;
GRHD_perturb_pressure(CCTK_REAL &press, CCTK_REAL &eps, CCTK_REAL const &rho,
CCTK_REAL const random_number_between_min_and_max,
CCTK_REAL const gamma) {

press = press * (1.0 + random_number_between_min_and_max);

// Add 1e-300 to rho to avoid division by zero when density is zero.
eps = press / ((rho + 1e-300) * (gamma - 1.0));
}

// CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
// GRMHD_set_A(CCTK_REAL const &press, CCTK_REAL const &xtilde, CCTK_REAL const
// &ytilde, CCTK_REAL &Ax, CCTK_REAL &Ay, CCTK_REAL &Az) {
//
// DECLARE_CCTK_PARAMETERS;
//
// Ax = - A_b * pow(fmax(press-press_cut,0.),A_n) * ytilde;
// Ay = A_b * pow(fmax(press-press_cut,0.),A_n) * xtilde;
// Az = 0.;
//}

CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
GRMHD_set_A(CCTK_REAL const &press, CCTK_REAL const &rho,
CCTK_REAL const &xtilde, CCTK_REAL const &ytilde, CCTK_REAL &Ax,
CCTK_REAL &Ay, CCTK_REAL &Az) {

DECLARE_CCTK_PARAMETERS;
CCTK_REAL &Ay, CCTK_REAL &Az,
bool const use_pressure,
CCTK_REAL const A_b,
CCTK_REAL const A_n,
CCTK_REAL const A_c,
CCTK_REAL const press_cut,
CCTK_REAL const rho_cut) {

const CCTK_REAL rcyl = fmax(sqrt(xtilde * xtilde + ytilde * ytilde), 1e-15);

Expand Down
Loading

0 comments on commit ec368fc

Please sign in to comment.