fix divu when theta=0 or pi (#2975)

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

@@ -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;