Merge branch 'develop' into lroberts36/refactor-meshblockdata-initial…

…ization
parthenon-hpc-lab · Jun 15, 2024 · 11324a9 · 11324a9
2 parents 04d8e1e + 3ad98d1
commit 11324a9
Show file tree

Hide file tree

Showing 3 changed files with 7 additions and 3 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -27,6 +27,7 @@
 
 ### Fixed (not changing behavior/API/variables/...)
 - [[PR 1092]](https://github.com/parthenon-hpc-lab/parthenon/pull/1092) Updates to DataCollection and MeshData to remove requirement of predefining MeshBlockData
+- [[PR 1113]](https://github.com/parthenon-hpc-lab/parthenon/pull/1113) Prevent division by zero
 - [[PR 1112]](https://github.com/parthenon-hpc-lab/parthenon/pull/1112) Remove shared_ptr cycle in forest::Tree
 - [[PR 1104]](https://github.com/parthenon-hpc-lab/parthenon/pull/1104) Fix reading restarts due to hidden ghost var
 - [[PR 1098]](https://github.com/parthenon-hpc-lab/parthenon/pull/1098) Move to symmetrized logical coordinates and fix SMR bug

diff --git a/example/poisson/poisson_package.cpp b/example/poisson/poisson_package.cpp
@@ -22,6 +22,7 @@
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 #include <solvers/solver_utils.hpp>
+#include <utils/robust.hpp>
 
 #include "defs.hpp"
 #include "kokkos_abstraction.hpp"
@@ -296,7 +297,8 @@ TaskStatus CheckConvergence(T *u, T *du) {
       parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0,
       v.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &eps) {
-        Real reps = std::abs(dv(b, idphi, k, j, i) / v(b, iphi, k, j, i));
+        Real reps = std::abs(
+            parthenon::robust::ratio(dv(b, idphi, k, j, i), v(b, iphi, k, j, i)));
         Real aeps = std::abs(dv(b, idphi, k, j, i));
         eps = std::max(eps, std::min(reps, aeps));
       },

diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
@@ -25,8 +25,8 @@
 #include "interface/state_descriptor.hpp"
 #include "kokkos_abstraction.hpp"
 #include "solvers/solver_utils.hpp"
-
 #include "tasks/tasks.hpp"
+#include "utils/robust.hpp"
 
 namespace parthenon {
 
@@ -246,7 +246,8 @@ class MGSolver {
                   DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints - 1, [&](const int idx) {
                     const Real off_diag = Ax[idx] - diag[idx] * xo[idx];
                     const Real val = prhs[idx] - off_diag;
-                    xn[idx] = weight * val / diag[idx] + (1.0 - weight) * xo[idx];
+                    xn[idx] =
+                        weight * robust::ratio(val, diag[idx]) + (1.0 - weight) * xo[idx];
                   });
             }
           });