Hybrid Advection Scheme

Preface

March 1, 2014 by Joel E. Tohline

<Click here for full preface> … Throughout my research career, I have sought new or modified techniques and algorithms that would allow my group to perform more accurate numerical simulations of astrophysical fluid flows. At the beginning — following the advice of my dissertation co-advisors, Peter Bodenheimer and David Black — I adopted a cylindrical, rather than Cartesian, computational grid. When a cylindrical coordinate system is used, one of the components of the equation of motion can be written in a form that fairly naturally conserves angular momentum, and this was quite a desirable feature, given that our investigations were focusing on analyzing the stability of rotating configurations.

During the mid-90s, when Kimberly Barker New was conducting her dissertation research, we started using a rotating cylindrical coordinate mesh. On a grid that was spinning at a suitably chosen frequency, advection of fluid through the grid could be minimized and this, in turn, reduced the undesirable effects of numerical diffusion. We adopted a fairly standard algorithmic approach very similar to the one used by Norman & Wilson (1980) and acknowledged that we were making a tradeoff: While the shift to a rotating cylindrical coordinate mesh reduced the effects of numerical diffusion, the shift introduced a rather ugly Coriolis "source" term into two components of the equation of motion. This made it more difficult to ensure conservation of angular momentum.

During my continuing efforts to develop an improved computational fluid dynamics algorithm, one day I noticed that, while imposing some fairly reasonable constraints, the Coriolis term could be removed from the "source" term and folded into the divergence term on the left-hand side of the angular momentum conservation equation. This manipulation of terms seemed to be saying that the undesirable Coriolis term would disappear while employing a rotating coordinate mesh if the variable that was advected through the grid was the inertial-frame angular momentum density, rather than the rotating-frame angular momentum density. 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. My notes on this topic have been preserved, as they were included in my earliest version of this web-based H_Book; the relevant page can be accessed here, which is an html file whose linux time stamp is August 27, 2000. The "symbol" fonts utilized throughout this old html page seem now to be readable only through Microsoft's Internet Explorer web browser. Hence, for posterity sake, I have retyped this "year 2000" set of notes into an accompanying page of this wiki.

Setting the Stage

Recognizing Statements of Conservation

When dealing with the compressible fluid equations, we will often encounter hyperbolic PDEs of the following form:

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{v} </math>

where we are using <math>~\vec{v}</math> 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 time-rate of change of <math>~\psi</math> is being measured in a so-called 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, <math>~S</math>, on the right-hand side is zero, we will be able to identify the scalar variable, <math>~\psi</math>, as the volume density of some conserved quantity. For example, the continuity equation — which is a mathematical statement of mass conservation — has the form,

<math>\frac{d\rho}{dt} + \rho \nabla \cdot \vec{v} = 0</math>

or, equivalently,

<math> ~- \nabla\cdot \vec{v} \, , </math>

where, <math>~\rho</math> 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,

<math> \frac{ds}{dt} + s\nabla\cdot \vec{v} </math>

or, equivalently,

<math> ~- \nabla\cdot \vec{v} \, , </math>

where, <math>~s</math> 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,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{v} </math>

where, <math>~\varpi</math> is the Lagrangian fluid element's (cylindrical radial) distance measured from the symmetry axis of the underlying potential and <math>~v_\phi = \varpi\dot\phi</math> is the azimuthal component of the inertial velocity field, <math>~\vec{v}</math>, at the location of the fluid element.

Alternative Reference Frames

Now, we might want to examine the time-dependent 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 <math>~\vec{v}</math> with,

<math> \vec{v} - \vec{v}_\mathrm{frame} \, . </math>

(Here, we will consider only time-independent functional expressions for the frame velocity, <math>~\vec{v}_\mathrm{frame}</math>.) 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 <math>~\vec{v}_\mathrm{frame}</math> is a divergence-free velocity field, then the transformation is trivial. For example, if we choose a frame of reference that is rotating uniformly with angular velocity, <math>~\Omega_0</math>, then,

<math> \vec{v}_\mathrm{frame} </math>

<math> \boldsymbol{\hat{e}}_\phi (\varpi \Omega_0) \, , </math>

and, utilizing cylindrical coordinates,

<math> \nabla\cdot\vec{v}_\mathrm{frame} </math>

<math> \frac{\partial(0)}{\partial \varpi} + \frac{1}{\varpi}\frac{\partial(\varpi \Omega_0)}{\partial \phi} + \frac{\partial(0)}{\partial z} = 0 \, . </math>

Hence,

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{u} </math>

<math> \frac{d\psi}{dt} + \psi \nabla\cdot [\vec{v} - \vec{v}_\mathrm{frame}] </math>

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{v} \, , </math>

so the new generic hyperbolic PDE becomes,

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{u} </math>

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 time-rate 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:

<math>~~~\rightarrow~~~</math>

<math> \frac{\partial \psi}{\partial t} + \vec{u} \cdot \nabla\psi \, . </math>

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

<math> \frac{\partial\psi}{\partial t} + \vec{u} \cdot \nabla\psi + \psi \nabla\cdot \vec{u} </math>

or, more succinctly,

<math> \frac{\partial\psi}{\partial t} + \nabla\cdot (\psi \vec{u} ) </math>

This equation also is broadly recognized as a conservation statement because, when <math>~S = 0</math>, the variable <math>~\psi</math> will represent the volume density of a conserved quantity.

We should emphasize that the inertial-frame version of this Eulerian conservation equation can be retrieved straightforwardly by setting <math>~\Omega_0 = 0</math>, which is equivalent to setting <math>~\vec{u} = \vec{v}</math>. It is,

<math> \frac{\partial\psi}{\partial t} + \nabla\cdot (\psi \vec{v}) </math>

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 <math>~\vec{v}</math>) provides the time-rate of change of <math>~\psi</math> at a fixed point in inertial space, while the "new" version (explicitly containing <math>~\vec{u}</math>) provides the time-rate of change of <math>~\psi</math> 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,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{v} </math>

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, . </math>

For this equation, the source term is identified as,

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, , </math>

and <math>~\psi = (\rho\varpi v_\phi)</math> is the inertial-frame angular momentum density, as measured with respect to the <math>~z</math>-coordinate axis. This corresponds the scalar equation and representation referred to as "Case B (<math>~\eta=3</math>)" in CTL (2010).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010) Case B <math>~(\eta = 3)</math> with the following replacements: <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math> ; <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math> ; <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>
<math>~\psi_{(3)}</math>	<math>~S_{(3)}</math>
<math>~\rho \varpi v_\phi</math>	<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi}</math>

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, <math>~S = 0</math>. This situation will arise if, at the fluid element's location, the azimuthal pressure variation, <math>~\partial P/\partial\phi</math>, and the azimuthal variation in the gravitational potential, <math>~\partial \Phi/\partial\phi</math>, are both zero, or if the two balance one another (i.e.,<math>~\partial P/\partial\phi=-\rho\partial\Phi/\partial\phi</math>).

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, <math>~\Omega_0</math>, and write,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{u} </math>

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, , </math>

where, as before,

<math>~\equiv~</math>

<math> \vec{v} - \boldsymbol{\hat{e}}_\phi \varpi\Omega_0 \, . </math>

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

<math> \frac{\partial (\rho \varpi v_\phi)}{\partial t} + \nabla\cdot [(\rho \varpi v_\phi) \vec{v}] </math>

If, however, one wants to follow the time-variation of the fluid's inertial-frame 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, <math>~\vec{v}</math> with <math>~\vec{u}</math>; specifically,

<math> \frac{\partial (\rho \varpi v_\phi)}{\partial t} + \nabla\cdot [(\rho \varpi v_\phi) \vec{u}] </math>

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, <math>~v_\phi</math>, that appears in the formulation of the relevant conserved quantity — the inertial-frame angular momentum density — is drawn from the velocity vector, <math>~\vec{v}</math>, which is different from the transport velocity vector, <math>~\vec{u}</math>, 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,

<math> ~u_\phi + \varpi\Omega_0 \, . </math>

So we can write,

<math> \frac{\partial [\rho \varpi (u_\phi + \varpi\Omega_0 ) ]}{\partial t} + \nabla\cdot \{[\rho \varpi ( u_\phi + \varpi\Omega_0)] \vec{u} \} </math>

where, as shorthand, we have used,

<math>~\equiv~</math>

<math> - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} \, . </math>

This implies,

<math> \frac{\partial (\rho \varpi u_\phi )}{\partial t} + \nabla\cdot [ (\rho \varpi u_\phi) \vec{u} ] </math>	<math>~=~</math>	<math> S_{\phi i} - \frac{\partial [\rho \varpi (\varpi\Omega_0 ) ]}{\partial t} - \nabla\cdot \{[\rho \varpi (\varpi\Omega_0)] \vec{u} \} </math>
	<math>~=~</math>	<math> S_{\phi i} - \varpi^2\Omega_0 \biggl\{ \frac{\partial \rho}{\partial t} + \nabla\cdot (\rho \vec{u} ) \biggr\} - \rho \vec{u}\cdot \nabla(\varpi^2 \Omega_0) </math>
	<math>~=~</math>	<math> S_{\phi i} - 2\rho \varpi u_\varpi \Omega_0 \, . </math>

As we see, all terms involving the velocity now explicitly refer to <math>~\vec{u}</math> 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 (<math>~\eta=3'</math>)" in CTL (2010).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010) Case B <math>~(\eta = 3')</math> as before: <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math> ; <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math> ; <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math> additional replacements: <math>~(\bar\omega u^{t'})_\mathrm{CTL} \rightarrow \Omega_0</math> ; <math>~u^R \rightarrow v_\varpi = u_\varpi</math>
<math>~\psi_{(3')}</math>	<math>~S_{(3')}</math>
<math>~\rho \varpi (v_\phi - \varpi\Omega_0) = \rho \varpi u_\phi </math>	<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} - 2\rho\varpi u_\varpi \Omega_0</math>

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 <math>~(\rho\varpi v_\phi)</math> — and the "grid" frame — in which the measured angular momentum density is <math>~\rho\varpi u_\phi = \rho\varpi(v_\phi - \varpi\Omega_0)</math>. Quite generally, we can choose to measure the angular momentum with respect to a separate "primed" frame that is rotating with angular velocity <math>~\omega_0</math> and in which the measured azimuthal component of the fluid velocity is,

<math> ~v'_\phi ~= ~v_\phi - \varpi \omega_0 \, . </math>

With this definition in hand, we also recognize that,

<math> ~u_\phi = v_\phi - \varpi\Omega_0 = v'_\phi + \varpi (\omega_0 - \Omega_0) \, . </math>

These two substitutions allow us to rewrite the angular momentum evolution equation in the forms that Call, Tohline, & Lehner (2010) label as Case C (<math>~\eta=3</math>) and Case C (<math>~\eta=3'</math>).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010) Case C <math>~(\eta = 3)</math> with the following replacements: <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math> ; <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math> ; <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>
<math>~\psi_{(3)}</math>	<math>~S_{(3)}</math>
<math>~\rho \varpi (v'_\phi +\varpi\omega_0) = \rho \varpi v_\phi</math>	<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi}</math>

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010) Case C <math>~(\eta = 3')</math> as before: <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math> ; <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math> ; <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math> additional replacements: <math>~(\bar\omega u^{t'})_\mathrm{CTL} \rightarrow \Omega_0</math> ; <math>~u^R \rightarrow v_\varpi = u_\varpi</math>
<math>~\psi_{(3')}</math>	<math>~S_{(3')}</math>
<math>~\rho \varpi[ v'_\phi + \varpi (\omega_0 - \Omega_0)] = \rho \varpi u_\phi </math>	<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} - 2\rho\varpi u_\varpi \Omega_0</math>

In the latter case, <math>~v'_\phi</math> becomes <math>~u_\phi</math>, 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 <math>~\omega_0 = \Omega_0</math>.

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 so-called "Eulerian form" of the Euler equation can be straightforwardly derived from the standard Lagrangian representation to obtain,

Eulerian Representation
of the Euler Equation,

<math>~\frac{\partial\vec{v}}{\partial t} + (\vec{v}\cdot \nabla) \vec{v}= - \frac{1}{\rho} \nabla P - \nabla \Phi</math>

in terms of momentum density:

Also, we can multiply this expression through by <math>~\rho</math> and combine it with the continuity equation to derive what is commonly referred to as the,

Conservative Form
of the Euler Equation,

<math>~\frac{\partial(\rho\vec{v})}{\partial t} + \nabla\cdot [(\rho\vec{v})\vec{v}]= - \nabla P - \rho \nabla \Phi</math>

The second term on the left-hand-side of this last expression represents the divergence of the "dyadic product" of the vector momentum density (<math>~\rho</math><math>~\vec{v}</math>) and the velocity vector <math>~\vec{v}</math> and is sometimes written as, <math>\nabla\cdot [(\rho \vec{v}) \otimes \vec{v}]</math>.

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
Example #	Basis	Rotating?	Basis	Frame
1	Cartesian	No	Cartesian	Inertial
2	Cylindrical	Yes <math>~(\Omega_0)</math>	Cylindrical	Rotating <math>~(\Omega_0)</math>
3	Cylindrical	Yes <math>~(\Omega_0)</math>	Cylindrical	Rotating <math>~(\omega_0)</math>

In the following expressions, we will use <math>~\vec{v}</math> to denote the fluid velocity when it is associated with the rate of fluid transport across the coordinate grid, and we will use <math>~\vec{u}</math> 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 <math>~\vec{v} = \vec{u}</math>, 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 <math>~z</math>-axis of the coordinate system, that is,

<math>~v_\phi - R\Omega_0 \, ,</math>

and,

<math>~u_\phi - R\omega_0 \, ,</math>

but we will not insist that the two rotation frequencies, <math>~\Omega_0</math> and <math>~\omega_0</math>, have the same value. Hence, in general, <math>~(\vec{u})' \ne (\vec{v})'</math>. It is worth emphasizing that, because we will only be considering frame rotation about the <math>z</math>-axis, the cylindrical <math>R</math> and <math>z</math> components of the velocity are interchangeable, that is: <math>~u'_R = v'_R = u_R = v_R</math>; and <math>~u'_z = v'_z = u_z = v_z</math>.

Example #1

This is certainly the most familiar component set.

<math>\boldsymbol{\hat{e}}_x: ~~~\frac{\partial (\rho v_x)}{\partial t} + \nabla\cdot[(\rho v_x) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial x} - \rho \frac{\partial \Phi}{\partial x} \, , </math>
<math>\boldsymbol{\hat{e}}_y: ~~~\frac{\partial (\rho v_y)}{\partial t} + \nabla\cdot[(\rho v_y) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial y} - \rho \frac{\partial \Phi}{\partial y} \, , </math>
<math>\boldsymbol{\hat{e}}_z: ~~~\frac{\partial (\rho v_z)}{\partial t} + \nabla\cdot[(\rho v_z) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial z} - \rho \frac{\partial \Phi}{\partial z} \, , </math>

where, for any one of the three scalar PDEs, advection of the relevant component of the momentum density, <math>~\psi_i</math>, is handled via the operation,

<math> \nabla\cdot[\psi_{i} \vec{v} ] </math>

<math> \frac{\partial (\psi_i v_x)}{\partial x} + \frac{\partial (\psi_i v_y)}{\partial y} + \frac{\partial (\psi_i v_z)}{\partial z} \, . </math>

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).

<math>\boldsymbol{\hat{e}}_R: ~~~~~~~\frac{\partial (\rho v_R)}{\partial t} + \nabla\cdot[(\rho v_R) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} + \frac{(\rho R v_\phi)^2}{\rho R^3} + \rho\Omega_0^2 R + \frac{2\Omega_0 (\rho R v_\phi)}{R} \, , </math>
	<math>~=~</math>	<math> -~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} + \frac{\rho}{R} (v_\phi + R\Omega_0)^2 \, , </math>
<math>\boldsymbol{\hat{e}}_\phi: ~~~\frac{\partial (\rho R v_\phi)}{\partial t} + \nabla\cdot[(\rho R v_\phi) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial \phi} - \rho \frac{\partial \Phi}{\partial \phi} - 2\rho (\Omega_0 R )v_R \, , </math>
<math>\boldsymbol{\hat{e}}_z: ~~~~~~~~\frac{\partial (\rho v_z)}{\partial t} + \nabla\cdot[(\rho v_z) \vec{v}~]</math>	<math>~=~</math>	<math> -~\frac{\partial P}{\partial z} - \rho \frac{\partial \Phi}{\partial z} \, , </math>

where, as noted above,