From 1655a056ec46b2095159052358b521653004a124 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 17 May 2024 15:50:20 -0400 Subject: [PATCH] fix openmp race conditions in diag (#455) there are a few regions where it appears we need a critical block --- Source/MaestroDiag.cpp | 136 +++++++++++++---------------------------- 1 file changed, 43 insertions(+), 93 deletions(-) diff --git a/Source/MaestroDiag.cpp b/Source/MaestroDiag.cpp index 28d648b0a..fa9d40b1c 100644 --- a/Source/MaestroDiag.cpp +++ b/Source/MaestroDiag.cpp @@ -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; @@ -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) @@ -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 } } @@ -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); } } } @@ -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) {