From b0ad43b545ae7d50a27ad2a2caf51b7094576c19 Mon Sep 17 00:00:00 2001 From: Zhi Date: Fri, 11 Oct 2024 23:54:56 -0400 Subject: [PATCH] fix divu when theta=0 or pi --- 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;