Skip to content

Commit

Permalink
fix openmp race conditions in diag (#455)
Browse files Browse the repository at this point in the history 
there are a few regions where it appears we need a critical block
  • Loading branch information
@zingale
zingale authored May 17, 2024
1 parent e8db3be commit 1655a05
Showing 1 changed file with 43 additions and 93 deletions.
136 changes: 43 additions & 93 deletions Source/MaestroDiag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,122 +178,84 @@ void Maestro::DiagFile(const int step, const Real t_in,

// For spherical, we only consider cells inside of where the
// sponging begins
if (cell_valid &&
(!spherical ||
scal(i, j, k, Rho) >=
sponge_start_factor * sponge_center_density)) {
if (cell_valid && (!spherical || scal(i, j, k, Rho) >=
sponge_start_factor * sponge_center_density)) {
Real vel = 0.0;
if (spherical) {
#if (AMREX_SPACEDIM == 3)
// is it one of the 8 zones surrounding the center?
if (amrex::Math::abs(x - center[0]) < dx[0] &&
amrex::Math::abs(y - center[1]) < dx[1] &&
amrex::Math::abs(z - center[2]) < dx[2]) {
#pragma omp critical
if (std::abs(x - center[0]) < dx[0] &&
std::abs(y - center[1]) < dx[1] &&
std::abs(z - center[2]) < dx[2]) {
ncenter_level++;

T_center_level += scal(i, j, k, Temp);

vel_center_level[0] +=
u(i, j, k, 0) +
0.5 * (w0macx(i, j, k) +
w0macx(i + 1, j, k));
u(i, j, k, 0) + 0.5 * (w0macx(i, j, k) + w0macx(i + 1, j, k));
vel_center_level[1] +=
u(i, j, k, 1) +
0.5 * (w0macy(i, j, k) +
w0macy(i, j + 1, k));
u(i, j, k, 1) + 0.5 * (w0macy(i, j, k) + w0macy(i, j + 1, k));
vel_center_level[2] +=
u(i, j, k, 2) +
0.5 * (w0macz(i, j, k) +
w0macz(i, j, k + 1));
u(i, j, k, 2) + 0.5 * (w0macz(i, j, k) + w0macz(i, j, k + 1));
}

// velr is the projection of the velocity (including w0) onto
// the radial unit vector
// Real velr = u(i,j,k,0)*normal_arr(i,j,k,0) + \ //
// u(i,j,k,1)*normal_arr(i,j,k,1) + \ //
// u(i,j,k,2)*normal_arr(i,j,k,2) + w0r(i,j,k);
// u(i,j,k,1)*normal_arr(i,j,k,1) + \ //
// u(i,j,k,2)*normal_arr(i,j,k,2) + w0r(i,j,k);

// vel is the magnitude of the velocity, including w0
vel = std::sqrt(
(u(i, j, k, 0) +
0.5 * (w0macx(i, j, k) +
w0macx(i + 1, j, k))) *
(u(i, j, k, 0) +
0.5 * (w0macx(i, j, k) +
w0macx(i + 1, j, k))) +
(u(i, j, k, 1) +
0.5 * (w0macy(i, j, k) +
w0macy(i, j + 1, k))) *
(u(i, j, k, 1) +
0.5 * (w0macy(i, j, k) +
w0macy(i, j + 1, k))) +
(u(i, j, k, 2) +
0.5 * (w0macz(i, j, k) +
w0macz(i, j, k + 1))) *
(u(i, j, k, 2) +
0.5 * (w0macz(i, j, k) +
w0macz(i, j, k + 1))));
amrex::Math::powi<2>(u(i, j, k, 0) + 0.5 * (w0macx(i, j, k) + w0macx(i + 1, j, k))) +
amrex::Math::powi<2>(u(i, j, k, 1) + 0.5 * (w0macy(i, j, k) + w0macy(i, j + 1, k))) +
amrex::Math::powi<2>(u(i, j, k, 2) + 0.5 * (w0macz(i, j, k) + w0macz(i, j, k + 1))));

// max T, location, and velocity at that location (including w0)
#pragma omp critical
if (scal(i, j, k, Temp) > T_max_local) {
T_max_local = scal(i, j, k, Temp);
coord_Tmax_local[0] = x;
coord_Tmax_local[1] = y;
coord_Tmax_local[2] = z;
vel_Tmax_local[0] =
u(i, j, k, 0) +
0.5 * (w0macx(i, j, k) +
w0macx(i + 1, j, k));
u(i, j, k, 0) + 0.5 * (w0macx(i, j, k) + w0macx(i + 1, j, k));
vel_Tmax_local[1] =
u(i, j, k, 1) +
0.5 * (w0macy(i, j, k) +
w0macy(i, j + 1, k));
u(i, j, k, 1) + 0.5 * (w0macy(i, j, k) + w0macy(i, j + 1, k));
vel_Tmax_local[2] =
u(i, j, k, 2) +
0.5 * (w0macz(i, j, k) +
w0macz(i, j, k + 1));
u(i, j, k, 2) + 0.5 * (w0macz(i, j, k) + w0macz(i, j, k + 1));
}

// max enuc
if (rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho) >
enuc_max_local) {
enuc_max_local = rho_Hnuc_arr(i, j, k) /
scal(i, j, k, Rho);
#pragma omp critical
if (rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho) > enuc_max_local) {
enuc_max_local = rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho);
coord_enucmax_local[0] = x;
coord_enucmax_local[1] = y;
coord_enucmax_local[2] = z;
vel_enucmax_local[0] =
u(i, j, k, 0) +
0.5 * (w0macx(i, j, k) +
w0macx(i + 1, j, k));
u(i, j, k, 0) + 0.5 * (w0macx(i, j, k) + w0macx(i + 1, j, k));
vel_enucmax_local[1] =
u(i, j, k, 1) +
0.5 * (w0macy(i, j, k) +
w0macy(i, j + 1, k));
u(i, j, k, 1) + 0.5 * (w0macy(i, j, k) + w0macy(i, j + 1, k));
vel_enucmax_local[2] =
u(i, j, k, 2) +
0.5 * (w0macz(i, j, k) +
w0macz(i, j, k + 1));
u(i, j, k, 2) + 0.5 * (w0macz(i, j, k) + w0macz(i, j, k + 1));
}
#endif
} else {
// vel is the magnitude of the velocity, including w0
#if (AMREX_SPACEDIM == 2)
Real vert_vel =
u(i, j, k, 1) +
0.5 * (w0_arr(lev, j) + w0_arr(lev, j + 1));
vel = std::sqrt(u(i, j, k, 0) * u(i, j, k, 0) +
vert_vel * vert_vel);
Real vert_vel = u(i, j, k, 1) + 0.5 * (w0_arr(lev, j) + w0_arr(lev, j + 1));
vel = std::sqrt(u(i, j, k, 0) * u(i, j, k, 0) + vert_vel * vert_vel);
#else
Real vert_vel =
u(i, j, k, 2) +
0.5 * (w0_arr(lev, k) + w0_arr(lev, k + 1));
Real vert_vel = u(i, j, k, 2) + 0.5 * (w0_arr(lev, k) + w0_arr(lev, k + 1));
vel = std::sqrt(u(i, j, k, 0) * u(i, j, k, 0) +
u(i, j, k, 1) * u(i, j, k, 1) +
vert_vel * vert_vel);
#endif

// max T, location, and velocity at that location (including w0)
#pragma omp critical
if (scal(i, j, k, Temp) > T_max_local) {
T_max_local = scal(i, j, k, Temp);
coord_Tmax_local[0] = x;
Expand All @@ -305,21 +267,17 @@ void Maestro::DiagFile(const int step, const Real t_in,
vel_Tmax_local[1] = u(i, j, k, 1);
#if (AMREX_SPACEDIM == 2)
vel_Tmax_local[1] +=
0.5 *
(w0_arr(lev, j) + w0_arr(lev, j + 1));
0.5 * (w0_arr(lev, j) + w0_arr(lev, j + 1));
#else
vel_Tmax_local[2] =
u(i, j, k, 2) +
0.5 * (w0_arr(lev, k) +
w0_arr(lev, k + 1));
u(i, j, k, 2) + 0.5 * (w0_arr(lev, k) + w0_arr(lev, k + 1));
#endif
}

// max enuc
if (rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho) >
enuc_max_local) {
enuc_max_local = rho_Hnuc_arr(i, j, k) /
scal(i, j, k, Rho);
#pragma omp critical
if (rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho) > enuc_max_local) {
enuc_max_local = rho_Hnuc_arr(i, j, k) / scal(i, j, k, Rho);
coord_enucmax_local[0] = x;
coord_enucmax_local[1] = y;
#if (AMREX_SPACEDIM == 3)
Expand All @@ -329,13 +287,10 @@ void Maestro::DiagFile(const int step, const Real t_in,
vel_enucmax_local[1] = u(i, j, k, 1);
#if (AMREX_SPACEDIM == 2)
vel_enucmax_local[1] +=
0.5 *
(w0_arr(lev, j) + w0_arr(lev, j + 1));
0.5 * (w0_arr(lev, j) + w0_arr(lev, j + 1));
#else
vel_enucmax_local[2] =
u(i, j, k, 2) +
0.5 * (w0_arr(lev, k) +
w0_arr(lev, k + 1));
u(i, j, k, 2) + 0.5 * (w0_arr(lev, k) + w0_arr(lev, k + 1));
#endif
}
}
Expand All @@ -347,30 +302,25 @@ void Maestro::DiagFile(const int step, const Real t_in,
eos_state.rho = scal(i, j, k, Rho);
for (auto comp = 0; comp < NumSpec; ++comp) {
eos_state.xn[comp] =
scal(i, j, k, FirstSpec + comp) /
eos_state.rho;
scal(i, j, k, FirstSpec + comp) / eos_state.rho;
}
#if NAUX_NET > 0
for (auto comp = 0; comp < NumAux; ++comp) {
eos_state.aux[comp] =
scal(i, j, k, FirstAux + comp) /
eos_state.rho;
scal(i, j, k, FirstAux + comp) / eos_state.rho;
}
#endif

eos(eos_input_rt, eos_state);

// kinetic, internal, and nuclear energies
kin_ener_level +=
weight * scal(i, j, k, Rho) * vel * vel;
int_ener_level +=
weight * scal(i, j, k, Rho) * eos_state.e;
kin_ener_level += weight * scal(i, j, k, Rho) * vel * vel;
int_ener_level += weight * scal(i, j, k, Rho) * eos_state.e;
nuc_ener_level += weight * rho_Hnuc_arr(i, j, k);

// max vel and Mach number
U_max_level = amrex::max(U_max_level, vel);
Mach_max_level =
amrex::max(Mach_max_level, vel / eos_state.cs);
U_max_level = std::max(U_max_level, vel);
Mach_max_level = std::max(Mach_max_level, vel / eos_state.cs);
}
}
}
Expand Down Expand Up @@ -505,8 +455,8 @@ void Maestro::DiagFile(const int step, const Real t_in,
int_ener += int_ener_level;
nuc_ener += nuc_ener_level;

U_max = amrex::max(U_max, U_max_level);
Mach_max = amrex::max(Mach_max, Mach_max_level);
U_max = std::max(U_max, U_max_level);
Mach_max = std::max(Mach_max, Mach_max_level);

// if T_max_level is the new max, then copy the location as well
if (T_max_level > T_max) {
Expand Down

0 comments on commit 1655a05

Please sign in to comment.