Editing Appendix/Ramblings/HybridSchemeOld (section)

==Abbreviated Arguments==
This is the set of abbreviated arguments that I think should be used in &sect;2.3 of our "hybrid scheme" ApJ paper, that is, the paper by Z. D. Byerly, B. Adelstein-Lelbach, J. E. Tohline, &amp; D. C. Marcello (2014).

We evolve the following two coupled fluid equations,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} \rho + \nabla\cdot\rho\mathbf{u}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
0 \, ,
</math>
  </td>
</tr>
</table>
</div>
and,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} (\rho\mathbf{v}) + \nabla\cdot(\rho\mathbf{v}\mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
-\nabla p - \rho\nabla\Phi \, ,
</math>
  </td>
</tr>
</table>
</div>
where, <math>~\rho</math> is the mass density, <math>~p</math> is the gas pressure,  both <math>~\mathbf{v}</math> and <math>~\mathbf{u}</math> identify the same fluid velocity field (that is, <math>~\mathbf{v} = \mathbf{u}</math>), and <math>~\Phi</math> is the gravitational potential generated by a central point mass, <math>~M_\mathrm{pt}</math>, specifically,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\Phi
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~ -\frac{GM_\mathrm{pt}}{(x^2 + y^2 + z^2)^{1/2}} \, .
</math>
  </td>
</tr>
</table>
</div>
We use two different variables to represent the same velocity field to emphasize that, following Call et al. (2010), we have the freedom to choose different coordinate bases for each of the velocity terms that appear in the dyadic tensor product, <math>~\mathbf{v}\mathbf{u}</math>, in the nonlinear advection term of equation (7).  See Appendix A for further elaboration.

[<font color="red">NOTE:  We need to add an energy equation and an equation of state.  Make sure that all notation is completely consistent with the notation used in Dominic's Appendix B.</font>]

When rewriting the "momentum conservation" equation (7?) in terms of three orthogonal vector components, we begin by identifying two familiar sets of equations:  When advecting Cartesian momentum components &#8212; <math>~s_x \equiv \rho v_x</math>, <math>~ s_y \equiv \rho v_y</math>, <math>~ s_z \equiv \rho v_z</math> &#8212; we start with,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_x + \nabla\cdot (s_x \mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\hat{\mathbf{e}}_x \cdot \nabla p - ~\hat{\mathbf{e}}_x \cdot \rho \nabla \Phi \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_y + \nabla\cdot (s_y \mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\hat{\mathbf{e}}_y \cdot \nabla p - ~\hat{\mathbf{e}}_y \cdot \rho \nabla \Phi \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_z + \nabla\cdot (s_z \mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\hat{\mathbf{e}}_z \cdot \nabla p - ~\hat{\mathbf{e}}_z \cdot \rho \nabla \Phi \, ;
</math>
  </td>
</tr>
</table>
</div>
and, when advecting cylindrical momentum components &#8212; <math>~s_R \equiv \rho v_R</math>, <math>~ \ell_z \equiv R \rho v_\varphi</math>, <math>~ s_z \equiv \rho v_z</math> &#8212; the <math>~z</math>-component is identical to the Cartesian case but the other two orthogonal components are,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_R + \nabla\cdot (s_R \mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\hat{\mathbf{e}}_R \cdot \nabla p - ~\hat{\mathbf{e}}_R \cdot \rho \nabla \Phi 
+ \frac{\ell^2_z}{\rho R^3}\, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} \ell_z + \nabla\cdot (\ell_z \mathbf{u})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~R \hat{\mathbf{e}}_\varphi \cdot \nabla p - ~R\hat{\mathbf{e}}_\varphi \cdot \rho \nabla \Phi \, ,
</math>
  </td>
</tr>
</table>
</div>
where, <math>~R \equiv (x^2 + y^2)^{1/2}</math>.

These familiar sets of equations are morphed into the sets of equations used in our hybrid scheme by recognizing several things.  First, as is demonstrated in Appendix A, in all five identified momentum component equations, we can immediately replace <math>~\mathbf{u}</math> by the velocity field as viewed from a frame of reference that is rotating at angular frequency, <math>~\Omega</math>, namely,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\mathbf{u}'
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~\mathbf{u} - \hat\mathbf{e}_\varphi R\Omega \, ,
</math>
  </td>
</tr>
</table>
</div>
because <math>~\nabla\cdot (\hat\mathbf{e}_\varphi R\Omega) = 0</math>, that is, because the velocity field introduced by the frame rotation is divergence free.  All of the other elements of the five component equations remain unchanged when <math>~\mathbf{u}</math> is replaced by <math>~\mathbf{u}'</math> &#8212; in particular, all five advected quantities, <math>~s_x</math>, <math>~s_y</math>, <math>~s_z</math>, <math>~s_R</math>, and <math>~\ell_z</math>, still refer to components of the inertial-frame momentum (or angular-momentum) density. This recognition, alone, permits us to rewrite all three Cartesian components of the momentum conservation equation in the form that we have used for this project:
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_x + \nabla\cdot (s_x \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial x} p - \rho \frac{\partial}{\partial x} \Phi \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_y + \nabla\cdot (s_y \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial y} p - ~\rho \frac{\partial}{\partial y} \Phi \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_z + \nabla\cdot (s_z \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial z}  p - ~\rho \frac{\partial}{\partial z}  \Phi \, .
</math>
  </td>
</tr>
</table>
</div>

<table border="1" width="100%">
<tr>
  <th align="center">
Resolving Discrepancy Between Zach's Expressions and Mine
  </th>
</tr>
<tr><td align="left">
Now, let's see how first two of these equations gets modified if we shift things to advect linear momenta as measured in the rotating frame.  First, note that, in Cartesian Coordinates, 
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\mathbf{u}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~\mathbf{u}' + \hat\mathbf{e}_\varphi R\Omega = 
\mathbf{u}' + \biggl[ \hat\mathbf{j} \biggl(\frac{x}{R} \biggr) - \hat\mathbf{i} \biggl(\frac{y}{R} \biggr) \biggr] R\Omega \, .
</math>
  </td>
</tr>
</table>
</div>
Hence,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
~s_x = \rho u_x
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\rho( u_x' - y\Omega ) \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
~s_y = \rho u_y
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\rho( u_y' + x\Omega ) \, .
</math>
  </td>
</tr>
</table>
</div>
This means that the left-hand-side of the ''x''-momentum conservation equation becomes,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} [ \rho( u_x' - y\Omega ) ] + \nabla\cdot [ \rho( u_x' - y\Omega )  \mathbf{u}' ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\frac{\partial}{\partial t} [ \rho( u_x' ) ] + \nabla\cdot [ \rho( u_x')  \mathbf{u}' ]
- \frac{\partial}{\partial t} [ \rho( y\Omega ) ] - \nabla\cdot [ \rho( y\Omega )  \mathbf{u}' ]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\frac{\partial}{\partial t} s_x' + \nabla\cdot ( s_x'  \mathbf{u}' )
- y\Omega\biggl[ \frac{\partial\rho}{\partial t} + \nabla\cdot (\rho\mathbf{u}' )\biggr]
-\rho \Omega (\mathbf{u}' \cdot \nabla) y
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\frac{\partial}{\partial t} s_x' + \nabla\cdot ( s_x'  \mathbf{u}' ) -s_y' \Omega \, .
</math>
  </td>
</tr>
</table>
</div>
Similarly (but without detailing the steps), the left-hand-side of the ''y''-momentum conservation equation becomes,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} [ \rho( u_y' + x\Omega ) ] + \nabla\cdot [ \rho( u_y' + x\Omega )  \mathbf{u}' ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\frac{\partial}{\partial t} s_y' + \nabla\cdot ( s_y'  \mathbf{u}' ) + s_x' \Omega \, .
</math>
  </td>
</tr>
</table>
</div>
So, if rotating-frame quantities are advected on a rotating Cartesian grid, the governing set of equations is,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_x' + \nabla\cdot (s_x' \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial x} p - \rho \frac{\partial}{\partial x} \Phi  + s_y' \Omega\, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_y' + \nabla\cdot (s_y' \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial y} p - ~\rho \frac{\partial}{\partial y} \Phi - s_x' \Omega \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_z' + \nabla\cdot (s_z' \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial z}  p - ~\rho \frac{\partial}{\partial z}  \Phi \, .
</math>
  </td>
</tr>
</table>
</div>

</td></tr>
</table>

Second, we note that an evaluation of the advection term that appears on the left-hand-side of each component of the momentum equation, which is generically of the form,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\nabla\cdot(\Psi\mathbf{u}') \, ,
</math>
  </td>
</tr>
</table>
</div>
requires an assessment of the divergence of the three-dimensional flow field at each location on the computational grid.  But, in practice, it shouldn't matter whether this "assessment" is done on a Cartesian mesh or on a cylindrical mesh (or on any of a multitude of other mesh choices); the result should be the determination of the proper scalar value at every point on the chosen computational grid.  So, although the familiar form of the set of equations governing the time-rate-of-change of the cylindrical components of the momentum, presented above, was derived with the implicit assumption that each term would be evaluated on a cylindrical coordinate mesh, we can just as well evaluate the advection term on a Cartesian mesh.  This only requires that the divergence operator and the "transport" velocity, <math>~\mathbf{u}'</math>, be handled in exactly the same manner as they are handled when evaluating advection in the Cartesian set of equations. In the hybrid scheme being presented here, all simulations are conducted on a Cartesian mesh so, in all cases, the divergence operator and the transport velocity are broken down into Cartesian components before the advection term is evaluated.

Finally, because a Cartesian mesh is being adopted, the gradient operator on the right-hand-side of each component of the momentum equation is also explicitly broken down into its Cartesian components.  This means that, for our hybrid scheme, the right-hand-sides of equations (nn) and (mm) incorporate the operator projections,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\hat{\mathbf{e}}_R \cdot \nabla
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\biggl[ \hat{i}\biggr( \frac{x}{R} \biggr) + \hat{j} \biggr( \frac{y}{R} \biggr)\biggr] \cdot \nabla
= \frac{x}{R} \frac{\partial}{\partial x} + \frac{y}{R} \frac{\partial}{\partial y } \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
~\hat{\mathbf{e}}_\varphi \cdot \nabla
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\biggl[ \hat{j}\biggr( \frac{x}{R} \biggr) - \hat{i} \biggr( \frac{y}{R} \biggr)\biggr] \cdot \nabla
= \frac{y}{R} \frac{\partial}{\partial x} - \frac{x}{R} \frac{\partial}{\partial y } \, .
</math>
  </td>
</tr>
</table>
</div>

With all of these recognitions in hand, in our hybrid scheme the three components of the cylindrical momentum equations are,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_R + \nabla\cdot (s_R \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{1}{R} \biggl( x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} \biggr) p 
- \frac{\rho}{R} \biggl( x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} \biggr) \Phi 
+ \frac{\ell^2_z}{\rho R^3} \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} \ell_z + \nabla\cdot (\ell_z \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\biggl( y\frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr) p 
+ \rho \biggl( y\frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr)\Phi \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\frac{\partial}{\partial t} s_z + \nabla\cdot (s_z \mathbf{u}')
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- ~\frac{\partial}{\partial z}  p - ~\rho \frac{\partial}{\partial z}  \Phi \, .
</math>
  </td>
</tr>
</table>
</div>

As discussed by Call, Tohline, &amp; Lehner (2010), it is noteworthy that the right-hand-side of the hybrid-scheme equation that governs transport (and conservation) of angular momentum does not contain a Coriolis term.  This is because <math>~\ell_z</math>, the quantity being advected and tracked, is the angular momentum density as measured in the ''inertial'' frame of reference.  As is demonstrated in &sect;A.4 of Appendix A, a Coriolis term arises if the equation is written in a form where the quantity being advected is the ''rotating-frame'' angular momentum density.  This equation, which contains a Coriolis term, is more familiar to the astrophysics community &#8212; see, for example, Norman &amp; Wilson (1980) and New &amp; Tohline (1987).  But, for purposes of angular momentum conservation, we consider it to be far preferable to adopt a version of the equation in which the velocity does not explicitly appear in the source term.

Finally, we note that <math>~s_R</math> and <math>~\ell_z</math> can be straightforwardly expressed in terms of Cartesian components of <math>~\mathbf{u}</math> or <math>~\mathbf{u}'</math>.  Specifically, remembering that <math>~\mathbf{u} = \mathbf{v}</math>,
<div align="center">
<table border="0" cellpadding="3">

<tr>
  <td align="right">
<math>
~s_R \equiv \rho v_R
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\frac{\rho}{R} (xu_x + y u_y) = \frac{\rho}{R} (xu_x' + y u_y') \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
~\ell_z \equiv R\rho v_\varphi
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\rho (xu_y - y u_x) = \rho (xu_y' - y u_x')  + \rho\Omega_0 (x^2 + y^2) \, .
</math>
  </td>
</tr>
</table>
</div>