Skip to content

Commit

Permalink
Merge pull request cholla-hydro#356 from mabruzzo/dual-energy-fix
Browse files Browse the repository at this point in the history
fixed dual-energy formalism synchronization.
  • Loading branch information
evaneschneider authored Jan 3, 2024
2 parents b30468a + 0c6efd8 commit ec7e7c3
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 7 deletions.
23 changes: 19 additions & 4 deletions src/global/global.h
Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,25 @@ typedef double Real;
#define TEMP_FLOOR 1e-3
#define DENS_FLOOR 1e-5 // in code units

// Parameter for Enzo dual Energy Condition
#define DE_ETA_1 \
0.001 // Ratio of U to E for which Internal Energy is used to compute the
// Pressure
// Parameters for Enzo dual Energy Condition
// - Prior to GH PR #356, DE_ETA_1 nominally had a value of 0.001 in all
// simulations (in practice, the value of DE_ETA_1 had minimal significance
// in those simulations). In PR #356, we revised the internal-energy
// synchronization to account for the value of DE_ETA_1. This was necessary
// for non-cosmology simulations.
// - In Cosmological simulation, we set DE_ETA_1 to a large number (it doesn't
// really matter what, as long as its >=1) to maintain the older behavior
// - In the future, we run tests and revisit the choice of DE_ETA_1 in
// cosmological simulations
#ifdef COSMOLOGY
#define DE_ETA_1 10.0
#else
#define DE_ETA_1 \
0.001 // Ratio of U to E for which Internal Energy is used to compute the
// Pressure. This also affects when the Internal Energy is used for
// the update.
#endif

#define DE_ETA_2 \
0.035 // Ratio of U to max(E_local) used to select which Internal Energy is
// used for the update.
Expand Down
18 changes: 15 additions & 3 deletions src/hydro/hydro_cuda.cu
Original file line number Diff line number Diff line change
Expand Up @@ -840,6 +840,7 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho
int imo, ipo;
n_cells = nx;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand All @@ -865,7 +866,10 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho
Emax = fmax(dev_conserved[4 * n_cells + imo], E);
Emax = fmax(Emax, dev_conserved[4 * n_cells + ipo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand All @@ -888,6 +892,7 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i
int imo, ipo, jmo, jpo;
n_cells = nx * ny;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand Down Expand Up @@ -923,7 +928,10 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i
Emax = fmax(Emax, dev_conserved[4 * n_cells + jmo]);
Emax = fmax(Emax, dev_conserved[4 * n_cells + jpo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand All @@ -946,6 +954,7 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i
int imo, ipo, jmo, jpo, kmo, kpo;
n_cells = nx * ny * nz;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand Down Expand Up @@ -988,7 +997,10 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i
Emax = fmax(Emax, dev_conserved[4 * n_cells + kmo]);
Emax = fmax(Emax, dev_conserved[4 * n_cells + kpo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand Down

0 comments on commit ec7e7c3

Please sign in to comment.