Spherical Coordinate Source Term

[hyluoo]

Hi all,
In which file/part of the code can I see how the geometric source is done in spherical coordinate?
In particular, the equation 19 in the Athena++ method paper is confusing for me and I would like to see how it is implemented in practice.

[felker]

athena/src/hydro/srcterms/rotating_system_srcterms.cpp

Lines 58 to 74 in e58b549

	} else if(std::strcmp(COORDINATE_SYSTEM, "spherical_polar") == 0) {
	if (pmb->block_size.nx2 == 1) {
	std::stringstream msg;
	msg << "### FATAL ERROR in HydroSourceTerms::RotatingSystemSourceTerms" << std::endl
	<< "Rotating System does not support spherical polar coordinates in 1D."
	<< std::endl;
	ATHENA_ERROR(msg);
	}
	// dM1/dt = 2 \rho vc vp / r
	// +\rho (vc)^2 / r
	// dM2/dt = 2 \rho vp cot(\theta) vc / r
	// + \rho cot(\theta) (vc)^2 /r
	// dM3/dt = -\rho vr (2 vc)/r
	// -\rho vt (2 vc) cot(\theta) /r
	// dE/dt = \rho vr (vc)^2/r
	// + \rho vt (vc)^2 cot(\theta) /r
	// vc = r sin(\theta)\Omega

In the future, I recommend searching the source code for formulas using something like grep -nri "cot("

[c-white]

The above is for a rotating (i.e. noninertial) source term. The base case of the particular discretization of the spherical divergence operator, whether or not the reference frame is rotating, is found in the appropriate coordinates file:

athena/src/coordinates/spherical_polar.cpp

Lines 442 to 526 in e58b549

	void SphericalPolar::AddCoordTermsDivergence(const Real dt, const AthenaArray<Real> *flux,
	const AthenaArray<Real> &prim,
	const AthenaArray<Real> &bcc, AthenaArray<Real> &u) {
	Real iso_cs = pmy_block->peos->GetIsoSoundSpeed();
	bool use_x2_fluxes = pmy_block->block_size.nx2 > 1;
	HydroDiffusion &hd = pmy_block->phydro->hdif;
	bool do_hydro_diffusion = (hd.hydro_diffusion_defined &&
	(hd.nu_iso > 0.0 || hd.nu_aniso > 0.0));
	// Go through cells
	for (int k=pmy_block->ks; k<=pmy_block->ke; ++k) {
	for (int j=pmy_block->js; j<=pmy_block->je; ++j) {
	#pragma omp simd
	for (int i=pmy_block->is; i<=pmy_block->ie; ++i) {
	// src_1 = < M_{theta theta} + M_{phi phi} ><1/r>
	Real m_ii = prim(IDN,k,j,i)*(SQR(prim(IM2,k,j,i)) + SQR(prim(IM3,k,j,i)));
	if (NON_BAROTROPIC_EOS) {
	m_ii += 2.0*prim(IEN,k,j,i);
	} else {
	m_ii += 2.0*(iso_cs*iso_cs)*prim(IDN,k,j,i);
	}
	if (MAGNETIC_FIELDS_ENABLED) {
	m_ii += SQR(bcc(IB1,k,j,i));
	}
	if (!STS_ENABLED) {
	if (do_hydro_diffusion) {
	m_ii += 0.5*(hd.visflx[X2DIR](IM2,k,j+1,i) + hd.visflx[X2DIR](IM2,k,j,i));
	m_ii += 0.5*(hd.visflx[X3DIR](IM3,k+1,j,i) + hd.visflx[X3DIR](IM3,k,j,i));
	}
	}
	u(IM1,k,j,i) += dt*coord_src1_i_(i)*m_ii;
	// src_2 = -< M_{theta r} ><1/r>
	u(IM2,k,j,i) -= dt*coord_src2_i_(i)*
	(coord_area1_i_(i)*flux[X1DIR](IM2,k,j,i)
	+ coord_area1_i_(i+1)*flux[X1DIR](IM2,k,j,i+1));
	// src_3 = -< M_{phi r} ><1/r>
	u(IM3,k,j,i) -= dt*coord_src2_i_(i)*
	(coord_area1_i_(i)*flux[X1DIR](IM3,k,j,i)
	+ coord_area1_i_(i+1)*flux[X1DIR](IM3,k,j,i+1));
	// src_2 = < M_{phi phi} ><cot theta/r>
	Real m_pp = prim(IDN,k,j,i)*SQR(prim(IM3,k,j,i));
	if (NON_BAROTROPIC_EOS) {
	m_pp += prim(IEN,k,j,i);
	} else {
	m_pp += (iso_cs*iso_cs)*prim(IDN,k,j,i);
	}
	if (MAGNETIC_FIELDS_ENABLED) {
	m_pp += 0.5*( SQR(bcc(IB1,k,j,i)) + SQR(bcc(IB2,k,j,i))
	- SQR(bcc(IB3,k,j,i)) );
	}
	if (!STS_ENABLED) {
	if (do_hydro_diffusion)
	m_pp += 0.5*(hd.visflx[X3DIR](IM3,k+1,j,i) + hd.visflx[X3DIR](IM3,k,j,i));
	}
	u(IM2,k,j,i) += dt*coord_src1_i_(i)*coord_src1_j_(j)*m_pp;
	// src_3 = -< M_{phi theta} ><cot theta/r>
	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 (!STS_ENABLED) {
	if (do_hydro_diffusion)
	m_ph += 0.5*(hd.visflx[X2DIR](IM3,k,j+1,i) + hd.visflx[X2DIR](IM3,k,j,i));
	}
	u(IM3,k,j,i) -= dt*coord_src1_i_(i)*coord_src3_j_(j)*m_ph;
	}
	}
	}
	}
	return;
	}

[hyluoo]

Thank you Chris that's exactly what I was looking for

[hyluoo]

I checked the code, the source term should be like this, is the one in paper a typo?