-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathnumerics.html
More file actions
1640 lines (1627 loc) · 129 KB
/
numerics.html
File metadata and controls
1640 lines (1627 loc) · 129 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html lang="en" data-content_root="./">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Numerical technique — Magic 6.3 documentation</title>
<link rel="stylesheet" type="text/css" href="_static/pygments.css?v=03e43079" />
<link rel="stylesheet" type="text/css" href="_static/magic.css?v=a97dc138" />
<script src="_static/documentation_options.js?v=4b074805"></script>
<script src="_static/doctools.js?v=9bcbadda"></script>
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
<script async="async" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
<link rel="search" type="application/opensearchdescription+xml"
title="Search within Magic 6.3 documentation"
href="_static/opensearch.xml"/>
<link rel="icon" href="_static/favicon.ico"/>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="Contributing to the code" href="contribute.html" />
<link rel="prev" title="Formulation of the (magneto)-hydrodynamics problem" href="equations.html" />
<link href='http://fonts.googleapis.com/css?family=Open+Sans:300,400,700'
rel='stylesheet' type='text/css' />
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1/jquery.js"></script>
<script src="galleria/galleria-1.4.2.min.js"></script>
<style type="text/css">
table.right { float: right; margin-left: 20px; }
table.right td { border: 1px solid #ccc; }
</style>
<script type="text/javascript">
// intelligent scrolling of the sidebar content
$(window).scroll(function() {
var sb = $('.sphinxsidebarwrapper');
var win = $(window);
var sbh = sb.height();
var offset = $('.sphinxsidebar').position()['top'];
var wintop = win.scrollTop();
var winbot = wintop + win.innerHeight();
var curtop = sb.position()['top'];
var curbot = curtop + sbh;
// does sidebar fit in window?
if (sbh < win.innerHeight()) {
// yes: easy case -- always keep at the top
sb.css('top', $u.min([$u.max([0, wintop - offset - 10]),
$(document).height() - sbh - 200]));
} else {
// no: only scroll if top/bottom edge of sidebar is at
// top/bottom edge of window
if (curtop > wintop && curbot > winbot) {
sb.css('top', $u.max([wintop - offset - 10, 0]));
} else if (curtop < wintop && curbot < winbot) {
sb.css('top', $u.min([winbot - sbh - offset - 20,
$(document).height() - sbh - 200]));
}
}
});
</script>
</head><body>
<div class="pageheader">
<ul>
<li><a href="index.html">Home</a></li>
<li><a href="install.html">Get it/Run it</a></li>
<li><a href="contribute.html">Contribute!</a></li>
<li><a href="#">Numerical methods</a></li>
<li><a href="contents.html">Contents</a></li>
</ul>
<div>
<a href="index.html">
<img src="_static/logo.png" alt="magic" height="120px" width="192px"/>
</a>
</div>
</div>
<div class="related" role="navigation" aria-label="Related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="f-modindex.html" title="Fortran Module Index"
>fortran modules</a> |</li>
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="contribute.html" title="Contributing to the code"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="equations.html" title="Formulation of the (magneto)-hydrodynamics problem"
accesskey="P">previous</a> |</li>
<li class="nav-item nav-item-0"><a href="contents.html">Magic 6.3 documentation</a> »</li>
<li class="nav-item nav-item-this"><a href="">Numerical technique</a></li>
</ul>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body" role="main">
<section id="numerical-technique">
<span id="secnumerics"></span><h1>Numerical technique<a class="headerlink" href="#numerical-technique" title="Link to this heading">¶</a></h1>
<p><strong>MagIC</strong> is a pseudo-spectral MHD code. This numerical technique was
originally developed by P. Gilman and G. Glatzmaier for the spherical geometry.
In this approach the unknowns are expanded into complete sets of functions in
radial and angular directions: <strong>Chebyshev polynomials or Finite differences in
the radial direction</strong> and <strong>spherical harmonic functions in the azimuthal and
latitudinal directions</strong>. This allows to express all partial derivatives
analytically. Employing orthogonality relations of spherical harmonic
functions and using collocation in radius then lead to algebraic equations that
are integrated in time with a <strong>mixed implicit/explicit time stepping scheme</strong>.
The nonlinear terms and the Coriolis force are evaluated in the physical (or
grid) space rather than in spectral space. Although this approach requires
costly numerical transformations between the two representations (from spatial
to spectral using Legendre and Fourier transforms), the resulting decoupling of
all spherical harmonic modes leads to a net gain in computational speed.
Before explaining these methods in more detail, we introduce the
poloidal/toroidal decomposition.</p>
<section id="poloidal-toroidal-decomposition">
<h2>Poloidal/toroidal decomposition<a class="headerlink" href="#poloidal-toroidal-decomposition" title="Link to this heading">¶</a></h2>
<p>Any vector <span class="math notranslate nohighlight">\(\vec{v}\)</span> that fulfills <span class="math notranslate nohighlight">\(\vec{\nabla}\cdot\vec{v}=0\)</span>, i.e.
a so-called <em>solenoidal field</em>,
can be decomposed in a poloidal and a toroidal part <span class="math notranslate nohighlight">\(W\)</span> and <span class="math notranslate nohighlight">\(Z\)</span>,
respectively</p>
<div class="math notranslate nohighlight">
\[\vec{v} = \vec{\nabla}\times\left(\vec{\nabla}\times W\,\vec{e_r}\right) +
\vec{\nabla}\times Z\,\vec{e_r}.\]</div>
<p>Three unknown vector components are thus replaced by two scalar fields,
the poloidal potential <span class="math notranslate nohighlight">\(W\)</span> and the toroidal potential <span class="math notranslate nohighlight">\(Z\)</span>.
This decomposition is unique, aside from an arbitrary radial function <span class="math notranslate nohighlight">\(f(r)\)</span>
that can be added to <span class="math notranslate nohighlight">\(W\)</span> or <span class="math notranslate nohighlight">\(Z\)</span> without affecting <span class="math notranslate nohighlight">\(\vec{v}\)</span>.</p>
<p>In the anelastic approximation, such a decomposition can be used for the
mass flux <span class="math notranslate nohighlight">\(\tilde{\rho}\vec{u}\)</span> and the magnetic field <span class="math notranslate nohighlight">\(\vec{B}\)</span>. This yields</p>
<div class="math notranslate nohighlight" id="equation-eqtoropolo">
<span class="eqno">(1)<a class="headerlink" href="#equation-eqtoropolo" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
\tilde{\rho} \vec{u} & = \vec{\nabla}\times(\vec{\nabla}\times W\,\vec{e_r}) +
\vec{\nabla}\times Z\,\vec{e_r}\,, \\
\vec{B} & = \vec{\nabla}\times(\vec{\nabla}\times g\,\vec{e_r}) +
\vec{\nabla}\times h\,\vec{e_r}\,.
\end{aligned}}\end{split}\]</div>
<p>The two scalar potentials of a divergence free vector field can be extracted
from its radial component and the radial component of its curl using the fact that
the toroidal field has not radial component:</p>
<div class="math notranslate nohighlight" id="equation-eqdeltah">
<span class="eqno">(2)<a class="headerlink" href="#equation-eqdeltah" title="Link to this equation">¶</a></span>\[\begin{split}\begin{aligned}
\vec{e_r}\cdot \tilde{\rho}\vec{u} &= - \nabla_H^2\,W, \\
\vec{e_r}\cdot\left(\vec{\nabla}\times\vec{u}\right) & = - \nabla_H^2\,Z,
\end{aligned}\end{split}\]</div>
<p>where the operator <span class="math notranslate nohighlight">\(\nabla_H^2\)</span> denotes the horizontal part of the Laplacian:</p>
<div class="math notranslate nohighlight" id="equation-eqlaplaceh">
<span class="eqno">(3)<a class="headerlink" href="#equation-eqlaplaceh" title="Link to this equation">¶</a></span>\[\nabla_H^2= \frac{1}{r^{2}\sin{\theta}}
\frac{\partial}{\partial\theta}\left(\sin{\theta}\frac{\partial}{\partial\theta}\right)
+ \frac{1}{r^{2}\sin^2{\theta}} \frac{\partial^{2}}{\partial^{2}\phi}.\]</div>
<p>The equation <a class="reference internal" href="#equation-eqtoropolo">(1)</a> can be expanded in spherical coordinates.
The three components of <span class="math notranslate nohighlight">\(\tilde{\rho}\vec{u}\)</span>
are given by</p>
<div class="math notranslate nohighlight" id="equation-eqtoropolo1">
<span class="eqno">(4)<a class="headerlink" href="#equation-eqtoropolo1" title="Link to this equation">¶</a></span>\[\tilde{\rho}\vec{u} = -(\nabla_H^2 W)\,\vec{e_r} + \left( \dfrac{1}{r}
\dfrac{\partial^2 W}{\partial r \partial \theta} +
\dfrac{1}{r\sin\theta}\dfrac{\partial Z}{\partial \phi}\right)\,\vec{e_\theta}
+\left(\dfrac{1}{r\sin\theta}\dfrac{\partial^2 W}{\partial r\partial \phi}-
\dfrac{1}{r}\dfrac{\partial Z}{\partial\theta} \right)\,\vec{e_\phi},\]</div>
<p>while the curl of <span class="math notranslate nohighlight">\(\tilde{\rho}\vec{u}\)</span> is expressed by</p>
<div class="math notranslate nohighlight" id="equation-eqtoropolo2">
<span class="eqno">(5)<a class="headerlink" href="#equation-eqtoropolo2" title="Link to this equation">¶</a></span>\[\begin{split}\begin{aligned}
\vec{\nabla}\times\tilde{\rho}\vec{u} = &-(\nabla_H^2 Z)\,\vec{e_r}+
\left[-\dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\phi}\left(
\dfrac{\partial^2}{\partial r^2}+\nabla_H^2 \right) W + \dfrac{1}{r}
\dfrac{\partial^2 Z}{\partial r\partial\theta}\right]\vec{e_\theta} \\
&+\left[\dfrac{1}{r}\dfrac{\partial}{\partial\theta}\left(
\dfrac{\partial^2}{\partial r^2}+\nabla_H^2\right) W + \dfrac{1}{r \sin\theta}
\dfrac{\partial^2 Z}{\partial r\partial\phi}\right]\vec{e_\phi},
\end{aligned}\end{split}\]</div>
<p>Using the horizontal part of the divergence operator</p>
<div class="math notranslate nohighlight">
\[\vec{\nabla}_H = \dfrac{1}{r\sin\theta} \dfrac{\partial}{\partial \theta}\sin\theta\;\vec{e}_\theta
+ \dfrac{1}{r\sin\theta} \dfrac{\partial}{\partial \phi}\;\vec{e}_\phi\]</div>
<p>above expressions can be simplified to</p>
<div class="math notranslate nohighlight">
\[\tilde{\rho}\vec{u} = -\nabla_H^2\;\vec{e_r}\; W + \vec{\nabla}_H \dfrac{\partial}{\partial r}\;W
+ \vec{\nabla}_H\times\vec{e}_r\;Z\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[\nabla\times\tilde{\rho}\vec{u} = -\nabla_H^2\;\vec{e}_r\;Z + \vec{\nabla}_H \dfrac{\partial}{\partial r}\;Z
- \vec{\nabla}_H\times\nabla_H^2\vec{e}_r\;W\;\;.\]</div>
<p>Below we will use the fact that the horizontal components of the poloidal field depend
on the radial derivative of the poloidal potential.</p>
</section>
<section id="spherical-harmonic-representation">
<h2>Spherical harmonic representation<a class="headerlink" href="#spherical-harmonic-representation" title="Link to this heading">¶</a></h2>
<p>Spherical harmonic functions <span class="math notranslate nohighlight">\(Y_\ell^m\)</span> are a natural choice for the
horizontal expansion in colatitude <span class="math notranslate nohighlight">\(\theta\)</span> and longitude <span class="math notranslate nohighlight">\(\phi\)</span>:</p>
<div class="math notranslate nohighlight">
\[Y_\ell^m(\theta,\phi) = P_{\ell}^m(\cos{\theta})\,e^{i m \phi},\]</div>
<p>where <span class="math notranslate nohighlight">\(\ell\)</span> and <span class="math notranslate nohighlight">\(m\)</span> denote spherical harmonic degree and order, respectively,
<span class="math notranslate nohighlight">\(P_\ell^m\)</span> is an associated Legendre function. Different normalization are in
use. Here we adopt a complete normalization so that the orthogonality relation
reads</p>
<div class="math notranslate nohighlight" id="equation-eqorthoylm">
<span class="eqno">(6)<a class="headerlink" href="#equation-eqorthoylm" title="Link to this equation">¶</a></span>\[\int_{0}^{2\pi} d\,\phi \int_{0}^{\pi}
\sin{\theta}\, d\theta\; Y_\ell^m(\theta,\phi)\,Y_{\ell^\prime}^{m^\prime}
(\theta,\phi) \; = \; \delta_{\ell \ell^\prime}\delta^{m m^\prime}.\]</div>
<p>This means that</p>
<div class="math notranslate nohighlight">
\[Y_{\ell}^{m}(\theta,\phi) = \left(\dfrac{(2\ell+1)}{4\pi}\dfrac{(\ell-|m|)!}{(\ell+|m|)!}\right)^{1/2}
P_\ell^m(\cos{\theta})\,e^{i m \phi},\]</div>
<p>As an example, the spherical harmonic representation of the
magnetic poloidal potential <span class="math notranslate nohighlight">\(g(r,\theta,\phi)\)</span>, truncated at degree and order
<span class="math notranslate nohighlight">\(\ell_{max}\)</span>, then reads</p>
<div class="math notranslate nohighlight" id="equation-eqspatspec">
<span class="eqno">(7)<a class="headerlink" href="#equation-eqspatspec" title="Link to this equation">¶</a></span>\[g(r,\theta,\phi) = \sum_{\ell=0}^{\ell_{max}}\sum_{m=-\ell}^{\ell} g_{\ell m}(r)\,Y_{\ell}^{m}(\theta,\phi),\]</div>
<p>with</p>
<div class="math notranslate nohighlight" id="equation-eqlegtf1">
<span class="eqno">(8)<a class="headerlink" href="#equation-eqlegtf1" title="Link to this equation">¶</a></span>\[g_{\ell m}(r) = \frac{1}{\pi}\,\int_{0}^{\pi} d \theta \sin{\theta}\; g_m(r,\theta)\;
P_\ell^m(\cos{\theta}),\]</div>
<div class="math notranslate nohighlight" id="equation-eqlegtf2">
<span class="eqno">(9)<a class="headerlink" href="#equation-eqlegtf2" title="Link to this equation">¶</a></span>\[g_{m}(r,\theta) = \frac{1}{2\pi}\,\int_{0}^{2\pi} d \phi\; g(r,\theta,\phi)\; e^{- i m \phi} .\]</div>
<p>The potential <span class="math notranslate nohighlight">\(g(r,\theta,\phi)\)</span> is a real function so that
<span class="math notranslate nohighlight">\(g_{\ell m}^\star(r)=g_{\ell,-m}(r)\)</span>, where the asterisk denotes the complex conjugate.
Thus, only coefficients with <span class="math notranslate nohighlight">\(m \ge 0\)</span> have to be considered. The same kind of
expansion is made for the toroidal magnetic potential, the mass flux potentials,
pressure, entropy (or temperature) and chemical composition.</p>
<p>The equations <a class="reference internal" href="#equation-eqlegtf1">(8)</a> and <a class="reference internal" href="#equation-eqlegtf2">(9)</a> define a two-step transform
from the longitude/latitude representation to the spherical harmonic
representation <span class="math notranslate nohighlight">\((r,\theta,\phi)\longrightarrow(r,\ell,m)\)</span>. The equation
<a class="reference internal" href="#equation-eqspatspec">(7)</a> formulates the inverse procedure
<span class="math notranslate nohighlight">\((r,\ell,m)\longrightarrow(r,\theta,\phi)\)</span>. Fast-Fourier transforms are
employed in the longitudinal direction, requiring (at least) <span class="math notranslate nohighlight">\(N_\phi = 2 \ell_{max}+1\)</span>
evenly spaced grid points <span class="math notranslate nohighlight">\(\phi_i\)</span>.
MagIC relies on the Gauss-Legendre quadrature for evaluating the integral
<a class="reference internal" href="#equation-eqlegtf1">(8)</a></p>
<div class="math notranslate nohighlight">
\[g_{\ell m}(r) = \frac{1}{N_{\theta}}
\sum_{j=1}^{N_{\theta}}\,w_j\,g_m(r,\theta_j)\; P_\ell^m(\cos{\theta_j}),\]</div>
<p>where <span class="math notranslate nohighlight">\(\theta_j\)</span> are the <span class="math notranslate nohighlight">\(N_{\theta}\)</span> Gaussian quadrature points
defining the latitudinal grid, and <span class="math notranslate nohighlight">\(w_j\)</span> are the respective weights. Pre-stored
values of the associated Legendre functions at grid points <span class="math notranslate nohighlight">\(\theta_j\)</span> in
combination with a FFT in <span class="math notranslate nohighlight">\(\phi\)</span> provide the inverse transform <a class="reference internal" href="#equation-eqspatspec">(7)</a>.
Generally, <span class="math notranslate nohighlight">\(N_\phi= 2 N_\theta\)</span> is chosen, which provides
isotropic resolution in the equatorial region. Choosing
<span class="math notranslate nohighlight">\(\ell_{max}= [ \min(2 N_\theta,N_\phi)-1]/3\)</span> prevents aliasing errors.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>In MagIC, the Legendre functions are defined in the subroutine
<a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/plms_theta/plm_theta" title="f/plms_theta/plm_theta"><code class="xref f f-subr docutils literal notranslate"><span class="pre">plm_theta</span></code></a>. The Legendre transforms
from spectral to grid space are computed in the module
<code class="xref f f-mod docutils literal notranslate"><span class="pre">legendre_spec_to_grid</span></code>, while the backward transform (from grid
space to spectral space) are computed in the module
<code class="xref f f-mod docutils literal notranslate"><span class="pre">legendre_grid_to_spec</span></code>. The fast Fourier transforms are computed
in the module <a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/fft" title="f/fft: This module contains the native subroutines used to compute FFT's. They are based on the FFT99 package from Temperton: http://www.cesm.ucar.edu/models/cesm1.2/cesm/cesmBbrowser/html_code/cam/fft99.F90.html"><code class="xref f f-mod docutils literal notranslate"><span class="pre">fft</span></code></a>.</p>
</div>
<section id="special-recurrence-relations">
<h3>Special recurrence relations<a class="headerlink" href="#special-recurrence-relations" title="Link to this heading">¶</a></h3>
<p>The action of a horizontal Laplacian <a class="reference internal" href="#equation-eqlaplaceh">(3)</a> on spherical harmonics can be
analytically expressed by</p>
<div class="math notranslate nohighlight" id="equation-eqhorizlaplylm">
<span class="eqno">(10)<a class="headerlink" href="#equation-eqhorizlaplylm" title="Link to this equation">¶</a></span>\[\boxed{
\nabla_H^2 Y_{\ell}^{m} = -\dfrac{\ell(\ell+1)}{r^2}\,Y_{\ell}^{m}\,.
}\]</div>
<p>They are several useful recurrence relations for the Legendre polynomials that will
be further employed to compute Coriolis forces and the <span class="math notranslate nohighlight">\(\theta\)</span> and <span class="math notranslate nohighlight">\(\phi\)</span>
derivatives of advection and Lorentz forces.
Four of these operators are used in <strong>MagIC</strong>. The first one is defined by</p>
<div class="math notranslate nohighlight">
\[\vartheta_1 = \dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\sin^2\theta
=\sin\theta\dfrac{\partial}{\partial\theta}+2\cos\theta\,.\]</div>
<p>The action of this operator on a Legendre polynomials is given by</p>
<div class="math notranslate nohighlight">
\[\vartheta_1 = (\ell+2)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta)
-(\ell-1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta),\]</div>
<p>where <span class="math notranslate nohighlight">\(c_\ell^m\)</span> is defined by</p>
<div class="math notranslate nohighlight" id="equation-eqclmop">
<span class="eqno">(11)<a class="headerlink" href="#equation-eqclmop" title="Link to this equation">¶</a></span>\[c_\ell^m = \sqrt{\dfrac{(\ell+m)(\ell-m)}{(2\ell-1)(2\ell+1)}}\,.\]</div>
<p><em>How is that implemented in the code?</em> Let’s assume we want the spherical harmonic contribution
of degree <span class="math notranslate nohighlight">\(\ell\)</span> and order <cite>m</cite> for the expression</p>
<div class="math notranslate nohighlight">
\[\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}(\sin\theta\,f(\theta))\,.\]</div>
<p>In order to employ the operator <span class="math notranslate nohighlight">\(\vartheta_1\)</span> for the derivative, we thus define a
new function</p>
<div class="math notranslate nohighlight">
\[F(\theta)=f(\theta)/\sin\theta\,,\]</div>
<p>so that</p>
<div class="math notranslate nohighlight">
\[\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}[\sin\theta\,f(\theta)]
=\vartheta_1 F(\theta)\,.\]</div>
<p>Expanding <span class="math notranslate nohighlight">\(F(\theta)\)</span> in Legendre polynomials and using the respective
orthogonality relation we can then map out the required contribution in the following way:</p>
<div class="math notranslate nohighlight" id="equation-eqoptheta1">
<span class="eqno">(12)<a class="headerlink" href="#equation-eqoptheta1" title="Link to this equation">¶</a></span>\[\boxed{
\int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_1\sum_{\ell'}F_{\ell'}^m P_{\ell'}^m
=(\ell+1)\,c_{\ell}^m\,F_{\ell-1}^m-\ell\,c_{\ell+1}^m\,F_{\ell+1}^m\,.}\]</div>
<p>Here, we have assumed that the Legendre functions are completely normalized such that</p>
<div class="math notranslate nohighlight">
\[\int_0^\pi d\theta\,\sin\theta\,P_\ell^m P_{\ell'}^m = \delta_{\ell \ell'}\,.\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>This operator is defined in the module <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data" title="f/horizontal_data: Module containing functions depending on longitude and latitude plus help arrays depending on degree and order"><code class="xref f f-mod docutils literal notranslate"><span class="pre">horizontal_data</span></code></a> by the variables
<a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta1s" title="f/horizontal_data/dtheta1s"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta1S</span></code></a> for the first part of the right-hand
side of <a class="reference internal" href="#equation-eqoptheta1">(12)</a> and <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta1a" title="f/horizontal_data/dtheta1a"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta1A</span></code></a> for the
second part.</p>
</div>
<p>The second operator used to formulate colatitude derivatives is</p>
<div class="math notranslate nohighlight">
\[\vartheta_2 = \sin\theta\dfrac{\partial}{\partial\theta}\,.\]</div>
<p>The action of this operator on the Legendre polynomials reads</p>
<div class="math notranslate nohighlight">
\[\vartheta_2 P_\ell^m(\cos\theta)=\ell\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta)
-(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]</div>
<p>so that</p>
<div class="math notranslate nohighlight" id="equation-eqoptheta2">
<span class="eqno">(13)<a class="headerlink" href="#equation-eqoptheta2" title="Link to this equation">¶</a></span>\[\boxed{
\int_0^\pi d\theta\,\sin\theta \,P_\ell^m\vartheta_2\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m
=(\ell-1)\,c_{\ell}^m\,f_{\ell-1}^m-(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>This operator is defined in the module <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data" title="f/horizontal_data: Module containing functions depending on longitude and latitude plus help arrays depending on degree and order"><code class="xref f f-mod docutils literal notranslate"><span class="pre">horizontal_data</span></code></a> by the variables
<a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta2s" title="f/horizontal_data/dtheta2s"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta2S</span></code></a> for the first part of the right-hand
side of <a class="reference internal" href="#equation-eqoptheta2">(13)</a> and <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta2a" title="f/horizontal_data/dtheta2a"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta2A</span></code></a> for the
second part.</p>
</div>
<p>The third combined operator is defined by:</p>
<div class="math notranslate nohighlight">
\[\vartheta_3 = \sin\theta\dfrac{\partial}{\partial\theta}+\cos\theta\,L_H\,,\]</div>
<p>where <span class="math notranslate nohighlight">\(-L_H/r^2=\nabla_H^2\)</span>.</p>
<p>Acting with <span class="math notranslate nohighlight">\(\vartheta_3\)</span> on a Legendre function gives:</p>
<div class="math notranslate nohighlight">
\[\vartheta_3 P_\ell^m(\cos\theta)=\ell(\ell+1)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta)
+(\ell-1)(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]</div>
<p>which results into:</p>
<div class="math notranslate nohighlight" id="equation-eqoptheta3">
<span class="eqno">(14)<a class="headerlink" href="#equation-eqoptheta3" title="Link to this equation">¶</a></span>\[\boxed{
\int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_3\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m
=(\ell-1)(\ell+1)\,c_{\ell}^m\,f_{\ell-1}^m+\ell(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>This operator is defined in the module <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data" title="f/horizontal_data: Module containing functions depending on longitude and latitude plus help arrays depending on degree and order"><code class="xref f f-mod docutils literal notranslate"><span class="pre">horizontal_data</span></code></a> by the variables
<a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta3s" title="f/horizontal_data/dtheta3s"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta3S</span></code></a> for the first part of the right-hand
side of <a class="reference internal" href="#equation-eqoptheta3">(14)</a> and <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta3a" title="f/horizontal_data/dtheta3a"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta3A</span></code></a> for the
second part.</p>
</div>
<p>The fourth (and last) combined operator is defined by:</p>
<div class="math notranslate nohighlight">
\[\vartheta_4 = \dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\sin^2\theta\,L_H
=\vartheta1\,L_H\,,\]</div>
<p>Acting with <span class="math notranslate nohighlight">\(\vartheta_3\)</span> on a Legendre function gives:</p>
<div class="math notranslate nohighlight">
\[\vartheta_4 P_\ell^m(\cos\theta)=\ell(\ell+1)(\ell+2)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta)
-\ell(\ell-1)(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]</div>
<p>which results into:</p>
<div class="math notranslate nohighlight" id="equation-eqoptheta4">
<span class="eqno">(15)<a class="headerlink" href="#equation-eqoptheta4" title="Link to this equation">¶</a></span>\[\boxed{
\int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_4\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m
=\ell(\ell-1)(\ell+1)\,c_{\ell}^m\,f_{\ell-1}^m-\ell(\ell+1)(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>This operator is defined in the module <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data" title="f/horizontal_data: Module containing functions depending on longitude and latitude plus help arrays depending on degree and order"><code class="xref f f-mod docutils literal notranslate"><span class="pre">horizontal_data</span></code></a> by the variables
<a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta4s" title="f/horizontal_data/dtheta4s"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta4S</span></code></a> for the first part of the right-hand
side of <a class="reference internal" href="#equation-eqoptheta4">(15)</a> and <a class="reference internal" href="apiFortran/precalcModules.html#f/horizontal_data/dtheta4a" title="f/horizontal_data/dtheta4a"><code class="xref f f-var docutils literal notranslate"><span class="pre">dTheta4A</span></code></a> for the
second part.</p>
</div>
</section>
</section>
<section id="radial-representation">
<h2>Radial representation<a class="headerlink" href="#radial-representation" title="Link to this heading">¶</a></h2>
<p>In MagIC, the radial dependencies are either expanded into complete sets of functions, the
Chebyshev polynomials <span class="math notranslate nohighlight">\({\cal C}(x)\)</span>; or discretized using finite differences. For the former approach, the Chebyshev polynomial of degree <span class="math notranslate nohighlight">\(n\)</span> is defined by</p>
<div class="math notranslate nohighlight">
\[{\cal C}_n(x)\approx\cos\left[n\,\arccos(x)\right]\quad -1\leq x \leq 1\,.\]</div>
<p>When truncating at degree <span class="math notranslate nohighlight">\(N\)</span>, the radial expansion of the poloidal
magnetic potential reads</p>
<div class="math notranslate nohighlight" id="equation-eqgridcheb">
<span class="eqno">(16)<a class="headerlink" href="#equation-eqgridcheb" title="Link to this equation">¶</a></span>\[g_{\ell m}(r) = \sum_{n=0}^{N} g_{\ell mn}\;{\cal C}_n(r) ,\]</div>
<p>with</p>
<div class="math notranslate nohighlight" id="equation-eqspeccheb">
<span class="eqno">(17)<a class="headerlink" href="#equation-eqspeccheb" title="Link to this equation">¶</a></span>\[g_{\ell mn} = \frac{2-\delta_{n0}}{\pi}\int_{-1}^{1}
\frac{d x\, g_{\ell m}(r(x))\;{\cal C}_n(x)}{\sqrt{1-x^2}} .\]</div>
<p>The Chebyshev definition space <span class="math notranslate nohighlight">\((-1\leq x\leq 1)\)</span> is then linearly mapped
onto a radius range <span class="math notranslate nohighlight">\((r_i\leq r \leq r_o)\)</span> by</p>
<div class="math notranslate nohighlight" id="equation-eqchebmap">
<span class="eqno">(18)<a class="headerlink" href="#equation-eqchebmap" title="Link to this equation">¶</a></span>\[x(r)= 2 \frac{r-r_i}{r_o-r_i} - 1 .\]</div>
<p>In addition, nonlinear mapping can be defined to modify the radial dependence of the
grid-point density.</p>
<p>When choosing the <span class="math notranslate nohighlight">\(N_r\)</span> extrema of <span class="math notranslate nohighlight">\({\cal C}_{N_r-1}\)</span> as radial grid points,</p>
<div class="math notranslate nohighlight" id="equation-eqchebgrid">
<span class="eqno">(19)<a class="headerlink" href="#equation-eqchebgrid" title="Link to this equation">¶</a></span>\[x_k=\cos{\left[\pi \frac{(k-1)}{N_r-1}\right]}\;\;\;,\;\;\; k=1,2,\ldots,N_r ,\]</div>
<p>the values of the Chebyshev polynomials at these points are simply given by
the cosine functions:</p>
<div class="math notranslate nohighlight">
\[{\cal C}_{nk} = {\cal C}_n(x_k)=\cos{\left[\pi \frac{ n (k-1)}{N_r-1}\right]} .\]</div>
<p>This particular choice has two advantages.
For one, the grid points become denser toward the inner and outer
radius and better resolve potential thermal and viscous boundary layers.
In addition, type I Discrete Cosine Transforms (DCTs) can be employed to switch between
grid representation <a class="reference internal" href="#equation-eqgridcheb">(16)</a> and Chebyshev representations <a class="reference internal" href="#equation-eqspeccheb">(17)</a>,
rendering this procedure a fast-Chebyshev transform.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The Chebyshev (Gauss-Lobatto) grid is defined in the module
<a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/chebyshev_polynoms_mod" title="f/chebyshev_polynoms_mod"><code class="xref f f-mod docutils literal notranslate"><span class="pre">chebyshev_polynoms_mod</span></code></a>. The cosine transforms are computed in the
modules <a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/cosine_transform_even" title="f/cosine_transform_even"><code class="xref f f-mod docutils literal notranslate"><span class="pre">cosine_transform_even</span></code></a> and <a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/fft_fac_mod" title="f/fft_fac_mod"><code class="xref f f-mod docutils literal notranslate"><span class="pre">fft_fac_mod</span></code></a>.</p>
</div>
</section>
<section id="spectral-equations">
<h2>Spectral equations<a class="headerlink" href="#spectral-equations" title="Link to this heading">¶</a></h2>
<p>We have now introduced the necessary tools for deriving the
spectral equations.
Taking the <strong>radial components</strong> of the Navier-Stokes equation
and the induction equation provides the equations
for the poloidal potentials <span class="math notranslate nohighlight">\(W(r,\theta,\phi)\)</span> and <span class="math notranslate nohighlight">\(g(r,\theta,\phi)\)</span>.
The <strong>radial component of the curl</strong> of these equations provides
the equations for the toroidal counterparts
<span class="math notranslate nohighlight">\(Z(r,\theta,\phi)\)</span> and <span class="math notranslate nohighlight">\(h(r,\theta,\phi)\)</span>.
The pressure remains an additional unknown. Hence one more equation
involving <span class="math notranslate nohighlight">\(W_{\ell mn}\)</span> and <span class="math notranslate nohighlight">\(p_{\ell mn}\)</span>
is required. It is obtained by taking the
<strong>horizontal divergence</strong> of the Navier-Stokes equation.</p>
<p>Expanding all potentials in spherical harmonics and Chebyshev polynomials,
multiplying with <span class="math notranslate nohighlight">\({Y_{\ell}^{m}}^\star\)</span>, and integrating over spherical surfaces
(while making use of
the orthogonality relation <a class="reference internal" href="#equation-eqorthoylm">(6)</a> results in equations for the
coefficients <span class="math notranslate nohighlight">\(W_{\ell mn}\)</span>, <span class="math notranslate nohighlight">\(Z_{\ell mn}\)</span>, <span class="math notranslate nohighlight">\(g_{\ell mn}\)</span>,
<span class="math notranslate nohighlight">\(h_{\ell mn}\)</span>, <span class="math notranslate nohighlight">\(P_{\ell mn}\)</span> and <span class="math notranslate nohighlight">\(s_{\ell mn}\)</span>,
respectively.</p>
<section id="equation-for-the-poloidal-potential-w">
<h3>Equation for the poloidal potential <span class="math notranslate nohighlight">\(W\)</span><a class="headerlink" href="#equation-for-the-poloidal-potential-w" title="Link to this heading">¶</a></h3>
<p>The temporal evolution of <span class="math notranslate nohighlight">\(W\)</span> is obtained by taking <span class="math notranslate nohighlight">\(\vec{e_r}\cdot\)</span> of each
term entering the Navier-Stokes equation. For the
time-derivative, one gets using <a class="reference internal" href="#equation-eqdeltah">(2)</a>:</p>
<div class="math notranslate nohighlight">
\[\tilde{\rho}\vec{e_r}\cdot\left(\dfrac{\partial \vec{u}}{\partial t}\right) =
\dfrac{\partial}{\partial t}(\vec{e_r}\cdot\tilde{\rho}\vec{u})=-\nabla_H^2\dfrac{\partial
W}{\partial t}.\]</div>
<p>For the viscosity term, one gets</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{e_r}\cdot\vec{\nabla}\cdot \mathsf{S} = & -\nu\,\nabla_H^2\left[\dfrac{\partial^2 W}
{\partial r^2}
+\left\lbrace 2\dfrac{d\ln\nu}{dr}-\dfrac{1}{3}\dfrac{d\ln\tilde{\rho}}{dr}\right\rbrace
\dfrac{\partial W}{\partial r} \right. \\
& -\left. \left\lbrace -\nabla_H^2 + \dfrac{4}{3}\left(\dfrac{d^2\ln\tilde{\rho}}{dr^2}
+\dfrac{d\ln\nu}{dr} \dfrac{d\ln\tilde{\rho}}{dr} +
\dfrac{1}{r}\left[3\dfrac{d\ln\nu}{dr}+
\dfrac{d\ln\tilde{\rho}}{dr}\right] \right) \right\rbrace W\right],
\end{aligned}\end{split}\]</div>
<div class="admonition note">
<p class="admonition-title">Note</p>
<p>In case of a constant kinematic viscosity, the <span class="math notranslate nohighlight">\(d\ln\nu/dr\)</span>
terms vanish. If in addition,the background density is constant, the
<span class="math notranslate nohighlight">\(d\ln\tilde{\rho}/dr\)</span> terms also vanish. In that Boussinesq limit, this
viscosity term would then be simplified as</p>
<div class="math notranslate nohighlight">
\[\vec{e_r}\cdot\vec{\nabla}^2 \vec{u} = -\nabla_H^2\left[\dfrac{\partial^2 W}{\partial r^2}
+\nabla_H^2\,W\right]\,.\]</div>
</div>
<p>Using Eq. <a class="reference internal" href="#equation-eqhorizlaplylm">(10)</a> then allows to finally write the time-evolution equation
for the poloidal potential <span class="math notranslate nohighlight">\(W_{\ell m n}\)</span>:</p>
<div class="math notranslate nohighlight" id="equation-eqspecw">
<span class="eqno">(20)<a class="headerlink" href="#equation-eqspecw" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
E\,\dfrac{\ell(\ell+1)}{r^2}\left[\left\lbrace\dfrac{\partial}{\partial t} +
\nu\,\dfrac{\ell(\ell+1)}{r^2} + \dfrac{4}{3}\,\nu\,\left(\dfrac{d^2\ln\tilde{\rho}}{dr^2}
+\dfrac{d\ln\nu}{dr} \dfrac{d\ln\tilde{\rho}}{dr} +
\dfrac{1}{r}\left[3\dfrac{d\ln\nu}{dr}+
\dfrac{d\ln\tilde{\rho}}{dr}\right] \right)\right\rbrace\right. & \,{\cal C}_n & \\
-\nu\,\left\lbrace 2\dfrac{d\ln\nu}{dr}-\dfrac{1}{3}\dfrac{d\ln\tilde{\rho}}{dr}\right\rbrace
&\,{\cal C}'_n & \\
-\nu & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& W_{\ell m n} \\
+ \left[{\cal C}'_n -\dfrac{d\ln\tilde{\rho}}{dr}{\cal C}_n\right] & & P_{\ell m n} \\
- \left[\dfrac{Ra\,E}{Pr}\,\tilde{\rho}\,g(r)\right] & \,{\cal C}_n & s_{\ell m n} \\
- \left[\dfrac{Ra_\xi\,E}{Sc}\,\tilde{\rho}\,g(r)\right] & \,{\cal C}_n & \xi_{\ell m n} \\
= {\cal N}^W_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^W =\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \vec{F}\,. & &
\end{aligned}}\end{split}\]</div>
<p>Here, <span class="math notranslate nohighlight">\(d\Omega\)</span> is the spherical surface element. We use the summation convention
for the Chebyshev index <span class="math notranslate nohighlight">\(n\)</span>. The radial derivatives of Chebyshev
polynomials are denoted by primes.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecw">(20)</a> are coded in
the subroutines <a class="reference internal" href="apiFortran/timeAdvLinear.html#f/updatewp_mod/get_wpmat" title="f/updatewp_mod/get_wpmat"><code class="xref f f-subr docutils literal notranslate"><span class="pre">get_wpMat</span></code></a></p>
</div>
</section>
<section id="equation-for-the-toroidal-potential-z">
<h3>Equation for the toroidal potential <span class="math notranslate nohighlight">\(Z\)</span><a class="headerlink" href="#equation-for-the-toroidal-potential-z" title="Link to this heading">¶</a></h3>
<p>The temporal evolution of <span class="math notranslate nohighlight">\(Z\)</span> is obtained by taking the radial component of the
curl of the Navier-Stokes equation (i.e. <span class="math notranslate nohighlight">\(\vec{e_r}\cdot\vec{\nabla}\times\)</span>). For
the time derivative, one gets using <a class="reference internal" href="#equation-eqdeltah">(2)</a>:</p>
<div class="math notranslate nohighlight">
\[\vec{e_r}\cdot\vec\nabla\times\left(\dfrac{\partial\tilde{\rho}\vec{u}}{\partial t}\right)=
\dfrac{\partial}{\partial t}(\vec{e_r}\cdot\vec{\nabla}\times\tilde{\rho}
\vec{u})=-\dfrac{\partial}{\partial t}(\nabla_H^2 Z) =
-\nabla_H^2\dfrac{\partial Z}{\partial t}\,.\]</div>
<p>The pressure gradient, one has</p>
<div class="math notranslate nohighlight">
\[\vec{\nabla}\times \left[\tilde{\rho}\vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}}\right)\right] =
\vec{\nabla} \tilde{\rho} \times \vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}}\right) +
\underbrace{\tilde{\rho} \vec{\nabla} \times \left[\vec{\nabla}\left( \dfrac{p'}{\tilde{\rho}}
\right)\right]}_{=0}\,.\]</div>
<p>This expression has no component along <span class="math notranslate nohighlight">\(\vec{e_r}\)</span>, as a consequence, there is
no pressure gradient contribution here. The
gravity term also vanishes as <span class="math notranslate nohighlight">\(\vec{\nabla}\times(\tilde{\rho}g(r)\vec{e_r})\)</span> has no
radial component.</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{e_r}\cdot\vec{\nabla}\times\left[\vec{\nabla}\cdot\mathsf{S}\right] = &
-\nu\,\nabla_H^2\left[\dfrac{\partial^2 Z}{\partial r^2}
+\left(\dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr}\right)\,\dfrac{\partial Z}{\partial r} \right.\\
& \left. - \left(\dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+
\dfrac{2}{r}\dfrac{d\ln\nu}{dr}+
\dfrac{d^2\ln\tilde{\rho}}{dr^2}+\dfrac{2}{r}
\dfrac{d\ln\tilde{\rho}}{dr}-\nabla_H^2\right) Z \right].
\end{aligned}\end{split}\]</div>
<div class="admonition note">
<p class="admonition-title">Note</p>
<p>Once again, this viscous term can be greatly simplified in the Boussinesq limit:</p>
<div class="math notranslate nohighlight">
\[\vec{e_r}\cdot\vec{\nabla}\times\left(\vec{\nabla}^2 \vec{u}\right) =
-\nabla_H^2\left[\dfrac{\partial^2 Z}{\partial r^2}
+\nabla_H^2\,Z\right]\,.\]</div>
</div>
<p>Using Eq. <a class="reference internal" href="#equation-eqhorizlaplylm">(10)</a> then allows to finally write the time-evolution equation
for the poloidal potential <span class="math notranslate nohighlight">\(Z_{\ell m n}\)</span>:</p>
<div class="math notranslate nohighlight" id="equation-eqspecz">
<span class="eqno">(21)<a class="headerlink" href="#equation-eqspecz" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
E\,\dfrac{\ell(\ell+1)}{r^2}\left[\left\lbrace\dfrac{\partial}{\partial t} +
\nu\,\dfrac{\ell(\ell+1)}{r^2} + \nu\,\left(\dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+
\dfrac{2}{r}\dfrac{d\ln\nu}{dr}+ \dfrac{d^2\ln\tilde{\rho}}{dr^2}+\dfrac{2}{r}
\dfrac{d\ln\tilde{\rho}}{dr}\right)\right\rbrace\right. & \,{\cal C}_n & \\
-\nu\,\left(\dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr}\right) &\,{\cal C}'_n & \\
-\nu & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& Z_{\ell m n} \\
= {\cal N}^Z_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^Z =
\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \left(\vec{\nabla}
\times\vec{F}\right)\,. & &
\end{aligned}}\end{split}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecz">(21)</a> are coded in
the subroutines <a class="reference internal" href="apiFortran/timeAdvLinear.html#f/updatez_mod/get_zmat" title="f/updatez_mod/get_zmat"><code class="xref f f-subr docutils literal notranslate"><span class="pre">get_zMat</span></code></a></p>
</div>
</section>
<section id="equation-for-pressure-p">
<h3>Equation for pressure <span class="math notranslate nohighlight">\(P\)</span><a class="headerlink" href="#equation-for-pressure-p" title="Link to this heading">¶</a></h3>
<p>The evolution of equation for pressure is obtained by taking the horizontal
divergence (i.e. <span class="math notranslate nohighlight">\(\vec{\nabla}_H\cdot\)</span>)
of the Navier-Stokes equation. This operator is defined such
that</p>
<div class="math notranslate nohighlight">
\[\vec{\nabla}_H\cdot\vec{a} = r\sin \dfrac{\partial (\sin\theta\,a_\theta)}{\partial \theta}
+r\sin \dfrac{\partial a_\phi}{\partial \phi}.\]</div>
<p>This relates to the total divergence via:</p>
<div class="math notranslate nohighlight">
\[\vec{\nabla}\cdot\vec{a}= \dfrac{1}{r^2}\dfrac{\partial(r^2 a_r)}{\partial r}+
\vec{\nabla}_H\cdot\vec{a}.\]</div>
<p>The time-derivative term is thus expressed by</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{\nabla}_H\cdot\left(\tilde{\rho}\dfrac{\partial \vec{u}}{\partial t}\right)
&= \dfrac{\partial}{\partial t}\left[\vec{\nabla}_H\cdot(\tilde{\rho}\vec{u}
)\right], \\
& = \dfrac{\partial}{\partial t}\left[\vec{\nabla}\cdot(\tilde{\rho}\vec{u})
-\dfrac{1}{r^2}\dfrac{\partial(r^2\tilde{\rho} u_r)}{\partial r}\right], \\
& = -\dfrac{\partial}{\partial t}\left[\dfrac{\partial (\tilde{\rho} u_r)}{\partial r}
+\dfrac{2\tilde{\rho} u_r}{r}\right], \\
& = \dfrac{\partial}{\partial t}\left[\dfrac{\partial (\nabla_H^2 W)}{\partial r}
+\dfrac{2}{r}\nabla_H^2 W\right], \\
& = \nabla_H^2\dfrac{\partial}{\partial t}\left(\dfrac{\partial W}{\partial r}\right).
\end{aligned}\end{split}\]</div>
<p>We note that the gravity term vanishes since <span class="math notranslate nohighlight">\(\vec{\nabla}_H\cdot(\tilde{\rho}
g(r)\vec{e_r}) = 0\)</span>. Concerning the pressure gradient, one has</p>
<div class="math notranslate nohighlight">
\[-\vec{\nabla}_H\cdot\left[\tilde{\rho} \vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}}
\right)\right] = -\left\lbrace\vec{\nabla}\cdot\left[\tilde{\rho} \vec{\nabla}
\left(\dfrac{p'}{\tilde{\rho}}\right)\right]-
\dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left[ r^2 \tilde{\rho}
\dfrac{\partial}{\partial r}\left(\dfrac{p'}{\tilde{\rho}}\right)\right] \right\rbrace =
-\nabla_H^2 \, p'.\]</div>
<p>The viscosity term then reads</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{\nabla}_H\cdot \left( \vec{\nabla}\cdot\mathsf{S} \right) = & \nu\,\nabla_H^2\left[
\dfrac{\partial^3 W}{\partial r^3} + \left(\dfrac{d\ln\nu}{dr}-
\dfrac{d\ln\tilde{\rho}}{dr}\right) \dfrac{\partial^2 W}{\partial r^2} \right. \\
& - \left[\dfrac{d^2\ln\tilde{\rho}}{dr^2} + \dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+
\dfrac{2}{r}\left(\dfrac{d\ln\nu}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}\right)
-\nabla_H^2 \right]\dfrac{\partial W}{\partial r} \\
& \left. -\left( \dfrac{2}{3}\dfrac{d\ln\tilde{\rho}}{dr}+\dfrac{2}{r}+\dfrac{d\ln\nu}{dr}
\right)\nabla_H^2\,W \right].
\end{aligned}\end{split}\]</div>
<div class="admonition note">
<p class="admonition-title">Note</p>
<p>Once again, this viscous term can be greatly simplified in the Boussinesq limit:</p>
<div class="math notranslate nohighlight">
\[\vec{\nabla}_H\cdot\left(\nabla^2 \vec{u}\right) =
-\nabla_H^2\left[\dfrac{\partial^3 W}{\partial r^3}
+\nabla_H^2\,\dfrac{\partial W}{\partial r}-\dfrac{2}{r}\nabla_H^2\,W\right]\,.\]</div>
</div>
<p>Using Eq. <a class="reference internal" href="#equation-eqhorizlaplylm">(10)</a> then allows to finally write the equation for the pressure
<span class="math notranslate nohighlight">\(P_{\ell m n}\)</span>:</p>
<div class="math notranslate nohighlight" id="equation-eqspecp">
<span class="eqno">(22)<a class="headerlink" href="#equation-eqspecp" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
E\,\dfrac{\ell(\ell+1)}{r^2}\left[
-\nu\,\left( \dfrac{2}{3}\dfrac{d\ln\tilde{\rho}}{dr}+\dfrac{2}{r}+\dfrac{d\ln\nu}{dr}
\right)\dfrac{\ell(\ell+1)}{r^2} \right.
& \,{\cal C}_n & \\
\left\lbrace\dfrac{\partial}{\partial t} +
\nu\,\dfrac{\ell(\ell+1)}{r^2} + \nu\,\left[\dfrac{d^2\ln\tilde{\rho}}{dr^2}+
\dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+
\dfrac{2}{r}\left(\dfrac{d\ln\nu}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}\right)\right]\right\rbrace
& \,{\cal C}'_n & \\
-\nu\,\left( \dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr}
\right) &\,{\cal C}''_n & \\
-\nu & \,{\cal C}'''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& W_{\ell m n} \\
+ \left[\dfrac{\ell(\ell+1)}{r^2}\right] & \,{\cal C}_n & P_{\ell m n} \\
= {\cal N}^P_{\ell m} = -\int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^P=-\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{\nabla}_H\cdot\vec{F}\,. & &
\end{aligned}}\end{split}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecp">(22)</a> are coded in
the subroutines <code class="xref f f-subr docutils literal notranslate"><span class="pre">get_wpMat</span></code></p>
</div>
<div class="admonition note">
<p class="admonition-title">Note</p>
<p>We note that the terms on the left hand side of <a class="reference internal" href="#equation-eqspecw">(20)</a>, <a class="reference internal" href="#equation-eqspecz">(21)</a> and
<a class="reference internal" href="#equation-eqspecp">(22)</a> resulting from the viscous term, the pressure gradient,
the buoyancy term, and the explicit time derivative completely decouple
in spherical harmonic degree and order.</p>
<p>The terms that do not decouple, namely Coriolis force, Lorentz force and
advection of momentum, are collected on the right-hand side
of <a class="reference internal" href="#equation-eqspecw">(20)</a>, <a class="reference internal" href="#equation-eqspecz">(21)</a> and <a class="reference internal" href="#equation-eqspecp">(22)</a> into the forcing term
<span class="math notranslate nohighlight">\(\vec{F}\)</span>:</p>
<div class="math notranslate nohighlight" id="equation-eqforcing">
<span class="eqno">(23)<a class="headerlink" href="#equation-eqforcing" title="Link to this equation">¶</a></span>\[\vec{F}=-2\,\tilde{\rho}\,\vec{e_z}\times\vec{u} - E\,\tilde{\rho}\,
\vec{u}\cdot\vec{\nabla}\,\vec{u}
+\frac{1}{Pm}\left(\vec{\nabla}\times\vec{B}\right)\times\vec{B}\,.\]</div>
</div>
<p>Resolving <span class="math notranslate nohighlight">\(\vec{F}\)</span> into potential functions is not required. Its
numerical evaluation is discussed <a class="reference internal" href="#secnonlineareqs"><span class="std std-ref">below</span></a>.</p>
</section>
<section id="equation-for-entropy-s">
<h3>Equation for entropy <span class="math notranslate nohighlight">\(s\)</span><a class="headerlink" href="#equation-for-entropy-s" title="Link to this heading">¶</a></h3>
<p>The equation for the entropy (or temperature in the Boussinesq limit) is given by</p>
<div class="math notranslate nohighlight" id="equation-eqspecs">
<span class="eqno">(24)<a class="headerlink" href="#equation-eqspecs" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
\dfrac{1}{Pr}\left[\left(Pr\dfrac{\partial}{\partial t} +
\kappa\,\dfrac{\ell(\ell+1)}{r^2}
\right)\right. & \,{\cal C}_n & \\
-\kappa\,\left(\dfrac{d\ln\kappa}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}+
+\dfrac{dln\tilde{T}}{dr}+\dfrac{2}{r}\right)
&\,{\cal C}'_n & \\
-\kappa & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& s_{\ell m n} \\
= {\cal N}^S_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^S = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\left[-\vec{u}\cdot\vec{\nabla}s+
\dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}\left(\Phi_\nu+
\dfrac{\lambda}{Pm^2\,E}\,j^2\right) \right]\,. & &
\end{aligned}}\end{split}\]</div>
<p>In this expression, <span class="math notranslate nohighlight">\(j=\vec{\nabla}\times\vec{B}\)</span> is the current. Once again,
the numerical evaluation of the right-hand-side (i.e. the non-linear terms) is
discussed <a class="reference internal" href="#secnonlinears"><span class="std std-ref">below</span></a>.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecs">(24)</a> are coded in
the subroutines <code class="xref f f-subr docutils literal notranslate"><span class="pre">get_sMat</span></code></p>
</div>
</section>
<section id="equation-for-chemical-composition-xi">
<h3>Equation for chemical composition <span class="math notranslate nohighlight">\(\xi\)</span><a class="headerlink" href="#equation-for-chemical-composition-xi" title="Link to this heading">¶</a></h3>
<p>The equation for the chemical composition is given by</p>
<div class="math notranslate nohighlight" id="equation-eqspecxi">
<span class="eqno">(25)<a class="headerlink" href="#equation-eqspecxi" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
\dfrac{1}{Sc}\left[\left(Sc\dfrac{\partial}{\partial t} +
\kappa_\xi\,\dfrac{\ell(\ell+1)}{r^2}
\right)\right. & \,{\cal C}_n & \\
-\kappa_\xi\,\left(\dfrac{d\ln\kappa_\xi}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}+
+\dfrac{2}{r}\right)
&\,{\cal C}'_n & \\
-\kappa_\xi & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& \xi_{\ell m n} \\
= {\cal N}^\xi_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^\xi = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\left[-\vec{u}\cdot\vec{\nabla}\xi
\right]\,. & &
\end{aligned}}\end{split}\]</div>
<p>Once again, the numerical evaluation of the right-hand-side (i.e. the
non-linear term) is discussed <a class="reference internal" href="#secnonlinearxi"><span class="std std-ref">below</span></a>.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecxi">(25)</a> are coded in
the subroutines <a class="reference internal" href="apiFortran/timeAdvLinear.html#f/updatexi_mod/get_ximat" title="f/updatexi_mod/get_ximat"><code class="xref f f-subr docutils literal notranslate"><span class="pre">get_xiMat</span></code></a></p>
</div>
</section>
<section id="equation-for-the-poloidal-magnetic-potential-g">
<h3>Equation for the poloidal magnetic potential <span class="math notranslate nohighlight">\(g\)</span><a class="headerlink" href="#equation-for-the-poloidal-magnetic-potential-g" title="Link to this heading">¶</a></h3>
<p>The equation for the poloidal magnetic potential is the radial
component of the dynamo equation since</p>
<div class="math notranslate nohighlight">
\[\vec{e_r}\cdot\left(\dfrac{\partial \vec{B}}{\partial t}\right) =
\dfrac{\partial}{\partial t}(\vec{e_r}\cdot\vec{B})=-\nabla_H^2\dfrac{\partial
g}{\partial t}.\]</div>
<p>The spectral form then reads</p>
<div class="math notranslate nohighlight" id="equation-eqspecg">
<span class="eqno">(26)<a class="headerlink" href="#equation-eqspecg" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
\dfrac{\ell(\ell+1)}{r^2}\left[\left(\dfrac{\partial}{\partial t} +
\dfrac{1}{Pm}\lambda\,\dfrac{\ell(\ell+1)}{r^2}
\right)\right. & \,{\cal C}_n & \\
-\dfrac{1}{Pm}\,\lambda & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& g_{\ell m n} \\
= {\cal N}^g_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^g=\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \vec{D}\,. & &
\end{aligned}}\end{split}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspecg">(26)</a> are coded in
the subroutines <a class="reference internal" href="apiFortran/timeAdvLinear.html#f/updateb_mod/get_bmat" title="f/updateb_mod/get_bmat"><code class="xref f f-subr docutils literal notranslate"><span class="pre">get_bMat</span></code></a></p>
</div>
</section>
<section id="equation-for-the-toroidal-magnetic-potential-h">
<h3>Equation for the toroidal magnetic potential <span class="math notranslate nohighlight">\(h\)</span><a class="headerlink" href="#equation-for-the-toroidal-magnetic-potential-h" title="Link to this heading">¶</a></h3>
<p>The equation for the toroidal magnetic field coefficient reads</p>
<div class="math notranslate nohighlight" id="equation-eqspech">
<span class="eqno">(27)<a class="headerlink" href="#equation-eqspech" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
\dfrac{\ell(\ell+1)}{r^2}\left[\left(\dfrac{\partial}{\partial t} +
\dfrac{1}{Pm}\lambda\,\dfrac{\ell(\ell+1)}{r^2}
\right)\right. & \,{\cal C}_n & \\
-\dfrac{1}{Pm}\,\dfrac{d\lambda}{dr} &\,{\cal C}'_n & \\
-\dfrac{1}{Pm}\,\lambda & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& h_{\ell m n} \\
= {\cal N}^h_{\ell m}= \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^h = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \left(\vec{\nabla}\times \vec{D}\right)\,. & &
\end{aligned}}\end{split}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The exact computation of the linear terms of <a class="reference internal" href="#equation-eqspech">(27)</a> are coded in
the subroutines <a class="reference internal" href="apiFortran/timeAdvLinear.html#f/updateb_mod/get_bmat" title="f/updateb_mod/get_bmat"><code class="xref f f-subr docutils literal notranslate"><span class="pre">get_bMat</span></code></a></p>
</div>
<div class="admonition note">
<p class="admonition-title">Note</p>
<p>We note that the terms on the left hand side of <a class="reference internal" href="#equation-eqspecg">(26)</a> and <a class="reference internal" href="#equation-eqspech">(27)</a>
resulting from the magnetic diffusion term
and the explicit time derivative completely decouple
in spherical harmonic degree and order.</p>
<p>The dynamo term does not decouple:</p>
<div class="math notranslate nohighlight" id="equation-eqdynamoterm">
<span class="eqno">(28)<a class="headerlink" href="#equation-eqdynamoterm" title="Link to this equation">¶</a></span>\[\vec{D}=\vec{\nabla}\times\left(\vec{u}\times\vec{B}\right)\,.\]</div>
</div>
<p>We have now derived a full set of equations
<a class="reference internal" href="#equation-eqspecw">(20)</a>, <a class="reference internal" href="#equation-eqspecz">(21)</a>, <a class="reference internal" href="#equation-eqspecp">(22)</a>, <a class="reference internal" href="#equation-eqspecs">(24)</a>, <a class="reference internal" href="#equation-eqspecg">(26)</a> and
<a class="reference internal" href="#equation-eqspech">(27)</a>,
each describing the evolution of a single spherical harmonic mode of the
six unknown fields (assuming that the terms on the right hand side
are given). Each equation couples <span class="math notranslate nohighlight">\(N+1\)</span> Chebyshev coefficients
for a given spherical harmonic mode <span class="math notranslate nohighlight">\((\ell,m)\)</span>.
Typically, a collocation method is employed to solve for the Chebyshev coefficients.
This means that the equations are required to be exactly satisfied at <span class="math notranslate nohighlight">\(N-1\)</span>
grid points defined by the equations <a class="reference internal" href="#equation-eqchebmap">(18)</a> and <a class="reference internal" href="#equation-eqchebgrid">(19)</a>.
Excluded are the points <span class="math notranslate nohighlight">\(r=r_i\)</span> and <span class="math notranslate nohighlight">\(r=r_o\)</span>, where the
<a class="reference internal" href="#secboundaryconditions"><span class="std std-ref">boundary conditions</span></a> provide
additional constraints on the set of Chebyshev coefficients.</p>
</section>
</section>
<section id="time-stepping-schemes">
<h2>Time-stepping schemes<a class="headerlink" href="#time-stepping-schemes" title="Link to this heading">¶</a></h2>
<p>Implicit time stepping schemes theoretically offer increased stability and
allow for larger time steps.
However, fully implicit approaches have the disadvantage that
the nonlinear-terms couple all spherical harmonic modes.
The potential gain in computational speed is therefore lost at
higher resolution, where one very large matrix has to be dealt with
rather than a set of much smaller ones.
Similar considerations hold for the Coriolis force, one of
the dominating forces in the system and therefore a prime candidate for
implicit treatment. However, the Coriolis term couples modes <span class="math notranslate nohighlight">\((\ell,m,n)\)</span> with
<span class="math notranslate nohighlight">\((\ell+1,m,n)\)</span> and <span class="math notranslate nohighlight">\((\ell-1,m,n)\)</span> and also couples poloidal and
toroidal flow potentials. An implicit treatment of the Coriolis term therefore
also results in a much larger (albeit sparse) inversion matrix.</p>
<p>We consequently adopt in <strong>MagIC</strong> a mixed implicit/explicit algorithm.
The general differential equation in time can be written in the form</p>
<div class="math notranslate nohighlight" id="equation-eqtstep">
<span class="eqno">(29)<a class="headerlink" href="#equation-eqtstep" title="Link to this equation">¶</a></span>\[\dfrac{\partial }{\partial t} y = \mathcal{I}(y,t) + \mathcal{E}(y,t),\quad y(t_o)=y_o\;.\]</div>
<p>where <span class="math notranslate nohighlight">\(\mathcal{I}\)</span> denotes the terms treated in an implicit time step
and <span class="math notranslate nohighlight">\(\mathcal{E}\)</span> the terms treated explicitly, i.e. the nonlinear and Coriolis contributions. In MagIC, two families of time-stepping schemes are supported: IMEX multistep and IMEX multistage methods.</p>
<p>First of all, the IMEX multistep methods correspond to time schemes where the solution results from the combination of several previous steps (such as for instance the Crank-Nicolson/Adams-Bashforth scheme). In this case, a general <span class="math notranslate nohighlight">\(k\)</span>-step IMEX multistep scheme reads</p>
<div class="math notranslate nohighlight">
\[\left(I-b_o^{\mathcal I} \delta t\,\mathcal{I}\right)y_{n+1}=\sum_{j=1}^k a_j y_{n+1-j}+\delta t\sum_{j=1}^k \left(b_j^\mathcal{E} \mathcal{E}_{n+1-j}+b_{j}^\mathcal{I}\mathcal{I}_{n+1-j}\right)\,,\]</div>
<p>where <span class="math notranslate nohighlight">\(I\)</span> is the identity matrix. The vectors <span class="math notranslate nohighlight">\(\vec{a}\)</span>, <span class="math notranslate nohighlight">\(\vec{b}^\mathcal{E}\)</span> and <span class="math notranslate nohighlight">\(\vec{b}^\mathcal{I}\)</span> correspond to the weighting factors of an IMEX multistep scheme. For instance, the commonly-used second-order scheme assembled from the combination of a Crank-Nicolson for the implicit terms and a second-order Adams-Bashforth for the explicit terms (hereafter CNAB2) corresponds to the following vectors: <span class="math notranslate nohighlight">\(\vec{a}=(1,0)\)</span>, <span class="math notranslate nohighlight">\(\vec{b}^{\mathcal{I}}=(1/2,1/2)\)</span> and <span class="math notranslate nohighlight">\(\vec{b}^{\mathcal{E}}=(3/2,-1/2)\)</span> for a constant timestep size <span class="math notranslate nohighlight">\(\delta t\)</span>.</p>
<p>In addition to CNAB2, MagIC supports several semi-implicit backward differentiation schemes of second, third and fourth order that are known to have good stability properties (heareafter SBDF2, SBDF3 and SBDF4), a modified CNAB2 from <a class="reference external" href="https://doi.org/10.1137/0732037">Ascher et al. (1995)</a> (termed MODCNAB) and the CNLF scheme (combination of Crank-Nicolson and Leap-Frog for the explicit terms).</p>
<p>MagIC also supports several IMEX Runge-Kutta multistage methods, frequently
called DIRK, an acronym that stands for <em>Diagonally Implicit Runge Kutta</em>. For
such schemes, the equation <a class="reference internal" href="#equation-eqtstep">(29)</a> is time-advanced from <span class="math notranslate nohighlight">\(t_n\)</span> to <span class="math notranslate nohighlight">\(t_{n+1}\)</span> by solving <span class="math notranslate nohighlight">\(\nu\)</span> sub-stages</p>
<div class="math notranslate nohighlight">
\[(I-a_{ii}^\mathcal{I}\delta t \mathcal{I})y_{i}=y_n+\delta t \sum_{j=1}^{i-1}\left(a_{i,j}^{\mathcal{E}}\mathcal{E}_j+a_{i,j}^\mathcal{I}\mathcal{I}_j\right), \quad 1\leq i\leq \nu,\]</div>
<p>where <span class="math notranslate nohighlight">\(y_i\)</span> is the intermediate solution at the stage <span class="math notranslate nohighlight">\(i\)</span>. The matrices <span class="math notranslate nohighlight">\(a_{i,j}^\mathcal{E}\)</span> and <span class="math notranslate nohighlight">\(a_{i,j}^\mathcal{I}\)</span> constitute the so-called Butcher tables that correspond to a DIRK scheme. MagIC supports several second and third order schemes: ARS222 and ARS443 from <a class="reference external" href="https://doi.org/10.1016/S0168-9274(97)00056-1">Ascher et al. (1997)</a>, LZ232 from <a class="reference external" href="https://doi.org/10.1016/j.cam.2005.02.020">Liu and Zou (2006)</a>, PC2 from <a class="reference external" href="https://doi.org/10.2514/6.1981-1259">Jameson et al. (1981)</a> and BPR353 from <a class="reference external" href="https://doi.org/10.1137/110842855">Boscarino et al. (2013)</a>.</p>
<p>In the code the equation <a class="reference internal" href="#equation-eqtstep">(29)</a> is formulated for each unknown spectral coefficient
(expect pressure) of spherical harmonic degree <span class="math notranslate nohighlight">\(\ell\)</span> and order <span class="math notranslate nohighlight">\(m\)</span>
and for each radial grid point <span class="math notranslate nohighlight">\(r_k\)</span>.
Because non-linear terms and the Coriolis force are treated explicitly,
the equations decouple for all spherical harmonic modes.
The different radial grid points, however, couple via the
Chebyshev modes and form a linear algebraic system of equations that can
be solved with standard methods for the different spectral contributions.</p>
<p>For example the respective system of equations for the poloidal magnetic potential <span class="math notranslate nohighlight">\(g\)</span> time advanced with a CNAB2 reads</p>
<div class="math notranslate nohighlight" id="equation-imex">
<span class="eqno">(30)<a class="headerlink" href="#equation-imex" title="Link to this equation">¶</a></span>\[\left( \mathcal{A}_{kn} - \dfrac{1}{2}\,\delta t\,\mathcal{I}_{kn}\right)\;g_{\ell mn}(t+\delta t) =
\left( \mathcal{A}_{kn} + \dfrac{1}{2}\,\delta t\,\mathcal{I}_{kn} \right)\;g_{\ell mn}(t) +
\frac{3}{2}\,\delta t\,\mathcal{E}_{k\ell m}(t) - \frac{1}{2}\,\delta t\,\mathcal{E}_{k\ell m}(t-\delta t)\]</div>
<p>with</p>
<div class="math notranslate nohighlight">
\[\mathcal{A}_{kn} = \dfrac{\ell (\ell+1)}{r_k^2}\,{\cal C}_{nk}\;,\]</div>
<div class="math notranslate nohighlight">
\[\mathcal{I}_{kn}=\dfrac{\ell(\ell+1)}{r_k^2}\,\dfrac{1}{Pm}\left({\cal C}''_{nk} - \dfrac{\ell(\ell+1)}{r_k^2}\;
{\cal C}_{nk}\right)\;,\]</div>
<p>and <span class="math notranslate nohighlight">\({\cal C}_{nk}={\cal C}_n(r_k)\)</span>.
<span class="math notranslate nohighlight">\(\mathcal{A}_{kn}\)</span> is a matrix that converts the poloidal field modes <span class="math notranslate nohighlight">\(g_{\ell mn}\)</span>
to the radial magnetic field <span class="math notranslate nohighlight">\(B_r(r_k,\ell,m)\)</span> for a given spherical harmonic contribution.</p>
<p>Here <span class="math notranslate nohighlight">\(k\)</span> and <span class="math notranslate nohighlight">\(n\)</span> number the radial grid points and the Chebyshev coefficients, respectively.
Note that the Einstein sum convention is used for Chebyshev modes <span class="math notranslate nohighlight">\(n\)</span>.</p>
<p><span class="math notranslate nohighlight">\(\mathcal{I}_{kn}\)</span> is the matrix describing the implicit contribution which is purely diffusive here.
Neither <span class="math notranslate nohighlight">\(\mathcal{A}_{kn}\)</span> nor <span class="math notranslate nohighlight">\(\mathcal{I}_{kn}\)</span> depend on time but the former
needs to be updated when the time step <span class="math notranslate nohighlight">\(\delta t\)</span> is changed.
The only explicit contribution is the nonlinear dynamo term</p>
<div class="math notranslate nohighlight">
\[\mathcal{E}_{k\ell m}(t)= {\cal N}_{k\ell m}^g = \int d\Omega\; {Y_{\ell}^{m}}^\star\;
\vec{e_r} \cdot \vec{D}(t,r_k,\theta,\phi)\;\; .\]</div>
<p><span class="math notranslate nohighlight">\(\mathcal{E}_{k\ell m}\)</span> is a one dimensional vector for all spherical harmonic combinations
<span class="math notranslate nohighlight">\(\ell m\)</span>.</p>
<p><strong>Courant’s condition</strong> offers a guideline
concerning the value of <span class="math notranslate nohighlight">\(\delta t\)</span>, demanding that <span class="math notranslate nohighlight">\(\delta t\)</span> should be smaller
than the advection time between two grid points. Strong Lorentz forces require
an additional stability criterion that is obtained by replacing the flow speed
by Alfvén’s velocity in a modified Courant criterion.
The explicit treatment of the Coriolis force requires that the time step is
limited to a fraction of the rotation period, which may be the relevant
criterion at low Ekman number when flow and magnetic field remain weak.
Non-homogeneous grids and other numerical effects generally require an
additional safety factor in the choice of <span class="math notranslate nohighlight">\(\delta t\)</span>.</p>
</section>
<section id="coriolis-force-and-nonlinear-terms">
<span id="secnonlineareqs"></span><h2>Coriolis force and nonlinear terms<a class="headerlink" href="#coriolis-force-and-nonlinear-terms" title="Link to this heading">¶</a></h2>
<section id="nonlinear-terms-entering-the-equation-for-w">
<span id="secnonlinearw"></span><h3>Nonlinear terms entering the equation for <span class="math notranslate nohighlight">\(W\)</span><a class="headerlink" href="#nonlinear-terms-entering-the-equation-for-w" title="Link to this heading">¶</a></h3>
<p>The nonlinear term <span class="math notranslate nohighlight">\({\cal N}^W\)</span> that enters the equation for the poloidal potential
<a class="reference internal" href="#equation-eqspecw">(20)</a> contains the radial component of the advection, the Lorentz force
and Coriolis force. In spherical coordinate, the two first contributions read:</p>
<div class="math notranslate nohighlight" id="equation-eqadvection">
<span class="eqno">(31)<a class="headerlink" href="#equation-eqadvection" title="Link to this equation">¶</a></span>\[\begin{split}\tilde{\rho}\left(\vec{u}\cdot\vec{\nabla}\vec{u}\right)=
\left\lbrace
\begin{aligned}
{\cal A}_r \\
{\cal A}_\theta \\
{\cal A}_\phi
\end{aligned}
\right\rbrace
=
\left\lbrace
\begin{aligned}
-\tilde{\rho}\,E\,\left(
u_r\dfrac{\partial u_r}{\partial r}+
\dfrac{u_\theta}{r}\dfrac{\partial u_r}{\partial \theta}+
\dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_r}{\partial \phi}
-\dfrac{u_\theta^2+u_\phi^2}{r}\right)+
\dfrac{1}{Pm}\left(j_\theta\,B_\phi-j_\phi\,B_\theta\right)\, , \\
-\tilde{\rho}\,E\,\left(
u_r\dfrac{\partial u_\theta}{\partial r}+
\dfrac{u_\theta}{r}\dfrac{\partial u_\theta}{\partial \theta} +
\dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_\theta}{\partial \phi}+
\dfrac{u_r u_\theta}{r}-\dfrac{\cos\theta}{r\sin\theta}u_\phi^2\right)+
\dfrac{1}{Pm}\left(j_\phi\,B_r-j_r\,B_\phi\right)\, ,\\
-\tilde{\rho}\,E\,\left(
u_r\dfrac{\partial u_\phi}{\partial r}+
\dfrac{u_\theta}{r}\dfrac{\partial u_\phi}{\partial \theta} +
\dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_\phi}{\partial \phi}+
\dfrac{u_r u_\phi}{r} +\dfrac{\cos\theta}{r\sin\theta}u_\theta u_\phi\right)+
\dfrac{1}{Pm}\left(j_r\,B_\theta-j_\theta\,B_r\right)\, ,
\end{aligned}
\right\rbrace\end{split}\]</div>
<p>The Coriolis force can be expressed as a function of the potentials <span class="math notranslate nohighlight">\(W\)</span> and
<span class="math notranslate nohighlight">\(Z\)</span> using <a class="reference internal" href="#equation-eqtoropolo1">(4)</a></p>
<div class="math notranslate nohighlight">
\[2\tilde{\rho} \vec{e_r}\cdot(\vec{u}\times\vec{e_z})=2\sin\theta\,\tilde{\rho}
u_\phi=\dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial \phi}-\sin\theta
\dfrac{\partial Z}{\partial \theta}\right)\,.\]</div>
<p>The nonlinear terms that enter the equation for the poloidal potential <a class="reference internal" href="#equation-eqspecw">(20)</a> thus
reads:</p>
<div class="math notranslate nohighlight">
\[{\cal N}^W = \dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial \phi}-\sin\theta
\dfrac{\partial Z}{\partial \theta}\right)+{\cal A}_r\,.\]</div>
<p>The <span class="math notranslate nohighlight">\(\theta\)</span>-derivative entering the radial component of the Coriolis force is thus the
operator <span class="math notranslate nohighlight">\(\vartheta_2\)</span> defined in <a class="reference internal" href="#equation-eqoptheta1">(12)</a>. Using the recurrence
relation, one thus finally gets in spherical harmonic space:</p>
<div class="math notranslate nohighlight" id="equation-eqnlw">
<span class="eqno">(32)<a class="headerlink" href="#equation-eqnlw" title="Link to this equation">¶</a></span>\[\boxed{
{\cal N}^W_{\ell m} = \dfrac{2}{r}\left[i m \dfrac{\partial W_\ell^m}{\partial r}-(\ell-1)c_\ell^m
Z_{\ell-1}^m+(\ell+2)c_{\ell+1}^m Z_{\ell+1}^m\right]
+{{\cal A}_r}_\ell^m\, .
}\]</div>
<p>To get this expression, we need to first compute <span class="math notranslate nohighlight">\({\cal A}_r\)</span> in the physical space. This
term is computed in the subroutine <code class="xref f f-subr docutils literal notranslate"><span class="pre">get_nl</span></code> in
the module <code class="xref f f-mod docutils literal notranslate"><span class="pre">grid_space_arrays_mod</span></code>. <span class="math notranslate nohighlight">\({\cal A}_r\)</span> is then transformed to the
spectral space by using a Legendre and a Fourier transform to produce <span class="math notranslate nohighlight">\({{\cal A}_r}_\ell^m\)</span>.</p>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The final calculations of <a class="reference internal" href="#equation-eqnlw">(32)</a> are done in the subroutine
<code class="xref f f-subr docutils literal notranslate"><span class="pre">get_td</span></code>.</p>
</div>
</section>
<section id="nonlinear-terms-entering-the-equation-for-z">
<span id="secnonlinearz"></span><h3>Nonlinear terms entering the equation for <span class="math notranslate nohighlight">\(Z\)</span><a class="headerlink" href="#nonlinear-terms-entering-the-equation-for-z" title="Link to this heading">¶</a></h3>
<p>The nonlinear term <span class="math notranslate nohighlight">\({\cal N}^Z\)</span> that enters the equation for the toroidal potential
<a class="reference internal" href="#equation-eqspecz">(21)</a> contains the radial component of the curl of the advection and Coriolis force.
The Coriolis force can be rewritten as a function of <span class="math notranslate nohighlight">\(W\)</span> and <span class="math notranslate nohighlight">\(Z\)</span>:</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{e_r}\cdot\vec{\nabla}\times\left[(2\tilde{\rho}\vec{u})\times
\vec{e_z}\right] & =2\vec{e_r}\cdot\left[(\vec{e_z}\cdot\vec{\nabla})(\tilde{\rho}
\vec{u})\right], \\
& = 2\left[\cos\theta\dfrac{\partial (\tilde{\rho} u_r)}{\partial r}
-\dfrac{\sin\theta}{r}\dfrac{\partial (\tilde{\rho}
u_r)}{\partial \theta}+\dfrac{\tilde{\rho} u_\theta\sin\theta}{r}\right], \\
& = 2\left[-\cos\theta\dfrac{\partial}{\partial r}(\nabla_H^2 W)+
\dfrac{\sin\theta}{r}\dfrac{\partial}{\partial \theta}(\nabla_H^2
W)+\dfrac{\sin\theta}{r^2}\dfrac{\partial^2 W}{\partial r\partial \theta}+
\dfrac{1}{r^2}\dfrac{\partial Z}{\partial \phi}\right].
\end{aligned}\end{split}\]</div>
<p>Using the <span class="math notranslate nohighlight">\(\vartheta\)</span> operators defined in <a class="reference internal" href="#equation-eqoptheta1">(12)</a>-<a class="reference internal" href="#equation-eqoptheta4">(15)</a> then
allows to rewrite the Coriolis force in the following way:</p>
<div class="math notranslate nohighlight" id="equation-eqcorznl">
<span class="eqno">(33)<a class="headerlink" href="#equation-eqcorznl" title="Link to this equation">¶</a></span>\[\vec{e_r}\cdot\vec{\nabla}\times\left[(2\tilde{\rho}\vec{u})\times
\vec{e_z}\right]=\dfrac{2}{r^2}\left(\vartheta_3\,\dfrac{\partial W}{\partial r}
-\dfrac{1}{r}\,\vartheta_4\,W+ \dfrac{\partial Z}{\partial \phi} \right)\,.\]</div>
<p>The contributions of nonlinear advection and Lorentz forces that enter the equation
for the toroidal potential are written this way:</p>
<div class="math notranslate nohighlight">
\[\dfrac{1}{r\sin\theta}\left[
\dfrac{\partial (\sin\theta{\cal A}_\phi)}{\partial \theta} -
\dfrac{\partial {\cal A}_\theta}{
\partial\phi}\right]\,.\]</div>
<p>To make use of the recurrence relations <a class="reference internal" href="#equation-eqoptheta1">(12)</a>-<a class="reference internal" href="#equation-eqoptheta4">(15)</a>, the actual
strategy is to follow the following steps:</p>
<ol class="arabic">
<li><p>Compute the quantities <span class="math notranslate nohighlight">\({\cal A}_\phi/r\sin\theta\)</span>
and <span class="math notranslate nohighlight">\({\cal A}_\theta/r\sin\theta\)</span> in the physical space. In the code, this step
is computed in the subroutine <code class="xref f f-subr docutils literal notranslate"><span class="pre">get_nl</span></code> in
the module <code class="xref f f-mod docutils literal notranslate"><span class="pre">grid_space_arrays_mod</span></code>.</p></li>
<li><p>Transform <span class="math notranslate nohighlight">\({\cal A}_\phi/r\sin\theta\)</span> and <span class="math notranslate nohighlight">\({\cal A}_\theta/r\sin\theta\)</span> to
the spectral space (thanks to a Legendre and a Fourier transform). In MagIC, this step
is computed in the modules <code class="xref f f-mod docutils literal notranslate"><span class="pre">legendre_grid_to_spec</span></code> and <a class="reference internal" href="apiFortran/legFourierChebAlgebra.html#f/fft" title="f/fft: This module contains the native subroutines used to compute FFT's. They are based on the FFT99 package from Temperton: http://www.cesm.ucar.edu/models/cesm1.2/cesm/cesmBbrowser/html_code/cam/fft99.F90.html"><code class="xref f f-mod docutils literal notranslate"><span class="pre">fft</span></code></a>. After
this step <span class="math notranslate nohighlight">\({{\cal A}t}_{\ell}^m\)</span> and <span class="math notranslate nohighlight">\({{\cal A}p}_{\ell}^m\)</span> are defined.</p></li>
<li><p>Calculate the colatitude and theta derivatives using the recurrence relations:</p>
<div class="math notranslate nohighlight" id="equation-eqadvznl">
<span class="eqno">(34)<a class="headerlink" href="#equation-eqadvznl" title="Link to this equation">¶</a></span>\[\vartheta_1\,{{\cal A}p}_{\ell}^m-\dfrac{\partial {{\cal A}t}_{\ell}^m}{\partial \phi}\,.\]</div>
</li>
</ol>
<p>Using <a class="reference internal" href="#equation-eqcorznl">(33)</a> and <a class="reference internal" href="#equation-eqadvznl">(34)</a>, one thus finally gets</p>
<div class="math notranslate nohighlight" id="equation-eqnlz">
<span class="eqno">(35)<a class="headerlink" href="#equation-eqnlz" title="Link to this equation">¶</a></span>\[\begin{split}\boxed{
\begin{aligned}
{\cal N}^Z_{\ell m} = & \dfrac{2}{r^2}\left[(\ell-1)(\ell+1)\,c_\ell^m\,
\dfrac{\partial W_{\ell-1}^m}{\partial r}+\ell(\ell+2)\,c_{\ell+1}^m\,
\dfrac{\partial W_{\ell+1}^m}{\partial r} \right. \\
& \left. -\dfrac{\ell(\ell-1)(\ell+1)}{r}\,c_\ell^m\,W_{\ell-1}^m+
\dfrac{\ell(\ell+1)(\ell+2)}{r}\,c_{\ell+1}^m\,W_{\ell+1}^m+
im\,Z_\ell^m\right] \\
& + (\ell+1)\,c_\ell^m\,{{\cal A}p}_{\ell-1}^m-
\ell\,c_{\ell+1}^m\,{{\cal A}p}_{\ell+1}^m
-im\,{{\cal A}t}_{\ell}^m\,.
\end{aligned}
}\end{split}\]</div>
<div class="admonition seealso">
<p class="admonition-title">See also</p>
<p>The final calculations of <a class="reference internal" href="#equation-eqnlz">(35)</a> are done in the subroutine
<code class="xref f f-subr docutils literal notranslate"><span class="pre">get_td</span></code>.</p>
</div>
</section>
<section id="nonlinear-terms-entering-the-equation-for-p">
<span id="secnonlinearp"></span><h3>Nonlinear terms entering the equation for <span class="math notranslate nohighlight">\(P\)</span><a class="headerlink" href="#nonlinear-terms-entering-the-equation-for-p" title="Link to this heading">¶</a></h3>
<p>The nonlinear term <span class="math notranslate nohighlight">\({\cal N}^P\)</span> that enters the equation for the pressure
<a class="reference internal" href="#equation-eqspecp">(22)</a> contains the horizontal divergence of the advection and Coriolis force.
The Coriolis force can be rewritten as a function of <span class="math notranslate nohighlight">\(W\)</span> and <span class="math notranslate nohighlight">\(Z\)</span>:</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{aligned}
\vec{\nabla}_H\cdot\left[(2\tilde{\rho}\vec{u})\times
\vec{e_z}\right] & =2\vec{e_z}\cdot\left[\vec{\nabla}\times(\tilde{\rho}
\vec{u})\right] -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right)
\left[\vec{e_r}\cdot(2\tilde{\rho}\vec{u}\times\vec{e_z})\right],\\
& = -2\cos\theta\,\nabla_H^2 Z-2\sin\theta\left[-\dfrac{1}{r\sin\theta}
\dfrac{\partial}{\partial\phi}\left(
\dfrac{\partial^2}{\partial r^2}+\nabla_H^2 \right) W +
\dfrac{1}{r}\dfrac{\partial^2 Z}{\partial r\partial\theta}\right]
\\
& \phantom{=\cos\theta} -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right)
\left[2\sin\theta\tilde{\rho}u_\phi\right], \\
& = 2\left[\dfrac{1}{r}\left(\nabla_H^2+\dfrac{\partial^2}{\partial r^2}\right)
\dfrac{\partial W}{\partial \phi}-\cos\theta\nabla_H^2 Z -\dfrac{\sin\theta}{r}
\dfrac{\partial^2 Z}{\partial r \partial \theta}\right] \\
& \phantom{=\cos\theta} -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right)
\left[\dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial\phi}-\sin\theta