Remove division-by-zero bug for viscosity + 3d polar boundary condition

[apbailey]

Description

This is a simple fix for a divide-by-zero I stumbled upon. I think it is just a missed edge case when viscosity is used in conjunction with 3d spherical polar extending to the pole x2min=0.0. The offending code is line 398 of viscosity.cpp:

void HydroDiffusion::FaceYdz(const int k, const int j, const int il, const int iu,
                             const AthenaArray<Real> &prim, AthenaArray<Real> &len) {
  if (pmb_->pmy_mesh->f3) {
#pragma omp simd
    for (int i=il; i<=iu; ++i) {
      len(i) = pco_->h32f(j)
               * ( prim(IM3,k,j,i)/pco_->h32v(j) - prim(IM3,k,j-1,i)/pco_->h32v(j-1) )
               / pco_->h2v(i) / pco_->dx2v(j-1)
               // KGF: add the off-centered quantities first to preserve FP symmetry
               + 0.5*(    (prim(IM2,k+1,j,i) + prim(IM2,k+1,j-1,i))
                          - (prim(IM2,k-1,j,i) + prim(IM2,k-1,j-1,i)) )
               / pco_->h31v(i)
               / pco_->h32f(j) / (pco_->dx3v(k-1) + pco_->dx3v(k));
    }

When computing the FaceYdz viscous flux for active cells at the pole pco_->h32f(j) occurring in the denominator evaluates to zero at the pole. Adding the flux divergence propagates this NaN to the first active cell which spreads throughout the domain with each timestep.

This is just a simple fix that replaces pco_->h32f(j) in the denominator with a conditional that returns 1.0 if h32f(j)=0.0 and returns h32f(j) otherwise. The value the conditional returns when h32f(j)=0.0 is irrelevant since this flux gets multiplied by a zero area face later, it just needs to be non-NaN so as not to persist. The case for the conditional to return 1.0 is only triggered in 3D spherical polar as h32f is always 1.0 for cartesian and cylindrical and the scope of this block of code is limited to 3D.

Testing and validation

tested that a 3d polar run inviscid and vanishingly small nu_iso=1.0e-10 give similar and now non-NaN results. 10 out of 10 existing diffusion regression tests passed on laptop.

[buildbot-princeton]

Can one of the admins verify this patch?

[felker]

ok to test

[felker]

Nice find! I get a little nervous seeing a logical operator applied to a floating point value like this, and might prefer something more expressive, but I am not confident about a better alternative.

[apbailey]

feel similarly and can edit with any suggestions. mostly just calling attention to it with this PR and since this is a particular edge case tried to make my change as non-invasive as possible (which !float logic seemed not to be too terrible for -- needing to meet the most hyper-specific of criteria to do anything different from what the code is already doing)

[c-white]

Would the Coordinates::IsPole function be of use here? This could also catch the south pole, but maybe that's not a problem since sin(3.141592653589793) == 1.2246467991473532e-16, and the numerator will never be so large as to make the quotient overflow.

[apbailey]

Good suggestion, I think IsPole() is a better, more expressive way of doing things and looks to give reasonable results on my simple 3d viscous polar problem so I've updated my branch with it