Editing PGE/AStarScheme (section)

===The So-called A* Scheme===
Now, in a rotating frame of reference, this preferred form of the azimuthal component of the equation of motion takes the form,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t(\rho \varpi u_\theta) + \nabla\cdot(\rho\varpi u_\theta \boldsymbol{u} )
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~ - \varpi \nabla_\theta P - \varpi \rho \nabla_\theta \Phi -2\rho\varpi \Omega_f u_\varpi \, .
</math>
  </td>
</tr>
</table>
</div>
Notice that a new term has appeared due to the coriolis force.  Traditionally, in numerical hydrodynamics, this new term has been treated explicitly as a source term.  Hence, this component of the equation of motion no longer takes on a strictly conservative form, and the adopted "conservative" finite-difference is no longer a particularly useful tool to guarantee that the angular momentum is globally conserved even when the external forces due to pressure and gravity balance one another.

To derive a form of this equation that is a lot more suited to a "conservative" finite-difference implementation, note that,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t (\rho \varpi^2 \Omega_f ) + \nabla\cdot [(\rho \varpi^2 \Omega_f)\boldsymbol{u}]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
\varpi^2 \Omega_f [ \partial_t \rho + \nabla\cdot (\rho \boldsymbol{u})] + \rho \Omega_f (\boldsymbol{u}\cdot \nabla)\varpi^2
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
2\rho\varpi \Omega_f u_\varpi \, .
</math>
  </td>
</tr>

</table>
</div>
Hence, in a rotating frame of reference, the azimuthal component of the equation of motion can be written as,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t (\rho\varpi u_\theta) + \nabla\cdot(\rho\varpi u_\theta \boldsymbol{u} + \partial_t (\rho\varpi^2 \Omega_f) + 
\nabla\cdot [ (\rho\varpi^2 \Omega_f) \boldsymbol{u} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~ - \varpi \nabla_\theta P - \varpi \rho \nabla_\theta\Phi \, ,
</math>
  </td>
</tr>
</table>
</div>
or,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t [ \rho\varpi (u_\theta + \varpi\Omega_f) + \nabla\cdot \{ [ \rho\varpi (u_\theta + \varpi\Omega_f) ] \boldsymbol{u} \}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~ - \varpi \nabla_\theta P - \varpi \rho \nabla_\theta\Phi \, .
</math>
  </td>
</tr>
</table>
</div>
When advancing the angular momentum density (i.e., <math>~A_\mathrm{rot} = \rho\varpi u_\theta</math>) forward in time using a finite-difference scheme, I recommend that the "sourcing" step only include the terms on the right-hand-side of these last expressions (i.e., only the gradients in the pressure and gravitational potential), and the "fluxing" step should include the following terms:
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\partial_t A_\mathrm{rot}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \nabla\cdot (A_\mathrm{rot} \boldsymbol{u}) - \partial_t (\rho\varpi^2 \Omega_f) - \nabla\cdot [ ( \rho \varpi^2 \Omega_f) \boldsymbol{u} ]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \nabla\cdot (A_\mathrm{rot} \boldsymbol{u}) 
- \varpi^2 \Omega_f \partial_t (\rho ) - \Omega_f \nabla\cdot [(\rho\varpi^2) \boldsymbol{u} ]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \nabla\cdot (A_\mathrm{rot} \boldsymbol{u}) 
+ \varpi^2 \Omega_f \nabla\cdot [(\rho) \boldsymbol{u}] - \Omega_f \nabla\cdot [(\rho\varpi^2) \boldsymbol{u} ] \, .
</math>
  </td>
</tr>
</table>
</div>
This last expression is directly implementable using our standard fluxing scheme because all three terms have the form <math>~\nabla\cdot[Q\boldsymbol{u}]</math>.  (Note that the density that appears in the last two terms on the right-hand-side of this last expression must be taken from precisely the same point in time as the "<math>~A_\mathrm{rot}</math>" that appears in the first term on the right-hand-side.)

Whether you adopt precisely this final prescription or not, the primary point to keep in mind is that you want to advect the "intertial" [''sic''] angular momentum density <math>~[\rho\varpi(u_\theta + \varpi\Omega_f)] = [A_\mathrm{rot} + \rho\varpi^2 \Omega_f]</math> using the "rotating-frame velocity <math>~\boldsymbol{u}</math> at each grid cell face.  Hence, you might prefer to use one slightly earlier relation to guide your design of the fluxing step.  In the absence of "true" source terms (due to pressure and gravity), we have,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t [\rho \varpi (u_\theta + \varpi\Omega_f)]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \nabla\cdot\{[\rho\varpi ( u_\theta + \varpi\Omega_f)]\boldsymbol{u} \} \, ,
</math>
  </td>
</tr>
</table>
</div>
which is the same as,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t [A_\mathrm{rot} + \rho\varpi^2\Omega_f ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \nabla\cdot\{[A_\mathrm{rot} + \rho\varpi^2\Omega_f ]\boldsymbol{u} \} \, .
</math>
  </td>
</tr>
</table>
</div>
But this may also be written as,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\partial_t (A_\mathrm{rot})
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \partial_t(\rho\varpi^2\Omega_f) - \nabla\cdot\{[A_\mathrm{rot} + \rho\varpi^2\Omega_f ]\boldsymbol{u} \} \, .
</math>
  </td>
</tr>
</table>
</div>
So, you can carry out the calculation by first adding the quantity <math>~(\rho\varpi^2 \Omega_f)</math> to <math>~A_\mathrm{rot}</math> to get the value of the angular momentum density in the inertial frame; advect this "inertial" quantity; then subtract from this result ''the change in the quantity'' <math>~(\rho\varpi^2 \Omega_f)</math> as determined from the updated value of the density as derived via the continuity equation.