Is hydro diffusion doing the right thing in spherical coordinates?

[c-white]

I was looking at the 3 locations in SphericalPolar::CoordSrcTerms() where hydro diffusion is accounted for. I don't know much about the diffusion implementation, but it struck me that some of the fluxes being used might not exist in certain cases, and at other times terms might not be accounted for.

The r-momentum is apparently affected by theta-fluxes of theta-momentum and phi-fluxes of phi-momentum:

athena/src/coordinates/spherical_polar.cpp

Lines 526 to 529 in 5fd9ea5

	if (do_hydro_diffusion) {
	m_ii += 0.5*(phd->visflx[X2DIR](IM2,k,j+1,i)+phd->visflx[X2DIR](IM2,k,j,i));
	m_ii += 0.5*(phd->visflx[X3DIR](IM3,k+1,j,i)+phd->visflx[X3DIR](IM3,k,j,i));
	}

But in 1D or 2D will these always be set? Similarly for the theta-momentum:

athena/src/coordinates/spherical_polar.cpp

Lines 554 to 555 in 5fd9ea5

	if (do_hydro_diffusion)
	m_pp += 0.5*(phd->visflx[X3DIR](IM3,k+1,j,i)+phd->visflx[X3DIR](IM3,k,j,i));

I suppose those arrays could just be 0 if the dimensionality is low enough, so the right thing is always done. I'm more confused by the phi-momentum source term:

athena/src/coordinates/spherical_polar.cpp

Lines 560 to 573 in 5fd9ea5

	if (use_x2_fluxes) {
	u(IM3,k,j,i) -= dt*coord_src1_i_(i)*coord_src2_j_(j)*
	(coord_area2_j_(j)*flux[X2DIR](IM3,k,j,i)
	+ coord_area2_j_(j+1)*flux[X2DIR](IM3,k,j+1,i));
	} else {
	Real m_ph = prim(IDN,k,j,i) * prim(IM3,k,j,i) * prim(IM2,k,j,i);
	if (MAGNETIC_FIELDS_ENABLED) {
	m_ph -= bcc(IB3,k,j,i) * bcc(IB2,k,j,i);
	}
	if (do_hydro_diffusion)
	m_ph += 0.5*(phd->visflx[X2DIR](IM3,k,j+1,i)+phd->visflx[X2DIR](IM3,k,j,i));
	u(IM3,k,j,i) -= dt*coord_src1_i_(i)*coord_src3_j_(j)*m_ph;
	}

Since use_x2_fluxes is true if and only if we are in 2D or 3D, it seems this term only gets added in 1D, in which case the referenced fluxes are 0 (going by the above logic). And it seems like the theta-phi viscous stress is not properly being accounted for in 2D and 3D.

Again, I'm not really familiar with how diffusion works in the code, so maybe this is all fine. But I'm working on a project that involves taking a very close look at curvilinear source terms, so I'm trying to understand exactly what's happening in this function.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[pdmullen]

Can we re-open this issue @felker? Stale-bot closed this on Feb 6, 2019; but @c-white raises some good questions here that were never addressed. I'd raise another one:

• When computing the "geometric source term" when viscosity is enabled, hd.viscflx is used in some places:

  athena/src/coordinates/spherical_polar.cpp

  Lines 496 to 497 in 0092ae0

  	if (do_hydro_diffusion)
  	m_pp += 0.5*(hd.visflx[X3DIR](IM3,k+1,j,i) + hd.visflx[X3DIR](IM3,k,j,i));

  
  But in others terms flux is used:

  athena/src/coordinates/spherical_polar.cpp

  Lines 505 to 506 in 0092ae0

  	(coord_area2_j_(j)*flux[X2DIR](IM3,k,j,i)
  	+ coord_area2_j_(j+1)*flux[X2DIR](IM3,k,j+1,i));

Remember that in the explicit integration of diffusive physics, flux contains the sum of standard (M)HD fluxes and diffusive fluxes. Do these geometric source terms accurately reflect this? This also shows up in cylindrical coordinates.

Aside: When producing the STS extension, I assumed that all geometric source terms (when diffusion is enabled) were correct; so if any geometric issues with explicit diffusion exist, they would have also been carried over into STS. My testing of curvilinear STS was limited in #299.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.