Appendix/Ramblings/HybridSchemeOld
Hybrid Advection Scheme
Preface
March 1, 2014 by Joel E. Tohline
<full preface> … This seemed too good to be true. The discovered code modification would allow us to conserve angular momentum very accurately and, at the same time, allow us to use a rotating grid and thereby minimize numerical diffusion … <read more>
In his dissertation research, Jay Call (see Call, Tohline, & Lehner 2010) derived a complete description of this hybrid advection scheme in a fully relativistic and generalized coordinate framework. He showed that one can write the system of fluid equations in a manner that facilitates advection of inertialframe quantities across a rotating grid. In addition — and more importantly — he showed how to write the system of fluid equations to allow advection of inertialframe angular momentum (generally associated with a cylindrical coordinate mesh) across a rotating Cartesian grid. In the discussion that follows here, we derive the Newtonian version of Jay's hybrid scheme.
Setting the Stage
Recognizing Statements of Conservation
When dealing with the compressible fluid equations, we will often encounter hyperbolic PDEs of the following form:



where we are using to represent the velocity field of the fluid as viewed from an inertial frame of reference, and the total (as opposed to partial) time derivative indicates the timerate of change of is being measured in a socalled Lagrangian fashion, that is, at the location of some fluid element and moving along with that fluid element.
When we encounter a situation in which the "source" term, , on the righthand side is zero, we will be able to identify the scalar variable, , as the volume density of some conserved quantity. For example, the continuity equation — which is a mathematical statement of mass conservation — has the form,
where, is the mass per unit volume or, simply, the mass density of the fluid element. Clearly, when the mass of a Lagrangian fluid element is conserved, the fluid element's mass density changes only in accordance with the divergence of the local velocity field.
Similarly, if we are following the evolution of a fluid that expands and contracts adiabatically, we should expect to encounter an equation of the form,



or, equivalently, 



where, is the entropy density of a Lagrangian fluid element. Or, if an axisymmetric distribution of fluid is moving in an axisymmetric potential, we should expect the azimuthal component of the fluid's angular momentum to be conserved, in which case we should expect to encounter a dynamical equation of the form,



where, is the Lagrangian fluid element's (cylindrical radial) distance measured from the symmetry axis of the underlying potential and is the azimuthal component of the inertial velocity field, , at the location of the fluid element.
Alternative Reference Frames
Now, we might want to examine the timedependent behavior of a fluid system while viewing the flow from a reference frame that is more or less moving along with the fluid. This new frame of reference need not be an inertial frame; for example, when studying a rotating fluid, we may want to view the system's evolution from a rotating frame of reference. This will be accomplished mathematically by adjusting the dynamical equations so that the velocity that appears in the divergence term accounts for the new "frame" velocity field; specifically, we want to replace with,



(Here, we will consider only timeindependent functional expressions for the frame velocity, .) Of course, switching to the rotating frame must be done in such a way that the resulting, new PDE describes exactly the same physical behavior of the system as was described by the original equation; that is, the new equation must be derivable from the original one.
If is a divergencefree velocity field, then the transformation is trivial. For example, if we choose a frame of reference that is rotating uniformly with angular velocity, , then,



and, utilizing cylindrical coordinates,



Hence,





so the new generic hyperbolic PDE becomes,



and we can be confident that this new PDE represents the physics of the system just as well as the original PDE.
Eulerian Representation
We can shift any of the PDEs from a Lagrangian to an Eulerian representation — and thereby use them to follow the timerate of change of physical variables at a point in space that is fixed with respect to the chosen frame of reference — by using the following transformation to replace each total time derivative with a partial time derivative:



Hence, the "new" generic hyperbolic PDE derived above can be rewritten as,



or, more succinctly,



This equation also is broadly recognized as a conservation statement because, when , the variable will represent the volume density of a conserved quantity.
We should emphasize that the inertialframe version of this Eulerian conservation equation can be retrieved straightforwardly by setting , which is equivalent to setting . It is,



The physics of the flow that is being described by this PDE is identical to the physics that is described by the preceding PDE. But an important distinction must be made regarding how the two equations are interpreted. The "inertial frame" version of the equation (explicitly containing ) provides the timerate of change of at a fixed point in inertial space, while the "new" version (explicitly containing ) provides the timerate of change of at a fixed point in our "new" rotating coordinate frame.
Angular Momentum Conservation
When the three vector components of the Euler equation (of motion) are projected onto a nonrotating cylindrical coordinate grid, the azimuthal component of the Euler equation may be written as,



For this equation, the source term is identified as,



and is the inertialframe angular momentum density, as measured with respect to the coordinate axis. This corresponds the scalar equation and representation referred to as "Case B ()" in CTL (2010).
From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)






As foreshadowed above — see the subsection titled, Recognizing Statements of Conservation — the angular momentum of a Lagrangian fluid element will be conserved if the "source" term, . This situation will arise if, at the fluid element's location, the azimuthal pressure variation, , and the azimuthal variation in the gravitational potential, , are both zero, or if the two balance one another (i.e.,).
Based on the above discussion, we can equally well view the flow from a frame of reference that is rotating with a constant angular velocity, , and write,



where, as before,



Also, following the earlier discussion, if one wants to follow the timevariation of the fluid's inertialframe angular momentum at a fixed location in inertial space, then the appropriate Eulerian representation of this azimuthal component of the equation of motion is,



If, however, one wants to follow the timevariation of the fluid's inertialframe angular momentum at a fixed location on a rotating coordinate grid, then the appropriate Eulerian representation of this azimuthal component of the equation of motion is obtained by replacing the "transport" velocity, with ; specifically,



An Element of the Hybrid Scheme
This last equation displays one subtle, but valuable, element of the hybrid scheme developed by Call, Tohline, & Lehner (2010). The velocity component, , that appears in the formulation of the relevant conserved quantity — the inertialframe angular momentum density — is drawn from the velocity vector, , which is different from the transport velocity vector, , that defines the Eulerian frame from which the dynamical evolution of the system is being viewed. This equation is usually written, instead, in a form such that the angular momentum density is expressed in terms of the azimuthal component of the transport velocity; see, for example, equation (7) in Norman & Wilson (1978) and equation (12) in New & Tohline (1997). In this more familiar formulation, the momentum density and the transport velocity both directly refer to the same frame of reference. But, as a consequence, the source term is more complicated.
The more familiar formulation — including its modified source term — can be derived from our "hybrid" formulation by recognizing that,



So we can write,



where, as shorthand, we have used,



This implies,









As we see, all terms involving the velocity now explicitly refer to and, hence, to the velocity as measured in the rotating reference frame. But the source now includes a Coriolis term. This corresponds the scalar equation and representation referred to as "Case B ()" in CTL (2010).
From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)






Even Broader Generalization
As Call, Tohline, & Lehner (2010) point out, we are free to measure — and follow the evolution of — the angular momentum density with respect to any of a variety of different rotating frames of reference. Specifically, we are not constrained to choose between the inertial (nonrotating) frame — in which the measured angular momentum density is — and the "grid" frame — in which the measured angular momentum density is . Quite generally, we can choose to measure the angular momentum with respect to a separate "primed" frame that is rotating with angular velocity and in which the measured azimuthal component of the fluid velocity is,
With this definition in hand, we also recognize that,
These two substitutions allow us to rewrite the angular momentum evolution equation in the forms that Call, Tohline, & Lehner (2010) label as Case C () and Case C ().
From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)






From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)






In the latter case, becomes , as it should, when the choice is made to measure the angular momentum density in the "grid" frame, that is, when the choice is made to set .
Zach's Dissertation Paper
Section 2.3
Building on our introductory discussion of the Euler equation (see also Appendix 1.D, §3 of BT87), we begin with the,
Lagrangian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame
,
where we choose to define the frame rotation by the vector,
and,



is the velocity field as viewed from the frame of reference that is rotating at constant angular frequency, . (Note that we can retrieve the inertialframe Euler equation and inertialframe variables by setting at any point in the subsequent derivations.) Because the velocity field introduced by frame rotation is divergence free, that is, because,



the relevant (rotatingframe) continuity equation is identical in form to its inertialframe counterpart, specifically,



We can transform the above rotatingframe Euler equation to a momentum conservation equation by multiplying the equation through by and using the continuity equation to combine and simplify the lefthandside, obtaining,
,
In order to convert this generalpurpose vector equation into the specific set of scalar component equations that embody the desired elements of our hybrid scheme, we need to:
 Step 1: Choose the unitvector basis set associated with the momentum components that we want to track — if tracking linear momentum, or, if tracking radial and angular momentum — and break the vector equation into these specified components;
 Step 2: In all relevant equations, replace the scalar components of the "rotatingframe" momentum density with the scalar components of the "intertialframe" momentum density, drawing each component relation from the vector transformation, ;
 Step 3: Write all of the spatial operators in terms of spatial derivatives that are associated with the unitvector basis set of the desired computational mesh.
Focus on Tracking Linear Momentum
Step 1
If the focus is on tracking linear momentum components, then we need to break the vector momentum equation into its components. This is done by, in turn, "dotting" each unit vector into the vector equation. It is straightforward once we appreciate that the orientation of these Cartesian unit vectors does not vary in space and that, within the context of the rotating frame on which they are defined, these unit vectors do not vary in time. Hence, the first term in the vector equation — the material time derivative — can be written as,
and the process of "dotting" each unit vector into the equation leads to the following set of scalar momentumcomponent equations:
























In deriving these expressions, we also (a) have recognized from the start that, when expressed in Cartesian coordinates,









and (b) have used the familiar operator mapping,
to shift from the total (Lagrangian) time derivative to the partial (Eulerian) time derivative, which is usually the more desirable representation for computational simulations.
Step 2
Next, throughout this set of scalar equations, we replace each component of with the corresponding component of , that is, we perform the following mappings:









As a result, the first of the three "hybrid" momentumcomponent equations becomes,






Referencing the continuity equation, the middle bracketed term on the lefthand side can be set to zero; and the last term on the lefthand side,



can be combined with terms on the righthand side — cutting the Coriolis term in half and canceling the centrifugal acceleration term — to give,



Following a similar sequence of steps, the other two "hybrid" momentum conservation relations become,






Step 3
Cartesian Grid: Assuming the numerical simulation will be conducted on a Cartesian coordinate mesh, the divergence (advection) term on the lefthandside should be evaluated by breaking the transport velocity into its three Cartesian components,
and, on the righthandside, the projection of the spatial operators should be written in the familiar form,
In summary, then, the relevant set of momentum conservation equations is,
Cartesian Components of the InertialFrame Momentum 

 
This is the set of equations that has served as the foundation of the Cartesian simulations reported in Byerly, AdelsteinLelbach, Tohline, & Marcello (2014). 
Cylindrical Grid: If, instead, the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindricalcoordinate components. In concert with this, the divergence (advection) term on the lefthandside should be evaluated by breaking the transport velocity into its three cylindrical components,
Furthermore, recognizing that, when written in cylindrical coordinates, the gradient operator is,
and that the unit vectors in cylindrical coordinates can be related to their Cartesian counterparts via the mappings,
the relevant projections of the gradient operator on the righthandsides of the governing equations should take the form,









In summary, then, the relevant set of momentum conservation equations is,
Cartesian Components of the InertialFrame Momentum 


Focus on Tracking Angular Momentum
Step 1
If the focus is on tracking angular momentum, then we need to break the vector momentum equation into its components. As before, this is done by "dotting" each unit vector into the vector equation. This is less straightforward than in the Cartesian case because the orientation of both the and unit vectors vary in space. As a result, the first term in the vector equation — the material time derivative — generates a couple of extra terms, viz.,









We also recognize that, when expressed in cylindrical coordinates,









Hence, the process of "dotting" each unit vector into the equation leads to the following set of scalar momentumcomponent equations:
















(mult. thru by R) 















ASIDE: If we pause our discussion here and map this set of component equations onto a (rotating) cylindrical coordinate mesh — that is, if on the righthandsides we implement the straightforward operator projections,
we obtain a formulation that is familiar to the astrophysics community. For example, as the following table of equations illustrates, it is the component set that has been spelled out in equations (5)  (7) of Norman & Wilson (1978) and equations (11), (12), & (3) of New & Tohline (1997).
Cylindrical Components of the RotatingFrame Momentum 

 
Note: For complete correspondence, set in all three component equations. 

Note: When comparing this set of equations to the set presented by Norman & Wilson (1978), the definitions of the variables, and , must be swapped. 
[Comment by J. E. Tohline (April 7, 2014)] This is the set of equations that my research group has been using to simulate a wide variety of astrophysical fluid flows over the past twenty years. This is no longer our method of choice, however. A numerical algorithm based on the hybrid scheme, as summarized below, is far preferable to an algorithm that is based on this more familiar, traditional set of equations for several reasons:
 In the hybrid scheme, the Coriolis term disappears from the source term, so it is much easier to design and implement a computational algorithm that conserves angular momentum conservation.
 Although the hybrid scheme advects inertialframe angular momentum, it retains all of the advantages associated with using a rotating frame of reference; for example, numerical diffusion is less severe and, in general, the Courantlimited timestep is larger than would be the case if you were forced to transport fluid in the inertial frame of reference.
 The hybrid scheme facilitates transport (and conservation) of angular momentum across a (rotating) Cartesian mesh. This facilitates the use of adaptivemesh refinement (AMR) techniques and simplifies loadbalancing on distributed memory, highperformance computers.
Step 2
Next, throughout this set of scalar equations, we replace each component of with the corresponding component of , that is, we perform the following mappings:









As a result, the first and third of the three "hybrid" momentumcomponent equations become, respectively,








The second of the three "hybrid" momentum component equations — the one governing conservation of angular momentum — becomes,






Referencing the continuity equation, as before, the middle bracketed term on the lefthand side can be set to zero; and the last term on the lefthand side,



matches and, hence, exactly cancels the Coriolis term on the righthand side to give,




Step 3
Cylindrical Grid: If the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindricalcoordinate components. That is, as in the above "ASIDE," the appropriate operator projections on the righthandsides of the equations are,
In concert with this, the divergence (advection) term on the lefthandside should be evaluated by breaking the transport velocity into its three cylindrical components,
The relevant set of momentum conservation equations is, therefore,
Cylindrical Components of the InertialFrame Momentum 


Cartesian Grid: If, instead, the numerical simulation is to be conducted on a Cartesian coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their Cartesiancoordinate components. In concert with this, the divergence (advection) term on the lefthandside should be evaluated by breaking the transport velocity into its three Cartesian components,
Furthermore, recognizing that, when written in Cartesian coordinates, the gradient operator is,
and that the unit vectors in Cartesian coordinates can be related to their cylindrical counterparts via the mappings,
the relevant projections of the gradient operator on the righthandsides of the governing equations should take the form,









In summary, then, the relevant set of momentum conservation equations is,
Cylindrical Components of the InertialFrame Momentum 

 
This is the set of equations that has served as the foundation of the socalled Hybrid simulations reported in Byerly, AdelsteinLelbach, Tohline, & Marcello (2014). 
Abbreviated Arguments
This is the set of abbreviated arguments that I think should be used in §2.3 of our "hybrid scheme" ApJ paper, that is, the paper by Z. D. Byerly, B. AdelsteinLelbach, J. E. Tohline, & D. C. Marcello (2014).
We evolve the following two coupled fluid equations,



and,



where, is the mass density, is the gas pressure, both and identify the same fluid velocity field (that is, ), and is the gravitational potential generated by a central point mass, , specifically,



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, , in the nonlinear advection term of equation (7). See Appendix A for further elaboration.
[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.]
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 — , , — we start with,









and, when advecting cylindrical momentum components — , , — the component is identical to the Cartesian case but the other two orthogonal components are,






where, .
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 by the velocity field as viewed from a frame of reference that is rotating at angular frequency, , namely,



because , 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 is replaced by — in particular, all five advected quantities, , , , , and , still refer to components of the inertialframe momentum (or angularmomentum) 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:









Resolving Discrepancy Between Zach's Expressions and Mine 


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,
Hence,
This means that the lefthandside of the xmomentum conservation equation becomes,
Similarly (but without detailing the steps), the lefthandside of the ymomentum conservation equation becomes,
So, if rotatingframe quantities are advected on a rotating Cartesian grid, the governing set of equations is,

Second, we note that an evaluation of the advection term that appears on the lefthandside of each component of the momentum equation, which is generically of the form,

requires an assessment of the divergence of the threedimensional 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 timerateofchange 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, , 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 righthandside of each component of the momentum equation is also explicitly broken down into its Cartesian components. This means that, for our hybrid scheme, the righthandsides of equations (nn) and (mm) incorporate the operator projections,






With all of these recognitions in hand, in our hybrid scheme the three components of the cylindrical momentum equations are,









As discussed by Call, Tohline, & Lehner (2010), it is noteworthy that the righthandside of the hybridscheme equation that governs transport (and conservation) of angular momentum does not contain a Coriolis term. This is because , the quantity being advected and tracked, is the angular momentum density as measured in the inertial frame of reference. As is demonstrated in §A.4 of Appendix A, a Coriolis term arises if the equation is written in a form where the quantity being advected is the rotatingframe angular momentum density. This equation, which contains a Coriolis term, is more familiar to the astrophysics community — see, for example, Norman & Wilson (1980) and New & 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 and can be straightforwardly expressed in terms of Cartesian components of or . Specifically, remembering that ,






Text Following Equation A17
I recommend that the text following equation A17 read as follows:
Equation (A17)
[Do not indent] When comparing equations (A16) and (A17), notice that the conserved quantity is the same — the zcomponent of the angular momentum measure in the inertial frame. The only difference in the two equations is the "transport" velocity ( for the nonrotating frame, for the rotating reference frame).
[Start new paragraph] Equation (A17) is different from the more familiar formulation, where the angular momentum density as well as the transport velocity is measure with respect to the rotating frame, i.e., where the angular momentum density is expressed in terms of the azimuthal component of the transport velocity, . But, as a consequence, the source term in the more familiar formulation is more complicated. We can derive the more familiar formulation from equation (A17) by recognizing that,
Equation (A18) thru Equation (A23)
where the last step is accomplished by making use of the continuity relation, Notice that all velocities now refer to , the velocity as measured in the rotating frame, which is the more familiar formulation. The appearance of a Coriolis term is the result of choosing to measure angular momentum in the rotating frame rather than in the inertial frame. This corresponds to "Case B " in Call, Tohline, & Lehner (2010). In our hybrid scheme we have chosen to use equation (A17) instead of (A23) primarily because equation (A17) presents a simpler source term.
TRASH (do not read below this line)
Traditional Eulerian Representation (Review)
Here we review the traditional Eulerian representation of the Euler Equation, as has been discussed in detail earlier.
in terms of velocity:
The socalled "Eulerian form" of the Euler equation can be straightforwardly derived from the standard Lagrangian representation to obtain,
Eulerian Representation
of the Euler Equation,
in terms of momentum density:
Also, we can multiply this expression through by and combine it with the continuity equation to derive what is commonly referred to as the,
Conservative Form
of the Euler Equation,
The second term on the lefthandside of this last expression represents the divergence of the "dyadic product" of the vector momentum density () and the velocity vector and is sometimes written as, .
Component Forms
Let's split the vector Euler equation into its three scalar components; various examples are identified in Table 1.
Example # 
Grid 
Momentum Vector 


Basis 
Rotating? 
Basis 
Frame 

1 
Cartesian 
No 
Cartesian 
Inertial 
2 
Cylindrical 
Yes 
Cylindrical 
Rotating 
3 
Cylindrical 
Yes 
Cylindrical 
Rotating 
In the following expressions, we will use to denote the fluid velocity when it is associated with the rate of fluid transport across the coordinate grid, and we will use to denote the fluid velocity when it is associated with the momentum density that is being advected. In all cases, it should be understood that , as both vectors refer to the same fluid velocity. In addition, we will use a "prime" notation to indicate when a velocity is being viewed from a rotating frame of reference; specifically, we will consider rotation about the axis of the coordinate system, that is,



and,



but we will not insist that the two rotation frequencies, and , have the same value. Hence, in general, . It is worth emphasizing that, because we will only be considering frame rotation about the axis, the cylindrical and components of the velocity are interchangeable, that is: ; and .
Example #1
This is certainly the most familiar component set.









where, for any one of the three scalar PDEs, advection of the relevant component of the momentum density, , is handled via the operation,



Example #2
This component set has been spelled out in, for example, equations (5)  (7) of Norman & Wilson (1978) and equations (11), (12), & (3) of New & Tohline (1997).












where, as noted above,



Example #3
















where, as noted above,


