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

fix openmp race conditions in diag #455

Merged
merged 4 commits into from
May 17, 2024
Merged

fix openmp race conditions in diag #455

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
976bbc0
fix openmp race conditions in diag
zingale May 15, 2024
9c639a1
Merge branch 'development' into fix_omp_diag
zingale May 16, 2024
e4325f1
Merge branch 'development' into fix_omp_diag
zingale May 16, 2024
0df99a3
Merge branch 'development' into fix_omp_diag
biboyd May 17, 2024
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
136 changes: 43 additions & 93 deletions Source/MaestroDiag.cpp
View file
Open in desktop
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