Cell Centered vs Face Centered Potentials

[hwhitehead]

I have some questions concerning the implementation of energy changes in Athena++.

When calculating energy changes associated with the movement of the fluid in some field, the work done is calculated by averaging over the faces of a cell. For self-gravity, I believe Athena++ uses Mullen et al in which the derivatives of the potential are evaluated at the faces, using the cell-centered potential values in both the parent cell $(i,j,k)$ and its neighbour $(i+1,k,j)$.

These derivatives $g$ are then used along with the Riemann mass fluxes at the cell faces to calculate the work done. Below is the source term due to self-gravity in Mullen et al, where the last row indicates the source to energy.

See the use of $g$ evaluated at the cell faces throughout. This form appears to be supported by self_gravity.cpp, if it is the case that the array phi is defined at the cell centers (hard to confirm without better knowledge of the code).

athena/src/hydro/srcterms/self_gravity.cpp

Lines 43 to 50 in a94f18f

	Real dx1 = pmb->pcoord->dx1v(i);
	Real hdtodx1 = 0.5*dt/dx1;
	Real dpl = -(pgrav->phi(k,j,i ) - pgrav->phi(k,j,i-1));
	Real dpr = -(pgrav->phi(k,j,i+1) - pgrav->phi(k,j,i ));
	cons(IM1,k,j,i) += hdtodx1 * prim(IDN,k,j,i) * (dpl + dpr);
	if (NON_BAROTROPIC_EOS)
	cons(IEN,k,j,i) += hdtodx1 * (flux[X1DIR](IDN,k,j,i ) * dpl
	+ flux[X1DIR](IDN,k,j,i+1) * dpr);

What I am trying to work out is why the energy method in the shearing routine shearing_box.cpp appears to be different. Notably, the evaluation of the acceleration appears to be done not with the cell-centered potentials, but the face centered values.

athena/src/hydro/srcterms/shearing_box.cpp

Lines 56 to 62 in a94f18f

	if (NON_BAROTROPIC_EOS) {
	const Real phic = qO2*SQR(xc);
	const Real phil = qO2*SQR(pmb->pcoord->x1f(i));
	const Real phir = qO2*SQR(pmb->pcoord->x1f(i+1));
	cons(IEN,k,j,i) += dt*(flux[X1DIR](IDN,k,j,i)*(phic-phil)+
	flux[X1DIR](IDN,k,j,i+1)*(phir-phic))
	/pmb->pcoord->dx1f(i);

Here phil is not the potential evaluated in the center of the left-side neighbour (which would use x1v(i-1)) but instead uses the left-face-centered value of the parent cell. This is why a factor $\frac{1}{2}$ does not feature in this calculation, but did in the self-gravity method, as the points at which the potentials are being evaluated are only $\frac{1}{2}\Delta x$ apart, not $\Delta x$. However this means that the derivative of the potential is being evaluated not at the cell faces, but at some intermediate point between the face and cell center.

Is there a reason for this change in evaluation location? Both methods should produce similar results, but as they are technically different it would be interesting to know why such differences are required.

[tomidakn]

The self-gravity potential is defined at cell centers.

Regarding the shearing box, all I know is this is inherited from Athena, and its original paper is here:
https://iopscience.iop.org/article/10.1088/0067-0049/189/1/142/pdf
Maybe @jmstone remembers its rationale.

[apbailey]

I can't speak to self-gravity, but to add, the shearing box source is formulated in such a way as to conserve the total energy with cell centered-phi. As you posted the update to $E$ in the code is $\Delta t\left[F_L(\Phi_C - \Phi_L) + F_R(\Phi_R - \Phi_C)\right]/\Delta x$ (L=left face, C=center, R=right face) and the update to $\rho$ is given by the riemann solver $=\Delta t(F_R - F_L)/\Delta x$. So that the update to $E+\rho \Phi_C$ is $\Delta t\left[F_R\Phi_R - F_L\Phi_L\right]/\Delta x$ and the quantity $E+\rho \Phi_C$ is only changed by flux at the boundaries as one would expect for a conservative formulation. Note the similarity to Mullen eq 53. Similar to your question I was curious why pointmass.cpp does not use this conservative potential source formulation that was used by athena c-version, I suppose it just works well enough regardless