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.

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

Making static gravity more user friendly with custom gravity field flags #340

Merged
merged 16 commits into from
Oct 27, 2023
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
clang format
  • Loading branch information
evazlimen committed Sep 21, 2023
commit 11ae4813c2b98072a81136b211dc40f4b916e6a1
4 changes: 2 additions & 2 deletions src/global/global.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ void parse_params(char *param_file, struct parameters *parms, int argc, char **a
#endif /*ROTATED_PROJECTION*/

#ifdef STATIC_GRAV
//initialize custom gravity flag to zero
// initialize custom gravity flag to zero
parms->custom_grav = 0;
bcaddy marked this conversation as resolved.
Show resolved Hide resolved
#endif

Expand Down Expand Up @@ -219,7 +219,7 @@ void parse_param(char *name, char *value, struct parameters *parms)
#ifdef STATIC_GRAV
} else if (strcmp(name, "custom_grav") == 0) {
parms->custom_grav = atoi(value);
#endif
#endif
} else if (strcmp(name, "tout") == 0) {
parms->tout = atof(value);
} else if (strcmp(name, "outstep") == 0) {
Expand Down
2 changes: 1 addition & 1 deletion src/global/global.h
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ struct parameters {
int out_float32_GasEnergy = 0;
#endif
#ifdef STATIC_GRAV
int custom_grav; //flag to set specific static gravity field
int custom_grav; // flag to set specific static gravity field
bcaddy marked this conversation as resolved.
Show resolved Hide resolved
#endif
#ifdef MHD
int out_float32_magnetic_x = 0;
Expand Down
217 changes: 108 additions & 109 deletions src/gravity/static_grav.h
Original file line number Diff line number Diff line change
Expand Up @@ -18,87 +18,85 @@ inline __device__ void calc_g_1D(int xid, int x_off, int n_ghost, int custom_gra
{
Real x_pos, r_disk, r_halo;
x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
switch(custom_grav){
case 1:
// for disk components, calculate polar r
// r_disk = 0.220970869121;
// r_disk = 6.85009694274;
r_disk = 13.9211647546;
// r_disk = 20.9922325665;
// for halo, calculate spherical r
r_halo = sqrt(x_pos * x_pos + r_disk * r_disk);

// set properties of halo and disk (these must match initial conditions)
Real a_disk_z, a_halo, M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
M_vir = 1.0e12; // viral mass of MW in M_sun
M_d = 6.5e10; // mass of disk in M_sun
M_h = M_vir - M_d; // halo mass in M_sun
R_vir = 261; // viral radius in kpc
c_vir = 20.0; // halo concentration
R_h = R_vir / c_vir; // halo scale length in kpc
R_d = 3.5; // disk scale length in kpc
z_d = 3.5 / 5.0; // disk scale height in kpc
phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
x = r_halo / R_h;

// calculate acceleration due to NFW halo & Miyamoto-Nagai disk
a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
a_disk_z =
-GN * M_d * x_pos * (R_d + sqrt(x_pos * x_pos + z_d * z_d)) /
(pow(r_disk * r_disk + pow2(R_d + sqrt(x_pos * x_pos + z_d * z_d)), 1.5) * sqrt(x_pos * x_pos + z_d * z_d));

// total acceleration is the sum of the halo + disk components
*gx = (x_pos / r_halo) * a_halo + a_disk_z;
break;
default:
*gx = 0;
switch (custom_grav) {
case 1:
// for disk components, calculate polar r
// r_disk = 0.220970869121;
// r_disk = 6.85009694274;
r_disk = 13.9211647546;
// r_disk = 20.9922325665;
// for halo, calculate spherical r
r_halo = sqrt(x_pos * x_pos + r_disk * r_disk);

// set properties of halo and disk (these must match initial conditions)
Real a_disk_z, a_halo, M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
M_vir = 1.0e12; // viral mass of MW in M_sun
M_d = 6.5e10; // mass of disk in M_sun
M_h = M_vir - M_d; // halo mass in M_sun
R_vir = 261; // viral radius in kpc
c_vir = 20.0; // halo concentration
R_h = R_vir / c_vir; // halo scale length in kpc
R_d = 3.5; // disk scale length in kpc
z_d = 3.5 / 5.0; // disk scale height in kpc
phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
x = r_halo / R_h;

// calculate acceleration due to NFW halo & Miyamoto-Nagai disk
a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
a_disk_z =
-GN * M_d * x_pos * (R_d + sqrt(x_pos * x_pos + z_d * z_d)) /
(pow(r_disk * r_disk + pow2(R_d + sqrt(x_pos * x_pos + z_d * z_d)), 1.5) * sqrt(x_pos * x_pos + z_d * z_d));

// total acceleration is the sum of the halo + disk components
*gx = (x_pos / r_halo) * a_halo + a_disk_z;
break;
default:
*gx = 0;
}
return;
}

inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_ghost, int custom_grav, Real dx, Real dy, Real xbound,
Real ybound, Real *gx, Real *gy)
inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_ghost, int custom_grav, Real dx, Real dy,
Real xbound, Real ybound, Real *gx, Real *gy)
{
Real x_pos, y_pos, r, phi;
// use the subgrid offset and global boundaries to calculate absolute
// positions on the grid
x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
y_pos = (y_off + yid - n_ghost + 0.5) * dy + ybound;
// for Gresho, also need r & phi
r = sqrt(x_pos * x_pos + y_pos * y_pos);
phi = atan2(y_pos, x_pos);
switch(custom_grav){
// for Gresho, also need r & phi
r = sqrt(x_pos * x_pos + y_pos * y_pos);
phi = atan2(y_pos, x_pos);
switch (custom_grav) {
case 1:
// printf("gresho\n");
// set acceleration to balance v_phi in Gresho problem
if (r < 0.2) {
*gx = -cos(phi)*25.0*r;
*gy = -sin(phi)*25.0*r;
}
else if (r >= 0.2 && r < 0.4) {
*gx = -cos(phi)*(4.0 - 20.0*r + 25.0*r*r)/r;
*gy = -sin(phi)*(4.0 - 20.0*r + 25.0*r*r)/r;
}
else {
*gx = 0.0;
*gy = 0.0;
}
break;
// set acceleration to balance v_phi in Gresho problem
if (r < 0.2) {
*gx = -cos(phi) * 25.0 * r;
*gy = -sin(phi) * 25.0 * r;
} else if (r >= 0.2 && r < 0.4) {
*gx = -cos(phi) * (4.0 - 20.0 * r + 25.0 * r * r) / r;
*gy = -sin(phi) * (4.0 - 20.0 * r + 25.0 * r * r) / r;
} else {
*gx = 0.0;
*gy = 0.0;
}
break;
case 2:
//printf("rayleigh talor\n");
// printf("rayleigh talor\n");
*gx = 0;
*gy = -1;
break;
case 3:
//printf("keplerian\n");
// printf("keplerian\n");
Real M;
M = 1*MSUN_CGS;
*gx = -cos(phi)*GN*M/(r*r);
*gy = -sin(phi)*GN*M/(r*r);
break;
M = 1 * MSUN_CGS;
*gx = -cos(phi) * GN * M / (r * r);
*gy = -sin(phi) * GN * M / (r * r);
break;
case 4:
//printf("disk\n");
// set gravitational acceleration for Kuzmin disk + NFW halo
// printf("disk\n");
// set gravitational acceleration for Kuzmin disk + NFW halo
Real a_d, a_h, a, M_vir, M_d, R_vir, R_d, R_s, M_h, c_vir, x;
M_vir = 1.0e12; // viral mass of MW in M_sun
M_d = 6.5e10; // mass of disk in M_sun (assume all gas)
Expand All @@ -118,16 +116,17 @@ inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_g
*gy = -sin(phi) * a;
break;
default:
//printf("default\n");
// printf("default\n");
*gx = 0;
*gy = 0;
}
}

return;
}

inline __device__ void calc_g_3D(int xid, int yid, int zid, int x_off, int y_off, int z_off, int n_ghost, int custom_grav, Real dx,
Real dy, Real dz, Real xbound, Real ybound, Real zbound, Real *gx, Real *gy, Real *gz)
inline __device__ void calc_g_3D(int xid, int yid, int zid, int x_off, int y_off, int z_off, int n_ghost,
int custom_grav, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
Real *gx, Real *gy, Real *gz)
{
Real x_pos, y_pos, z_pos, r_disk, r_halo;
// use the subgrid offset and global boundaries to calculate absolute
Expand All @@ -140,49 +139,49 @@ inline __device__ void calc_g_3D(int xid, int yid, int zid, int x_off, int y_off
r_disk = sqrt(x_pos * x_pos + y_pos * y_pos);
// for halo, calculate spherical r
r_halo = sqrt(x_pos * x_pos + y_pos * y_pos + z_pos * z_pos);
switch(custom_grav){
case 1:
// set properties of halo and disk (these must match initial conditions)
Real a_disk_r, a_disk_z, a_halo, a_halo_r, a_halo_z;
Real M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
// MW model
M_vir = 1.0e12; // viral mass of in M_sun
M_d = 6.5e10; // viral mass of in M_sun
R_d = 3.5; // disk scale length in kpc
z_d = 3.5 / 5.0; // disk scale height in kpc
R_vir = 261.; // virial radius in kpc
c_vir = 20.0; // halo concentration
// M82 model
// M_vir = 5.0e10; // viral mass of in M_sun
// M_d = 1.0e10; // mass of disk in M_sun
// R_d = 0.8; // disk scale length in kpc
// z_d = 0.15; // disk scale height in kpc
// R_vir = R_d/0.015; // viral radius in kpc
// c_vir = 10.0; // halo concentration

M_h = M_vir - M_d; // halo mass in M_sun
R_h = R_vir / c_vir; // halo scale length in kpc
phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
x = r_halo / R_h;

// calculate acceleration due to NFW halo & Miyamoto-Nagai disk
a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
a_halo_r = a_halo * (r_disk / r_halo);
a_halo_z = a_halo * (z_pos / r_halo);
a_disk_r = -GN * M_d * r_disk * pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), -1.5);
a_disk_z =
-GN * M_d * z_pos * (R_d + sqrt(z_pos * z_pos + z_d * z_d)) /
(pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), 1.5) * sqrt(z_pos * z_pos + z_d * z_d));

// total acceleration is the sum of the halo + disk components
*gx = (x_pos / r_disk) * (a_disk_r + a_halo_r);
*gy = (y_pos / r_disk) * (a_disk_r + a_halo_r);
*gz = a_disk_z + a_halo_z;
break;
default:
*gx = 0;
*gy = 0;
*gz = 0;
switch (custom_grav) {
case 1:
// set properties of halo and disk (these must match initial conditions)
Real a_disk_r, a_disk_z, a_halo, a_halo_r, a_halo_z;
Real M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
// MW model
M_vir = 1.0e12; // viral mass of in M_sun
M_d = 6.5e10; // viral mass of in M_sun
R_d = 3.5; // disk scale length in kpc
z_d = 3.5 / 5.0; // disk scale height in kpc
R_vir = 261.; // virial radius in kpc
c_vir = 20.0; // halo concentration
// M82 model
// M_vir = 5.0e10; // viral mass of in M_sun
// M_d = 1.0e10; // mass of disk in M_sun
// R_d = 0.8; // disk scale length in kpc
// z_d = 0.15; // disk scale height in kpc
// R_vir = R_d/0.015; // viral radius in kpc
// c_vir = 10.0; // halo concentration

M_h = M_vir - M_d; // halo mass in M_sun
R_h = R_vir / c_vir; // halo scale length in kpc
phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
x = r_halo / R_h;

// calculate acceleration due to NFW halo & Miyamoto-Nagai disk
a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
a_halo_r = a_halo * (r_disk / r_halo);
a_halo_z = a_halo * (z_pos / r_halo);
a_disk_r = -GN * M_d * r_disk * pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), -1.5);
a_disk_z =
-GN * M_d * z_pos * (R_d + sqrt(z_pos * z_pos + z_d * z_d)) /
(pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), 1.5) * sqrt(z_pos * z_pos + z_d * z_d));

// total acceleration is the sum of the halo + disk components
*gx = (x_pos / r_disk) * (a_disk_r + a_halo_r);
*gy = (y_pos / r_disk) * (a_disk_r + a_halo_r);
*gz = a_disk_z + a_halo_z;
break;
default:
*gx = 0;
*gy = 0;
*gz = 0;
}
return;
}
Expand Down
10 changes: 5 additions & 5 deletions src/grid/grid3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -148,11 +148,11 @@ void Grid3D::Initialize(struct parameters *P)
int nz_in = P->nz;

#ifdef STATIC_GRAV
H.custom_grav = P->custom_grav; //Initialize the custom static gravity flag
H.custom_grav = P->custom_grav; // Initialize the custom static gravity flag
printf("H.custom_grav is %d\n", H.custom_grav);
if (H.custom_grav == 0){
if (H.custom_grav == 0) {
printf("WARNING: No custom gravity field given. Gravity field will be set to zero.\n");
}
}
#endif

// Set the CFL coefficient (a global variable)
Expand Down Expand Up @@ -467,8 +467,8 @@ Real Grid3D::Update_Grid(void)
#endif // VL
#ifdef SIMPLE
Simple_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, density_floor, U_floor,
C.Grav_potential);
H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, density_floor,
U_floor, C.Grav_potential);
#endif // SIMPLE
#endif
} else {
Expand Down
6 changes: 4 additions & 2 deletions src/hydro/hydro_cuda.cu
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,8 @@ __global__ void Update_Conserved_Variables_3D(Real *dev_conserved, Real *Q_Lx, R
Real *Q_Lz, Real *Q_Rz, Real *dev_F_x, Real *dev_F_y, Real *dev_F_z,
int nx, int ny, int nz, int x_off, int y_off, int z_off, int n_ghost,
Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound, Real dt,
Real gamma, int n_fields, int custom_grav, Real density_floor, Real *dev_potential)
Real gamma, int n_fields, int custom_grav, Real density_floor,
Real *dev_potential)
{
int id, xid, yid, zid, n_cells;
int imo, jmo, kmo;
Expand Down Expand Up @@ -299,7 +300,8 @@ __global__ void Update_Conserved_Variables_3D(Real *dev_conserved, Real *Q_Lx, R
#endif // DENSITY_FLOOR

#ifdef STATIC_GRAV
calc_g_3D(xid, yid, zid, x_off, y_off, z_off, n_ghost, custom_grav, dx, dy, dz, xbound, ybound, zbound, &gx, &gy, &gz);
calc_g_3D(xid, yid, zid, x_off, y_off, z_off, n_ghost, custom_grav, dx, dy, dz, xbound, ybound, zbound, &gx, &gy,
&gz);
d_n = dev_conserved[id];
d_inv_n = 1.0 / d_n;
vx_n = dev_conserved[1 * n_cells + id] * d_inv_n;
Expand Down
3 changes: 2 additions & 1 deletion src/hydro/hydro_cuda.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@ __global__ void Update_Conserved_Variables_3D(Real *dev_conserved, Real *Q_Lx, R
Real *Q_Lz, Real *Q_Rz, Real *dev_F_x, Real *dev_F_y, Real *dev_F_z,
int nx, int ny, int nz, int x_off, int y_off, int z_off, int n_ghost,
Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound, Real dt,
Real gamma, int n_fields, int custom_grav, Real density_floor, Real *dev_potential);
Real gamma, int n_fields, int custom_grav, Real density_floor,
Real *dev_potential);

/*!
* \brief Determine the maximum inverse crossing time in a specific cell
Expand Down
3 changes: 2 additions & 1 deletion src/integrators/simple_3D_cuda.cu
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@

void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor, Real *host_grav_potential)
Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
Real *host_grav_potential)
{
// Here, *dev_conserved contains the entire
// set of conserved variables on the grid
Expand Down
3 changes: 2 additions & 1 deletion src/integrators/simple_3D_cuda.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@

void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor, Real *host_grav_potential);
Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
Real *host_grav_potential);

void Free_Memory_Simple_3D();

Expand Down
2 changes: 1 addition & 1 deletion src/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ int main(int argc, char *argv[])
Write_Message_To_Log_File(message.c_str());
message = "Macro Flags = " + std::string(MACRO_FLAGS);
Write_Message_To_Log_File(message.c_str());

// initialize the grid
G.Initialize(&P);
chprintf("Local number of grid cells: %d %d %d %d\n", G.H.nx_real, G.H.ny_real, G.H.nz_real, G.H.n_cells);
Expand Down
Loading