From 7b80efd066e7b82e567e5d95837341e01f319960 Mon Sep 17 00:00:00 2001
From: Zhi Chen <62574124+zhichen3@users.noreply.github.com>
Date: Mon, 14 Oct 2024 16:16:14 -0400
Subject: [PATCH] fix divu when theta=0 or pi (#2975)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

When calculating vy, when 𝜃 = 0 or or 𝜃 = 𝜋, sinc=0 which causes division by 0. But the numerator should cancel due to azimuthal symmetry. So numerator is 0. So set vy=0 when that happens.
---
 Source/hydro/advection_util.cpp | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp
index e17e4d0326..27d3befa3c 100644
--- a/Source/hydro/advection_util.cpp
+++ b/Source/hydro/advection_util.cpp
@@ -310,12 +310,18 @@ Castro::divu(const Box& bx,
         // Finite difference to get divergence. ux = 1/r^2 d/dr(r^2 * u)
         ux = (ur * rr * rr - ul * rl * rl) * dxinv / (rc * rc);
 
-        // These are transverse averages in the x-direction
-        Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV));
-        Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV));
+        // If sinc == 0, then vy goes inf.
+        // But due to Phi-symmetry, vt*sint = vb*sinb, so set to 0.
+        if (sinc == 0.0_rt) {
+            vy = 0.0_rt;
+        } else {
+            // These are transverse averages in the x-direction
+            Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV));
+            Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV));
 
-        // Finite difference to get divergence. vy = 1/(r sin) * d/dtheta(v sin)
-        vy = (vt * sint - vb * sinb) * dyinv / (rc * sinc);
+            // Finite difference to get divergence. vy = 1/(r sin) * d/dtheta(v sin)
+            vy = (vt * sint - vb * sinb) * dyinv / (rc * sinc);
+        }
     }
 
     div(i,j,k) = ux + vy;