Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fishbone moncrief #80

Merged
merged 3 commits into from
Dec 12, 2024
Merged

Fishbone moncrief #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9ffbb70
FishboneMoncriefIDX: Couple to AsterUtils, compile on GPU
Dec 10, 2024
0c721ac
FishboneMoncriefIDX: Remove comments
Dec 10, 2024
caa627f
FishboneMoncriefIDX: Add compilation of FM disk thorn
Dec 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion FishboneMoncriefIDX/configuration.ccl
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading