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

Jump to bottom

fixed dual-energy formalism synchronization. #356

Merged
merged 3 commits into from
Jan 3, 2024
Merged

fixed dual-energy formalism synchronization. #356

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
9bb2328
fixed dual-energy formalism synchronization.
mabruzzo Dec 7, 2023
fc8f8d7
revert de formalism behavior in cosmological simulations until future…
mabruzzo Dec 20, 2023
0c6efd8
Merge branch 'dev' into dual-energy-fix
mabruzzo Dec 20, 2023
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
23 changes: 19 additions & 4 deletions src/global/global.h
View file
Open in desktop
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
View file
Open in desktop
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