SphericalFluid

English | 日本語

A real-time FMX demo application that solves two-dimensional incompressible fluid flow on the sphere — the barotropic vorticity equation — by the spherical-harmonic spectral transform method, using the LUX.SphericalHarmonics library, and renders the vorticity field as a colored relief on a 3D globe.

利用ライブラリ

• LUX ：The LUXOPHIA standard library providing basic mathematical types such as vectors, matrices, and complex numbers.
• LUX.FMX.Graphics.D3 ：A 3D graphics helper library built on the FireMonkey framework.
• LUX.Sphere ：A spherical geometry library for shapes, curves, and barycenter operations on the spheres $S^2$ and $S^3$.
• LUX.SphericalHarmonics ：A library for spherical harmonics and the associated Legendre functions.

1. Overview / Features

• Spectral transform method: the vorticity field is expanded in spherical harmonics with triangular truncation $T_N$ ($N$ = 21 / 42 / 85 / 170, selectable at run time). The inverse Laplacian is exact in spectral space.
• Alias-free transform grid: the quadratic nonlinear term is evaluated on a longitude–Gauss-latitude grid with $I \ge 3N+1$ and $J \ge (3N+1)/2$, so no aliasing error enters the retained modes [3].
• Diffusion-free: the model contains no viscous or filter term (an optional hyperviscosity filter exists, default off). The flux form of the nonlinear term, evaluated by integration by parts, conserves both energy and enstrophy — numerical diffusion is zero, leaving only the truncation error of the RK4 time integration.
• Rendering: vorticity is displayed as a radial relief of the sphere's vertices (ridge = positive, trough = negative), colored blue (trough) → white (zero) → red (ridge). The display mesh is always 256 × 128, resynthesized spectrally from the simulation coefficients (band-limited interpolation, exact), so even low resolutions look smooth. Surface normals are computed from the analytic derivatives $\partial\zeta/\partial\theta$, $\partial\zeta/\partial\lambda$ of the spherical-harmonic expansion; the relief uses a $\tanh$ mapping (continuously differentiable). The camera orbits while the lights stay fixed in world space.
• The info panel continuously shows the relative drift of energy and enstrophy, so the conservation (diffusion-free) property can be verified in real time.
• A companion application that visualizes the individual spherical-harmonic basis functions is available as SphericalHarmonics [9].

1.1 Verification (console tests)

Item	Result
Transform round-trip error (T170)	~1e-13
Rossby-wave phase speed vs. analytic solution	agrees to machine precision
Energy / enstrophy conservation (T85, 50 steps)	~1e-12
Speed (one RK4 step, Ryzen-class CPU)	T85 ≈ 6.5 ms / T170 ≈ 21 ms

2. Mathematical Background

2.1 Governing equation

The model describes a non-divergent two-dimensional fluid covering the whole sphere (no land), with the vorticity $\zeta$ as the single prognostic variable — the barotropic vorticity equation [1] in flux form:

$$\frac{\partial \zeta}{\partial t} = -\nabla\cdot\bigl(\mathbf{v}\,(\zeta + f)\bigr), \qquad f = 2\Omega\mu \tag{1}$$ $$\nabla^{2}\Psi = \zeta, \qquad \mathbf{v} = \mathbf{k}\times\nabla\Psi \tag{2}$$

where $\Psi$ is the streamfunction, $\mathbf{v}$ the horizontal velocity, $\Omega$ the rotation rate, $\lambda$ the longitude, $\varphi$ the latitude, $\mu = \sin\varphi = \cos\theta$, and $\mathbf{k}$ the local vertical unit vector. The sphere radius is normalized to $1$.

2.2 Spherical-harmonic expansion

The vorticity is expanded in spherical harmonics [4] with triangular truncation $T_N$:

$$\zeta(\lambda,\mu,t) = \sum_{m=-N}^{N}\ \sum_{n=|m|}^{N} \zeta_{n}^{m}(t)\,Y_{n}^{m}(\lambda,\mu), \qquad Y_{n}^{m}(\lambda,\mu) = \tilde{P}_{n}^{m}(\mu)\,e^{im\lambda} \tag{3}$$

where $\tilde{P}_n^m$ are the normalized associated Legendre functions [5]

$$\tilde{P}_{n}^{m}(\mu) = \sqrt{\frac{2n+1}{2}\,\frac{(n-m)!}{(n+m)!}}\;P_{n}^{m}(\mu), \qquad \int_{-1}^{1} \tilde{P}_{n}^{m}(\mu)\,\tilde{P}_{n'}^{m}(\mu)\,d\mu = \delta_{nn'} \tag{4}$$

evaluated by the four-term recurrence engine TNALFsTerm4 of the LUX.SphericalHarmonics library. Since $\zeta$ is real, $\zeta_n^{-m} = (\zeta_n^m)^{*}$ and only $m \ge 0$ is stored.

The spherical harmonics are eigenfunctions of the Laplacian on $S^2$:

$$\nabla^{2} Y_{n}^{m} = -\,n(n+1)\,Y_{n}^{m} \tag{5}$$

so the Laplacian is diagonal in spectral space and the inversion $(2)$ is exact:

$$\Psi_{n}^{m} = -\frac{\zeta_{n}^{m}}{n(n+1)} \qquad (n \ge 1), \qquad \Psi_{0}^{0} = 0 \tag{6}$$

as implemented in TSHFluid.Tendency.

2.3 Spectral transform method

One evaluation of the right-hand side of $(1)$ (spectral transform method [2][8]) proceeds as:

1. Synthesis (spectral → grid): $\zeta$ and the pseudo-velocities

$$U = u\cos\varphi = -(1-\mu^{2})\frac{\partial\Psi}{\partial\mu}, \qquad V = v\cos\varphi = \frac{\partial\Psi}{\partial\lambda} \tag{7}$$

are synthesized on the transform grid (Synthesize, SynthesizeUV) — equally spaced longitudes with FFT, Gauss–Legendre latitudes $\mu_j$ (nodes of $P_J$, computed by Newton iteration).

2. Fluxes on the grid: $A = U(\zeta+f)$, $B = V(\zeta+f)$; the maximum wind speed is recorded for the CFL condition.

3. Analysis of the flux divergence (grid → spectral, AnalyzeFluxDiv): using integration by parts against $\tilde{P}_n^m$, the divergence

$$\nabla\cdot\bigl(\mathbf{v}\,\zeta_a\bigr) = \frac{1}{1-\mu^{2}}\frac{\partial A}{\partial\lambda} + \frac{\partial B}{\partial\mu} \tag{8}$$

is obtained exactly in spectral space; the latitude integrals use Gauss–Legendre quadrature [6], which is exact for band-limited integrands.

Because the nonlinearity is quadratic, choosing the grid size

$$I \ge 3N+1, \qquad J \ge \frac{3N+1}{2} \tag{9}$$

(the 3/2 de-aliasing rule [3], enforced by TSHTransform.FitGrid) removes all aliasing errors from the retained modes; combined with the flux form $(8)$ this makes the semi-discrete scheme conserve energy and enstrophy exactly.

2.4 Time integration and conserved quantities

Time stepping uses the classical fourth-order Runge–Kutta scheme [7] (TSHFluid.Step):

$$\zeta^{t+\Delta t} = \zeta^{t} + \frac{\Delta t}{6}\bigl(k_1 + 2k_2 + 2k_3 + k_4\bigr) \tag{10}$$

with the step size limited by the CFL-type condition (SuggestDt)

$$\Delta t = \frac{C}{N\,\max(V_{\max},\,0.1)}, \qquad C = 0.5 \tag{11}$$

The conserved diagnostics, evaluated directly from the spectral coefficients, are the kinetic energy and the enstrophy:

$$E = \tfrac{1}{2}\int_{S^2} |\nabla\Psi|^{2}\,dS = \pi \sum_{n,m} (2-\delta_{m0})\,\frac{|\zeta_{n}^{m}|^{2}}{n(n+1)}, \qquad Z = \tfrac{1}{2}\int_{S^2} \zeta^{2}\,dS = \pi \sum_{n,m} (2-\delta_{m0})\,|\zeta_{n}^{m}|^{2} \tag{12}$$

The optional hyperviscosity filter (default off = fully diffusion-free) relaxes the highest wavenumbers after each step:

$$\zeta_{n}^{m} \leftarrow \zeta_{n}^{m}\, \exp\!\left[-50\,\Delta t\left(\frac{n(n+1)}{N(N+1)}\right)^{8}\right] \tag{13}$$

The initial condition (InitRandom) puts Gaussian random vorticity into the large scales $3 \le n \le 8$, normalized to RMS wind speed $1$.

2.5 Rendering mathematics

The display synthesizer (TSHTransform with Loose mode, requiring only the sampling condition $I \ge 2N+2$) resynthesizes the coefficients on a fixed 256 × 128 mesh. The surface is the radial relief

$$\mathbf{P}(\theta,\lambda) = f\,\hat{\mathbf{r}}, \qquad f = R\Bigl(1 + a\tanh\frac{\zeta}{s}\Bigr) \tag{14}$$

and the exact surface normal, using the analytic derivatives $\partial\zeta/\partial\theta$ and $\partial\zeta/\partial\lambda$ synthesized from the coefficients (SynthesizeGrad), is

$$\mathbf{N} \propto f\,\hat{\mathbf{r}} - \frac{\partial f}{\partial\theta}\,\hat{\boldsymbol{\theta}} - \frac{1}{\sin\theta}\frac{\partial f}{\partial\lambda}\,\hat{\boldsymbol{\lambda}} \tag{15}$$

as implemented in TFluidSphere3D.MakeGeometry.

3. Architecture

[Ownership] main form → scene, solver, display synthesizer

・TFormMain (Main.pas)
  ┣・TWorld3D                        ･･･ 3D scene (LUX.FMX.Graphics.D3)
  ┃  ┣・TCamera3D
  ┃  ┗・TLight3D ×3
  ┣・TFluidSphere3D                  ･･･ relief mesh (FluidSphere.pas)
  ┃  ┗・SetGrid / SetField          ･･･ vertices, analytic normals, palette
  ┣・TSHFluid                        ･･･ BVE solver (LUX.SH.Fluid)
  ┃  ┣・Tendency / Step(RK4) / Energy / Enstrophy / SuggestDt
  ┃  ┗・TSHTransform                ･･･ SH transform (LUX.SH.Transform)
  ┃     ┣・Synthesize / SynthesizeUV / SynthesizeGrad
  ┃     ┣・Analyze / AnalyzeFluxDiv ･･･ (Gauss-Legendre grid + FFT, TTask)
  ┃     ┗・TNALFsTerm4              ･･･ normalized ALFs (LUX.NALFs.Term4)
  ┗・TSHTransform (Loose)            ･･･ display synthesizer, 256×128

[Inheritance] relief mesh shape

・TF3DShaper                          ･･･ (LUX.FMX.Graphics.D3)
  ┗・TFluidSphere3D                  ･･･ (FluidSphere.pas)

・SphericalFluid/
  ┣・SphericalFluid.dpr / .dproj ･･･ project (FMX application)
  ┣・Main.pas / Main.fmx         ･･･ main form: UI, 3D scene, animation loop
  ┣・FluidSphere.pas             ･･･ TFluidSphere3D: vorticity relief sphere
  ┗・_LIBRARY/LUXOPHIA/
     ┣・LUX.SH.Transform.pas     ･･･ SH transform engine (under development)
     ┣・LUX.SH.Fluid.pas         ･･･ BVE solver (under development)
     ┣・LUX/                     ･･･ base math (vectors, matrices, complex)
     ┣・LUX.FMX.Graphics.D3/     ･･･ FMX 3D helpers (TWorld3D, TF3DShaper)
     ┣・LUX.Sphere/              ･･･ sphere shapes
     ┗・LUX.SphericalHarmonics/  ･･･ ALFs/NALFs tables

Library subtrees on GitHub: LUX ・ LUX.FMX.Graphics.D3 ・ LUX.Sphere ・ LUX.SphericalHarmonics

The LUX.SH.Transform engine parallelizes the latitude loops with TTask (one chunk per core) and pre-allocates all scratch buffers, so no memory is allocated inside the workers.

4. Usage / Controls

UI	Function
Resolution	T21 / T42 / T85 (default) / T170 (triangular truncation / transform grid)
Reset	restart from a random large-scale ($3 \le n \le 8$) vorticity field
Rotation Rate $\Omega$	Coriolis term ($\beta$-effect → Rossby waves, zonal jets)
Relief Amplitude	height of the ridges and troughs
Time Scale	simulated time per real-time second
Hyperviscosity	high-wavenumber relaxation $(13)$ (default off = fully diffusion-free)
Left mouse drag	orbit the camera

The info panel continuously displays FPS, simulated time, maximum wind speed, and the relative drift of energy and enstrophy $(12)$.

5. Building

Open SphericalFluid.dproj in RAD Studio (Delphi 13) and build. Target platforms defined in the project: Win32 and Win64 (Win64 is the primary target).

Command line:

dcc64 -B -$O+ -U"%BDS%\lib\win64\release" -R"%BDS%\lib\win64\release" -NSSystem;System.Win;Winapi -EWin64\Release -NUWin64\Release SphericalFluid.dpr

6. References

1. Barotropic vorticity equation — Wikipedia
2. Bourke, W. (1972): An Efficient, One-Level, Primitive-Equation Spectral Model. Monthly Weather Review, 100, 683–689.
3. Orszag, S. A. (1970): Transform Method for the Calculation of Vector-Coupled Sums: Application to the Spectral Form of the Vorticity Equation. Journal of the Atmospheric Sciences, 27, 890–895.
4. Spherical harmonics — Wikipedia
5. Associated Legendre polynomials — Wikipedia
6. Gauss–Legendre quadrature — Wikipedia
7. Runge–Kutta methods — Wikipedia
8. Hack, J. J. and Jakob, R. (1992): Description of a Global Shallow Water Model Based on the Spectral Transform Method. NCAR Technical Note NCAR/TN-343+STR.
9. SphericalHarmonics — 3D viewer application for spherical harmonics (LUXOPHIA).

💖 Embarcadero Delphi

Integrated Development Environment (IDE) for Creating Native Cross-Platform Apps.

Free Download: Delphi Community Edition