Skip to main content
Sandbox Physics

L2 · Derivation and worked examples

Why a vertical mesh changes a linear channel

Sandbox Physics ·

A vertical mesh replaces a spatial derivative with a difference between neighbouring samples. That replacement slightly changes the magnetic coupling of a wave, so a real grid calculation need not grow at exactly the continuum rate. Refining the mesh should reduce the difference in a predictable way. The saved channel calculations let us check that prediction without generating a new trajectory.

Begin with the four real channel equations and the dispersion relation developed in the preceding chapters. You will need sine and cosine, derivatives as local slopes, and a sum as an average over samples. This chapter derives the grid's effective wave number, reads an actual spatial refinement, and explains why energy and stress should still be calculated from the saved fields.

Keep the fields as the evolving objects

The local coordinates are radial outward xx, orbital yy, and vertical zz. The positive rotation rate is Ω\Omega, and the background orbital velocity is −qΩx-q\Omega x. The density ρ\rho and vertical magnetic field are uniform. Write magnetic disturbances in velocity units, b=δB/μ0ρb=\delta B/\sqrt{\mu_0\rho}, with μ0\mu_0 the magnetic permeability. Both bb and the horizontal fluid velocity have units of length per time.

Distinguish a spatial field from its amplitude. Let Ux(z,t)U_x(z,t) and Uy(z,t)U_y(z,t) be the actual sampled horizontal velocity perturbations, and let Bx(z,t)B_x(z,t) and By(z,t)B_y(z,t) denote the horizontal magnetic disturbances in those velocity units. The upper-case magnetic symbols here do not denote SI magnetic-field strengths. The total orbital velocity in this rotating frame is

vy=−qΩx+Uy(z,t).v_y=-q\Omega x+U_y(z,t).

The local linear field equations are

∂tUx=2ΩUy+vA∂zBx,∂tUy=−(2−q)ΩUx+vA∂zBy,∂tBx=vA∂zUx,∂tBy=−qΩBx+vA∂zUy.\begin{aligned} \partial_t U_x&=2\Omega U_y+v_A\partial_z B_x,\\ \partial_t U_y&=-(2-q)\Omega U_x+v_A\partial_z B_y,\\ \partial_t B_x&=v_A\partial_z U_x,\\ \partial_t B_y&=-q\Omega B_x+v_A\partial_z U_y. \end{aligned}

Here vAv_A is the constant vertical Alfvén speed. Substituting Ux=ux(t)sin⁡(kz)U_x=u_x(t)\sin(kz) and Bx=bx(t)cos⁡(kz)B_x=b_x(t)\cos(kz), and similarly for the orbital components, recovers the four-amplitude equations. The magnetic force has a minus sign in the amplitude equation because the derivative of cosine is minus sine.

The independent offline calculation advances four arrays of sampled fields. Its spatial derivative reads neighbouring array values. A sine projection is made only when saving an observation; the projected amplitudes never replace the arrays or feed the next step. A single-mode experiment can have a four-dimensional description and still test a real-space derivative. These are two ways of representing the same restricted linear problem, with different numerical errors to examine.

A centred difference has its own wave number

Divide the periodic vertical length LL into NN equally spaced samples. Their spacing is h=L/Nh=L/N, with zj=jhz_j=jh and periodic indices: the neighbour after the last sample is the first one. A centred difference uses the slope across both neighbours,

(Dhf)j=fj+1−fj−12h.(D_h f)_j=\frac{f_{j+1}-f_{j-1}}{2h}.

For a sinusoid, write the local phase as θj=kzj\theta_j=kz_j and the one-step phase change as η=kh\eta=kh. Use the addition formulas before taking any small-spacing approximation:

sin⁡(θj+η)−sin⁡(θj−η)=2cos⁡(θj)sin⁡(η),Dhsin⁡(kzj)=khcos⁡(kzj),Dhcos⁡(kzj)=−khsin⁡(kzj),kh=sin⁡(kh)h.\begin{gathered} \sin(\theta_j+\eta)-\sin(\theta_j-\eta)\\ =2\cos(\theta_j)\sin(\eta),\\ D_h\sin(kz_j)=k_h\cos(kz_j),\\ D_h\cos(kz_j)=-k_h\sin(kz_j),\\ k_h=\frac{\sin(kh)}{h}. \end{gathered}

The quantity khk_h, called the modified or effective wave number, has the same inverse-length unit as kk. The grid differentiates this sampled wave exactly as though its derivative coupling used khk_h. For a well-resolved fundamental wave, 0<kh<π0<kh<\pi, it is smaller than kk. Expanding sine gives

khk=sin⁡(kh)kh=1−(kh)26+O((kh)4).\begin{aligned} \frac{k_h}{k}&=\frac{\sin(kh)}{kh}\\ &=1-\frac{(kh)^2}{6}+O((kh)^4). \end{aligned}

The leading derivative error therefore falls by about four when the spacing halves. For this mode, the spatially differenced equations, with time still continuous, use Kh=vAkhK_h=v_Ak_h in place of the continuum coupling K=vAkK=v_Ak. This statement describes the derivative. It does not assert that every possible diagnostic has exactly the same leading error.

Analytic modified-wave-number deficit and actual saved projected-vector errors on three spatial meshes

The upper panel is an analytic reference: the fractional wave-number deficit for the fundamental periodic wave. The lower panel uses the actual saved four-component projection errors at the same final time. The dashed lower curve is a second-order reference anchored to the coarsest saved result. The two panels measure different quantities. Both axes showing error are dimensionless; the horizontal coordinate is the number of vertical samples.

Worked example: the same seed on three meshes

Use box length L0L_0, time unit T0T_0, and the saved parameters q=3/2q=3/2, Ω=1/T0\Omega=1/T_0, K=0.6/T0K=0.6/T_0. The fundamental wave has k=2π/L0k=2\pi/L_0, so the background speed is vA=0.6L0/(2πT0)v_A=0.6L_0/(2\pi T_0). It is not 0.6L0/T00.6L_0/T_0.

At N=64N=64, the angle between neighbouring samples is 2π/642\pi/64. Thus

khk≃0.998394393,KhT0≃0.599036636.\begin{aligned} \frac{k_h}{k}&\simeq0.998394393,\\ K_hT_0&\simeq0.599036636. \end{aligned}

The continuum growing rate is sT0=0.44≃0.663324958sT_0=\sqrt{0.44}\simeq0.663324958. The same continuum initial vector is used on all three meshes. With radial amplitude A=10−4L0/T0A=10^{-4}L_0/T_0, its four amplitudes are

a(0)A=(10.4/0.440.6/0.44−1.5).\frac{\mathbf a(0)}{A} =\begin{pmatrix} 1\\ 0.4/\sqrt{0.44}\\ 0.6/\sqrt{0.44}\\ -1.5 \end{pmatrix}.

To read the saved radial amplitude, multiply every actual velocity sample by the known sine pattern and average:

ax(t)=2N∑j=0N−1Ux,j(t) sj,sj=sin⁡ ⁣(2πjN).\begin{gathered} a_x(t)=\frac{2}{N}\sum_{j=0}^{N-1}U_{x,j}(t)\,s_j,\\ s_j=\sin\!\left(\frac{2\pi j}{N}\right). \end{gathered}

The factor two compensates for the mean squared sine, which is one-half on these meshes. The other velocity uses the same sine, and the two magnetic components use cosine. Compare this four-component observation aN\mathbf a_N with the continuum vector ac\mathbf a_{\rm c} at the same time:

eN=∥aN−ac∥2∥ac∥2.e_N=\frac{\|\mathbf a_N-\mathbf a_{\rm c}\|_2} {\|\mathbf a_{\rm c}\|_2}.

The norm squares the four differences, adds them and takes a square root. This is the saved modal-vector error; it is not a separately measured maximum error at every grid point. At t=3T0t=3T_0, with the same time step 0.0025T00.0025T_0, the actual results are

Vertical samplesProjected vector errorMeasured radial rate times the time unit
640.001287717890.001287717890.6629361790.662936179
1280.0003219424740.0003219424740.6632278860.663227886
2560.00008048641820.00008048641820.6633006980.663300698

The error ratios are 3.9998383.999838 and 3.9999603.999960. They support second-order spatial convergence for this experiment. The measured rate is s^N=t−1ln⁡(∣ax(t)∣/A)\widehat s_N=t^{-1}\ln(|a_x(t)|/A). Inserting KhK_h into the quartic gives shT0≃0.662855038s_hT_0\simeq0.662855038 on the 64-sample mesh, close to but distinct from the measured finite-time rate. The original seed is a continuum eigenvector, not an exact eigenvector of the grid equations. Its small additional mode components help explain that distinction.

For a prediction at a new resolution, suppose the same second-order range continues. Halving the finest spacing suggests an error near 2.01×10−52.01\times10^{-5}. This fourth resolution is a prediction for a future comparison.

Sum the actual fields when measuring transport

For arbitrary sampled fields, the vertically averaged perturbation energy and the two outward momentum fluxes are

Eh=Eu,h+Eb,h,Eu,h=ρ2N∑j(Ux,j2+Uy,j2),Eb,h=ρ2N∑j(Bx,j2+By,j2),TRe,h=ρN∑jUx,jUy,j,TMax,h=−ρN∑jBx,jBy,j.\begin{aligned} E_h&=E_{u,h}+E_{b,h},\\ E_{u,h}&=\frac{\rho}{2N}\sum_j\left(U_{x,j}^2+U_{y,j}^2\right),\\ E_{b,h}&=\frac{\rho}{2N}\sum_j\left(B_{x,j}^2+B_{y,j}^2\right),\\ T_{{\rm Re},h}&=\frac{\rho}{N}\sum_j U_{x,j}U_{y,j},\\ T_{{\rm Max},h}&=-\frac{\rho}{N}\sum_j B_{x,j}B_{y,j}. \end{aligned}

The two energy terms are the kinetic and magnetic pieces. Energy is energy per volume and each stress is force per area. When the fields are exactly the stated sine/cosine mode, averaging the squared patterns recovers the amplitude formulas with factors one-quarter and one-half. For a field containing other modes, products of four projected amplitudes need not equal these full-array sums. The saved diagnostics use the actual arrays.

Periodic centred differences also obey a useful cancellation. Shifting the summation index gives ∑jfj(Dhg)j=−∑jgj(Dhf)j\sum_j f_j(D_hg)_j=-\sum_j g_j(D_hf)_j. Magnetic tension therefore exchanges kinetic and magnetic energy without creating their sum in the continuous-time grid equations. The remaining source is qΩ(TRe,h+TMax,h)q\Omega(T_{{\rm Re},h}+T_{{\rm Max},h}). The next guide examines how a finite time step approximates that work.

Exercises with explanations

1. A wave alternates sign at every sample. Does the centred difference measure its steepest possible slope? No. For the alternating pattern, the two neighbours have the same value, so the difference is zero. This is the mesh's shortest represented wavelength. Increasing the wave number toward that limit does not keep increasing the effective coupling. A resolved fundamental mode and a near-grid-scale disturbance are different tests.

2. A calculation near the fastest continuum mode shows an unusually small growth-rate error. Does that prove fourth-order spatial accuracy? No. The growth curve is stationary at its maximum: its first derivative with respect to coupling vanishes there. A coupling error of second order can then produce a fourth-order change in the ideal growing eigenvalue. Finite-time projections can also contain other components. Check a nonstationary coupling and a field-vector error; the saved spatial series deliberately uses K=0.6/T0K=0.6/T_0.

3. Why use the same continuum seed on all three spatial meshes? It exposes the effect of changing the spatial derivative while holding the physical initial disturbance fixed. Seeding each mesh with its own grid eigenvector answers a different useful question, used to isolate temporal error in the next guide.

4. Could a small projected error hide a grid-scale disturbance? Yes. A component orthogonal to the fundamental sine/cosine pattern can disappear from these four projections. Inspect actual fields and their energy/stress sums as well. Convergence here qualifies the observed fundamental; other disturbances require their own tests.

Sources

Balbus and Hawley, A Powerful Local Shear Instability in Weakly Magnetized Disks. I. Linear Analysis, ApJ 376, 214–222 (1991), printed pages 216–218, equations 2.2 and 2.5–2.9, supply the continuum component equations and dispersion relation. Primary scan. Stone and Gardiner, Implementation of the Shearing Box Approximation in Athena, arXiv:1006.0139v1, Sections II–III, define the local shear and discuss numerical source terms. Pinned primary source.

The spatial tables and field observations come from the eleven-case linear-channel-replay-public.zip, specifically records/spatial-64, records/spatial-128 and records/spatial-256. Athena uses separate finite-volume and constrained-transport algorithms. The centred-difference derivation here describes the bundle's simple periodic channel calculation. Reading these existing arrays adds zero new evolution steps.

Try the prediction

Compare the saved spatial meshes

  1. Predict how halving the vertical spacing changes the projected vector error.
  2. Open the offline replay and compare the three saved spatial meshes at the same final time.
  3. Read actual field profiles, the four-component errors and successive error ratios.
  4. Use the centered derivative and its modified wave number to explain the convergence.

Keep reading

What this page links to directly

Magnetic tension and the growth of a disk disturbance · Reading path

Magnetic tension and the growth of a disk disturbance

  1. ArticleWho carries angular momentum through an accretion disk?
  2. ArticleFour real fields in a rotating patch of a disk
  3. ArticlePredict the growth, then follow its energy
  4. You are hereWhy a vertical mesh changes a linear channel
  5. ArticleChecking growth, transport and shear work
  6. ArticleCheck an MRI calculation from its saved fields
  7. ArticleRead beyond linear MRI growth
Something broken, confusing, or missing on this page?Tell us
What kind of feedback?

Sent with this page’s address, its Lab ID, the site language and your device type. Please leave out personal details. Feedback is kept for up to 12 months.Privacy policy