Update euler_grid.cc

src/core/diffusion/euler_grid.cc

@@ -244,12 +244,15 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) {
           b = c - nx * ny;
           t = c + nx * ny;
 
-          real_t left = c1_[c - 1];
-          real_t right = c1_[c + 1];
-          real_t bottom = c1_[b];
-          real_t top = c1_[t];
-          real_t north = c1_[n];
-          real_t south = c1_[s];
+          // Clamp to avoid out of bounds access. Clamped values are initialized
+          // to a wrong value but will be overwritten by the boundary condition
+          // evaluation. All other values are correct.
+          real_t left{c1_[std::clamp(l, size_t{0}, num_boxes - 1)]};
+          real_t right{c1_[std::clamp(r, size_t{0}, num_boxes - 1)]};
+          real_t bottom{c1_[std::clamp(b, size_t{0}, num_boxes - 1)]};
+          real_t top{c1_[std::clamp(t, size_t{0}, num_boxes - 1)]};
+          real_t north{c1_[std::clamp(n, size_t{0}, num_boxes - 1)]};
+          real_t south{c1_[std::clamp(s, size_t{0}, num_boxes - 1)]};
           real_t center_factor{6.0};
 
           if (x == 0 || x == (nx - 1) || y == 0 || y == (ny - 1) || z == 0 ||
@@ -335,12 +338,15 @@ void EulerGrid::DiffuseWithPeriodic(real_t dt) {
           b = c - nx * ny;
           t = c + nx * ny;
 
-          real_t left = c1_[l];
-          real_t right = c1_[r];
-          real_t bottom = c1_[b];
-          real_t top = c1_[t];
-          real_t north = c1_[n];
-          real_t south = c1_[s];
+           // Clamp to avoid out of bounds access. Clamped values are initialized
+          // to a wrong value but will be overwritten by the boundary condition
+          // evaluation. All other values are correct.
+          real_t left{c1_[std::clamp(l, size_t{0}, num_boxes - 1)]};
+          real_t right{c1_[std::clamp(r, size_t{0}, num_boxes - 1)]};
+          real_t bottom{c1_[std::clamp(b, size_t{0}, num_boxes - 1)]};
+          real_t top{c1_[std::clamp(t, size_t{0}, num_boxes - 1)]};
+          real_t north{c1_[std::clamp(n, size_t{0}, num_boxes - 1)]};
+          real_t south{c1_[std::clamp(s, size_t{0}, num_boxes - 1)]};
 
           if (x == 0 || x == (nx - 1) || y == 0 || y == (ny - 1) || z == 0 ||
               z == (nz - 1)) {