Model‐Physics‐and‐Mathematics - TSOR666/Atmospheric-and-radivative-transfer-suite GitHub Wiki
Model Physics and Mathematics
Table of Contents
- Overview
- MDNO v5.3: Multi-scale Diffusion Neural Operator
- RTNO v4.3: Radiative Transfer Neural Operator
- STRATUS: Modular Radiative Transfer Framework
- Numerical Methods
- Conservation Laws and Constraints
- Physical Constants
- Computational Complexity
- References
- Notation Summary
Overview
This repository implements three interconnected models for atmospheric dynamics and radiative transfer:
- MDNO v5.3: A multi-scale neural operator for atmospheric dynamics based on the Boltzmann equation and Navier-Stokes equations
- RTNO v4.3: A neural operator for radiative transfer with full polarization support
- STRATUS: A modular framework for volumetric ray marching, Monte Carlo transport, and physics-informed neural networks
All models enforce physical conservation laws and support both deterministic and stochastic solvers.
MDNO v5.3: Multi-scale Diffusion Neural Operator
1. Mathematical Foundation
MDNO operates on three spatial scales (micro, meso, macro) and couples kinetic theory with continuum fluid dynamics.
1.1 Boltzmann Equation (Micro-scale)
The fundamental kinetic equation governing the evolution of the particle distribution function $f(\mathbf{x}, \mathbf{v}, t)$:
$$\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f + \mathbf{F} \cdot \nabla_{\mathbf{v}} f = \mathcal{C}(f)$$
where:
- $f(\mathbf{x}, \mathbf{v}, t)$ is the distribution function in phase space
- $\mathbf{v} = (v_x, v_y, v_z)$ is the velocity vector
- $\mathbf{F}$ represents external forces
- $\mathcal{C}(f)$ is the collision operator
Discrete Implementation:
- Velocity space discretized on grid: $v \in [-v_{\max}, v_{\max}]^3$ with resolution $N_v = 8$
- Spatial advection computed using spectral methods with FFT
- Collision operator approximated by neural network with BGK relaxation:
$$\mathcal{C}(f) = -\frac{1}{\tau}(f - f^{\text{eq}})$$
where $f^{\text{eq}}$ is the Maxwell-Boltzmann equilibrium distribution:
$$f^{\text{eq}}(v) = n \left(\frac{m}{2\pi k_B T}\right)^{3/2} \exp\left(-\frac{m|\mathbf{v} - \mathbf{u}|^2}{2k_B T}\right)$$
1.2 Complete Primitive Equations (Meso-scale)
The meso-scale dynamics are governed by the full compressible Navier-Stokes equations with moisture and thermodynamics:
Momentum Equation (x-component):
$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + fv + \nu\nabla^2 u + F_u$$
Momentum Equation (y-component):
$$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial y} - fu + \nu\nabla^2 v + F_v$$
Momentum Equation (z-component):
$$\frac{\partial w}{\partial t} + u\frac{\partial w}{\partial x} + v\frac{\partial w}{\partial y} + w\frac{\partial w}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial z} - g + \nu\nabla^2 w + F_w$$
Thermodynamic Equation:
$$\frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} + w\frac{\partial T}{\partial z} = \kappa\nabla^2 T + \frac{\kappa}{c_p}T\frac{w}{p}\frac{\partial p}{\partial z} + \frac{Q}{c_p}$$
Moisture Equation:
$$\frac{\partial q}{\partial t} + u\frac{\partial q}{\partial x} + v\frac{\partial q}{\partial y} + w\frac{\partial q}{\partial z} = S$$
Continuity Equation:
$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$$
Equation of State:
$$p = \rho R_d T$$
Physical Parameters:
- $(u, v, w)$ are velocity components (m/s)
- $T$ is temperature (K)
- $p$ is pressure (Pa)
- $\rho$ is density (kg/m³)
- $q$ is specific humidity (kg/kg)
- $f$ is Coriolis parameter: $f = 2\Omega\sin\phi$ where $\Omega = 7.2921 \times 10^{-5}$ rad/s
- $\nu$ is kinematic viscosity (m²/s)
- $\kappa$ is thermal diffusivity (m²/s)
- $g = 9.80665$ m/s² (gravitational acceleration)
- $R_d = 287.058$ J/(kg·K) (dry air gas constant)
- $c_p = 1005.0$ J/(kg·K) (specific heat at constant pressure)
1.3 Hamiltonian Dynamics (Macro-scale)
For the large-scale flow, we enforce energy conservation through Hamiltonian mechanics:
$$\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial q_i}$$
The Hamiltonian is:
$$H = T + V = \frac{1}{2}\sum_i p_i^2 + \frac{1}{2}\sum_i q_i^2$$
where $q_i$ represents generalized coordinates and $p_i$ represents conjugate momenta.
Symplectic Integration: We use a symplectic integrator to preserve the Hamiltonian structure:
$$q^{n+1} = q^n + \Delta t \frac{\partial H}{\partial p}\bigg|_{p^{n+1/2}}$$
$$p^{n+1} = p^n - \Delta t \frac{\partial H}{\partial q}\bigg|_{q^{n+1}}$$
2. Turbulence Parameterization
2.1 Large Eddy Simulation (LES) Closure
The Smagorinsky subgrid-scale model:
$$\nu_t = (C_s \Delta)^2 |S|$$
where:
- $C_s = 0.17$ is the Smagorinsky constant
- $\Delta = (dx \cdot dy \cdot dz)^{1/3}$ is the filter width
- $|S| = \sqrt{2S_{ij}S_{ij}}$ is the strain rate magnitude
The strain rate tensor:
$$S_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)$$
Subgrid-scale stress tensor:
$$\tau_{ij} = -2\nu_t S_{ij}$$
2.2 k-ε Turbulence Model (RANS)
For Reynolds-Averaged Navier-Stokes:
Turbulent kinetic energy equation:
$$\frac{\partial k}{\partial t} + u_j\frac{\partial k}{\partial x_j} = P_k - \epsilon + \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_j}\right]$$
Dissipation rate equation:
$$\frac{\partial \epsilon}{\partial t} + u_j\frac{\partial \epsilon}{\partial x_j} = C_{1\epsilon}\frac{\epsilon}{k}P_k - C_{2\epsilon}\frac{\epsilon^2}{k} + \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_\epsilon}\right)\frac{\partial \epsilon}{\partial x_j}\right]$$
where:
- $k$ is turbulent kinetic energy
- $\epsilon$ is dissipation rate
- $\nu_t = C_\mu k^2/\epsilon$ is eddy viscosity
- Standard constants: $C_\mu = 0.09$, $C_{1\epsilon} = 1.44$, $C_{2\epsilon} = 1.92$, $\sigma_k = 1.0$, $\sigma_\epsilon = 1.3$
3. Fourier Neural Operator (FNO) Architecture
The core spatial operator uses spectral methods:
3.1 3D Fourier Layer
For input $u(\mathbf{x})$, the Fourier layer computes:
Forward Fourier Transform:
$(\mathcal{F}u)(\mathbf{k})=\int_{\Omega} u(\mathbf{x}),e^{-2\pi i,\mathbf{k}\cdot\mathbf{x}},\mathrm{d}\mathbf{x}$
Spectral Convolution:
$(\mathcal{K}u)_{\mathrm{spectral}}(\mathbf{x})=\mathcal{F}^{-1}!\bigl(W(\mathbf{k}),\mathcal{F}{u}(\mathbf{k})\bigr)(\mathbf{x})$
Total Operator:
$$(v)(\mathbf{x}) = \sigma(W_{\text{local}}u + (\mathcal{K}u)_{\text{spectral}})$$
where:
- $W(\mathbf{k})$ are learnable spectral weights (complex-valued)
- $W_{\text{local}}$ is a 1×1 convolution for local features
- $\sigma$ is a nonlinear activation (SiLU)
- Mode truncation: $(k_x, k_y, k_z) \in [0, 16] \times [0, 16] \times [0, 8]$
3.2 Spectral Anti-aliasing
To prevent aliasing, we apply dealiasing with the 2/3 rule:
$$k_{\max} = \frac{2}{3}k_{\text{Nyquist}}$$
For adaptive derivatives in advection terms, we compute spatial derivatives in Fourier space:
$$\frac{\partial u}{\partial x} = \mathcal{F}^{-1}[2\pi i k_x \mathcal{F}[u]]$$
4. Cloud Microphysics
Phase transitions are modeled using:
Condensation/Evaporation:
$$\frac{dq_v}{dt} = -C(q_v - q_{vs})$$
where:
- $q_v$ is water vapor mixing ratio
- $q_{vs}$ is saturation mixing ratio
- $C$ is condensation rate
Saturation mixing ratio (Clausius-Clapeyron equation):
$$q_{vs} = \frac{\epsilon e_s}{p - e_s}$$
$$e_s = e_0 \exp\left[\frac{L_v}{R_v}\left(\frac{1}{T_0} - \frac{1}{T}\right)\right]$$
where:
- $L_v = 2.5 \times 10^6$ J/kg is latent heat of vaporization
- $R_v = 461.495$ J/(kg·K) is water vapor gas constant
- $e_0 = 611$ Pa is reference saturation vapor pressure at $T_0 = 273.15$ K
5. Conservation Enforcement
5.1 Moment Conservation
For the Boltzmann equation, macroscopic quantities are moments of $f$:
Density (zeroth moment):
$$\rho = \int f(\mathbf{v}) d\mathbf{v}$$
Momentum (first moment):
$$\rho \mathbf{u} = \int \mathbf{v} f(\mathbf{v}) d\mathbf{v}$$
Energy (second moment):
$$\frac{3}{2}\rho k_B T = \int \frac{1}{2}m|\mathbf{v} - \mathbf{u}|^2 f(\mathbf{v}) d\mathbf{v}$$
Entropic Projection: To preserve moments while maximizing entropy:
$$ \begin{aligned} f^{*}=\arg\min_{f}\ & -\int f\ln f,\mathrm{d}\mathbf{v} \ \text{s.t. } & \int f,\mathrm{d}\mathbf{v}=\rho,\qquad \int \mathbf{v}f,\mathrm{d}\mathbf{v}=\rho\mathbf{u} \end{aligned} $$
5.2 CFL Condition
For numerical stability:
$$\text{CFL} = \max\left(\frac{|u|\Delta t}{\Delta x}, \frac{|v|\Delta t}{\Delta y}, \frac{|w|\Delta t}{\Delta z}\right) < 0.5$$
Adaptive time-stepping adjusts $\Delta t$ to maintain CFL < 0.5.
RTNO v4.3: Radiative Transfer Neural Operator
1. Radiative Transfer Equation (RTE)
The fundamental equation for radiative transfer in a scattering and absorbing medium:
Optical Depth Form:
$$\mu\frac{\partial I_\lambda(\tau_\lambda, \mu, \phi)}{\partial \tau_\lambda} = I_\lambda(\tau_\lambda, \mu, \phi) - S_\lambda(\tau_\lambda, \mu, \phi)$$
Geometric Form:
$$\frac{dI_\lambda}{ds} = -\kappa_\lambda I_\lambda + j_\lambda$$
where:
- $I_\lambda$ is spectral radiance (W/(m²·sr·nm))
- $\tau_\lambda = \int_0^s \kappa_\lambda ds'$ is optical depth
- $\mu = \cos\theta$ is the cosine of zenith angle
- $\phi$ is azimuthal angle
- $\kappa_\lambda$ is extinction coefficient (m⁻¹)
- $j_\lambda$ is source function
- $S_\lambda$ is total source term including scattering
2. Scattering Physics
2.1 Phase Function
The general scattering source term:
$$S_\lambda^{\text{scat}} = \frac{\omega_\lambda}{4\pi}\int_{4\pi} P(\Theta) I_\lambda(\mu', \phi') d\Omega'$$
where:
- $\omega_\lambda = \sigma_s/(\sigma_s + \sigma_a)$ is single scattering albedo
- $P(\Theta)$ is the phase function
- $\Theta$ is scattering angle
The scattering angle satisfies:
$$\cos\Theta = \mu\mu' + \sqrt{1-\mu^2}\sqrt{1-\mu'^2}\cos(\phi - \phi')$$
2.2 Rayleigh Scattering
For molecular scattering (particle size << wavelength):
Cross Section:
$$\sigma_R(\lambda) = \frac{8\pi^3(n^2-1)^2}{3N\lambda^4}$$
Rayleigh Phase Function:
$$P_R(\cos\Theta) = \frac{3}{4}(1 + \cos^2\Theta)$$
Wavelength Dependence:
$$\sigma_R \propto \lambda^{-4}$$
2.3 Mie Scattering
For aerosols and cloud droplets (particle size ≈ wavelength):
Size Parameter:
$$x = \frac{2\pi r}{\lambda}$$
Scattering Amplitudes:
The Mie solution expands in terms of vector spherical harmonics. The scattering amplitudes are:
$$S_1(\theta) = \sum_{n=1}^{N_{\max}} \frac{2n+1}{n(n+1)}(a_n \pi_n(\cos\theta) + b_n \tau_n(\cos\theta))$$
$$S_2(\theta) = \sum_{n=1}^{N_{\max}} \frac{2n+1}{n(n+1)}(a_n \tau_n(\cos\theta) + b_n \pi_n(\cos\theta))$$
where $a_n$ and $b_n$ are Mie coefficients computed from Riccati-Bessel functions.
Henyey-Greenstein Approximation:
$$P_{HG}(\cos\Theta) = \frac{1 - g^2}{(1 + g^2 - 2g\cos\Theta)^{3/2}}$$
where $g$ is the asymmetry parameter (typically 0.7-0.9 for clouds).
3. Polarization: Stokes Vector Formalism
3.1 Stokes Vector
The state of polarization is described by the Stokes vector I:
$$\mathbf{I} = \begin{pmatrix} I \ Q \ U \ V \end{pmatrix}$$
where:
- $I$ = total intensity
- $Q$ = linearly polarized component (horizontal - vertical)
- $U$ = linearly polarized component (45° - 135°)
- $V$ = circularly polarized component (right - left)
Physical Constraint:
$$I^2 \geq Q^2 + U^2 + V^2$$
Degree of Linear Polarization:
$$p_L = \frac{\sqrt{Q^2 + U^2}}{I}$$
Degree of Circular Polarization:
$$p_C = \frac{V}{I}$$
3.2 Mueller Matrix Calculus
Scattering transformations are represented by 4×4 Mueller matrices:
$$\mathbf{I}{\text{out}} = \mathbf{M}(\Theta) \cdot \mathbf{I}{\text{in}}$$
Rayleigh Mueller Matrix:
The Rayleigh scattering Mueller matrix is:
M_R = (3/4) × [
[1 + cos²Θ, -sin²Θ, 0, 0 ]
[-sin²Θ, 1 + cos²Θ, 0, 0 ]
[0, 0, 2cosΘ, 0 ]
[0, 0, 0, 2cosΘ ]
]
Mie Mueller Matrix:
The Mie scattering Mueller matrix is:
M_Mie = [
[S₁₁, S₁₂, 0, 0 ]
[S₁₂, S₁₁, 0, 0 ]
[0, 0, S₃₃, -S₃₄]
[0, 0, S₃₄, S₃₃ ]
]
where the elements are computed from the Mie scattering amplitudes:
$$S_{11} = \frac{1}{2}(|S_2|^2 + |S_1|^2)$$
$$S_{12} = \frac{1}{2}(|S_2|^2 - |S_1|^2)$$
$$S_{33} = \text{Re}(S_2 S_1^*)$$
$$S_{34} = \text{Im}(S_2 S_1^*)$$
3.3 Vector RTE
The full polarized RTE for Stokes vector:
$$\mu\frac{\partial \mathbf{I}}{\partial \tau} = \mathbf{I} - \frac{\omega}{4\pi}\int_{4\pi} \mathbf{M}(\Theta)\mathbf{I}(\mu', \phi') d\Omega'$$
4. Gas Absorption
4.1 Beer-Lambert Law
For pure absorption:
$$I(s) = I_0 e^{-\int_0^s \kappa_a(\lambda) ds'}$$
4.2 Line-by-Line Absorption
Absorption cross section for spectral line:
$$\sigma(\nu) = S f(\nu - \nu_0)$$
where $S$ is line strength and $f$ is line shape function.
Voigt Profile:
The Voigt profile is a convolution of Doppler (Gaussian) and pressure (Lorentzian) broadening.
4.3 Correlated-k Distribution
For computational efficiency, group absorption coefficients:
$$\overline{T}(\tau) = \int_0^1 e^{-k(g)\tau} dg$$
where $k(g)$ is the absorption coefficient at cumulative probability $g$.
5. Spherical Harmonics Expansion
For angle dependence, expand in spherical harmonics:
$$I(\mu, \phi) = \sum_{l=0}^{L}\sum_{m=-l}^{l} I_l^m Y_l^m(\mu, \phi)$$
where:
$$Y_l^m(\theta, \phi) = \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}} P_l^m(\cos\theta) e^{im\phi}$$
where $P_l^m$ are associated Legendre polynomials.
Discrete Ordinates (P_N method): Angular discretization on Gaussian quadrature points with $N$ streams.
6. Multiple Scattering
6.1 Source Function Iteration
Iterate the source function:
$$S^{(n+1)} = (1-\omega)B + \frac{\omega}{4\pi}\int_{4\pi} P(\Theta) I^{(n)} d\Omega'$$
Continue until convergence:
$$||S^{(n+1)} - S^{(n)}|| < \epsilon$$
6.2 Delta-Eddington Scaling
For highly forward-peaked phase functions, use similarity transformation:
$$\tau^* = (1 - \omega g^2)\tau$$
$$\omega^* = \frac{(1-g^2)\omega}{1 - \omega g^2}$$
$$g^* = \frac{g}{1+g}$$
This prevents excessive scattering orders in forward-scattering media.
7. Boundary Conditions
7.1 Lambertian Surface
Isotropic reflection:
$$I^+(\mu, \phi) = \frac{\rho}{\pi} \int_{\mu' > 0} I^-(\mu', \phi') \mu' d\Omega'$$
where $\rho$ is surface albedo.
7.2 Specular Reflection
$$I^+(\mu, \phi) = I^-(\mu, \phi + \pi)$$
7.3 Ocean Surface (Fresnel Reflection)
The Fresnel reflection coefficient:
$$R(\theta) = \frac{1}{2}\left[\frac{\sin^2(\theta - \theta')}{\sin^2(\theta + \theta')} + \frac{\tan^2(\theta - \theta')}{\tan^2(\theta + \theta')}\right]$$
where $\theta'$ is refraction angle from Snell's law:
$$\sin\theta = n\sin\theta'$$
8. Neural Architecture
8.1 Attention-Based Operator
Multi-head self-attention for angular coupling:
$$\text{Attention}(Q, K, V) = \text{softmax}\left(\frac{QK^T}{\sqrt{d_k}}\right)V$$
where $Q, K, V$ are query, key, value matrices.
8.2 Hybrid Physics-ML Pipeline
- Physical Forward Pass: Solve deterministic RTE using discrete ordinates
- Neural Correction: Learn residuals from high-fidelity reference
- Total Solution: $I_{\text{total}} = I_{\text{physics}} + I_{\text{neural}}$
STRATUS: Modular Radiative Transfer Framework
1. Volumetric Ray Marching
1.1 Ray-Volume Intersection
Ray parameterization:
$$\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$$
where $\mathbf{o}$ is origin and $\mathbf{d}$ is direction.
AABB (Axis-Aligned Bounding Box) Intersection:
$$t_{\min} = \max_i\left(\min\left(\frac{b_i^{\min} - o_i}{d_i}, \frac{b_i^{\max} - o_i}{d_i}\right)\right)$$
$$t_{\max} = \min_i\left(\max\left(\frac{b_i^{\min} - o_i}{d_i}, \frac{b_i^{\max} - o_i}{d_i}\right)\right)$$
1.2 Volumetric Integration
Discretize the RTE along ray path:
$$I(t_{\text{exit}}) = I(t_0) T(t_0, t_{\text{exit}}) + \int_{t_0}^{t_{\text{exit}}} S(t') T(t', t_{\text{exit}}) dt'$$
Transmittance:
$$T(t_1, t_2) = \exp\left(-\int_{t_1}^{t_2} \kappa(t') dt'\right)$$
Discrete Ray Marching:
$$I_{n+1} = I_n e^{-\kappa_n \Delta t} + \frac{S_n}{\kappa_n}(1 - e^{-\kappa_n \Delta t})$$
1.3 Adaptive Step Size
Optical depth control:
$$\Delta t = \min\left(\Delta t_{\text{base}}, \frac{\tau_{\max}}{\kappa + \epsilon}\right)$$
where $\tau_{\max} = 0.1$ ensures accuracy.
1.4 Trilinear Interpolation
For sampling at non-grid points, trilinear interpolation is used. For a point $\mathbf{p}$ between grid cells, the interpolated value is computed using the 8 neighboring grid values weighted by their distances.
2. Monte Carlo Radiative Transfer
2.1 Photon Path Sampling
Free path length: Sample from exponential distribution:
$$P(s) = \kappa e^{-\kappa s}$$
$$s = -\frac{\ln(\xi)}{\kappa}$$
where $\xi \sim U(0,1)$ is a uniform random variable.
Scattering angle: Sample from phase function (Henyey-Greenstein case):
$$\cos\Theta = \frac{1}{2g}\left[1 + g^2 - \left(\frac{1-g^2}{1+g(2\xi-1)}\right)^2\right]$$
For isotropic scattering:
$$\cos\Theta = 2\xi - 1$$
Azimuthal angle:
$$\phi = 2\pi\xi$$
2.2 Russian Roulette
Terminate low-weight photons probabilistically:
$$P_{\text{survive}} = \min\left(1, \frac{w}{w_{\text{threshold}}}\right)$$
If photon survives: $w' = w_{\text{threshold}}$
2.3 Variance Reduction
Importance Sampling: Sample preferentially in regions of high contribution.
Stratified Sampling: Divide phase space into strata and sample uniformly within each.
3. Physics-Informed Neural Networks (PINN)
3.1 PINN Formulation
Minimize combined loss:
$$\mathcal{L} = \mathcal{L}{\text{PDE}} + \lambda{\text{BC}}\mathcal{L}{\text{BC}} + \lambda{\text{IC}}\mathcal{L}{\text{IC}} + \lambda{\text{data}}\mathcal{L}_{\text{data}}$$
PDE Residual:
$$\mathcal{L}{\text{PDE}} = \frac{1}{N}\sum{i=1}^N \left|\mu\frac{\partial I}{\partial z} + \kappa I - S\right|^2$$
Derivatives computed via automatic differentiation:
$$\frac{\partial I}{\partial z} = \frac{\partial \mathcal{N}_\theta(\mathbf{x})}{\partial z}$$
where $\mathcal{N}_\theta$ is the neural network with parameters $\theta$.
3.2 Network Architecture
Multi-layer perceptron with $L$ layers. The network computes:
$$\mathcal{N}\theta(\mathbf{x}, t) = W_L \sigma(W{L-1} \sigma(\cdots \sigma(W_1 \mathbf{x} + b_1) \cdots) + b_{L-1}) + b_L$$
Activation: GELU or Tanh (smooth for derivatives).
Numerical Methods
1. Spectral Methods
1.1 Fourier Transform
Discrete Fourier Transform:
$$\hat{f}k = \sum{n=0}^{N-1} f_n e^{-2\pi ikn/N}$$
Fast Fourier Transform (FFT): Complexity: $O(N\log N)$ vs. $O(N^2)$ for naive DFT.
Convolution Theorem:
$$\mathcal{F}{f * g} = \mathcal{F}{f} \cdot \mathcal{F}{g}$$
1.2 Spectral Derivatives
Spatial derivatives computed in Fourier space:
$$\frac{\partial^n f}{\partial x^n} = \mathcal{F}^{-1}{(2\pi i k)^n \hat{f}(k)}$$
Spectral accuracy: exponential convergence for smooth functions.
2. Time Integration
2.1 Runge-Kutta Methods
RK4 (4th order):
$$k_1 = f(t_n, y_n)$$
$$k_2 = f(t_n + \Delta t/2, y_n + \Delta t k_1/2)$$
$$k_3 = f(t_n + \Delta t/2, y_n + \Delta t k_2/2)$$
$$k_4 = f(t_n + \Delta t, y_n + \Delta t k_3)$$
$$y_{n+1} = y_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$
2.2 Implicit Methods
For stiff equations (radiation, diffusion):
Backward Euler:
$$y_{n+1} = y_n + \Delta t f(t_{n+1}, y_{n+1})$$
Requires solving nonlinear system via Newton iteration.
3. Stabilization Techniques
3.1 Upwind Schemes
For advection-dominated flows, use upwind differencing based on flow direction:
- If $u > 0$: use backward difference
- If $u < 0$: use forward difference
3.2 Artificial Viscosity
Add small diffusive term to damp high-frequency oscillations:
$$\nu_{\text{art}} \nabla^2 f$$
Conservation Laws and Constraints
1. Mass Conservation
Integral form:
$$\frac{d}{dt}\int_V \rho , dV + \int_S \rho \mathbf{u} \cdot \mathbf{n} , dS = 0$$
Discrete Check:
$$\sum_{i,j,k} \rho_{i,j,k} = \text{const}$$
Tolerance: $|\Delta M / M| < 10^{-6}$
2. Momentum Conservation
$$\frac{d}{dt}\int_V \rho \mathbf{u} , dV + \int_S \rho \mathbf{u}(\mathbf{u} \cdot \mathbf{n}) , dS = \int_V \mathbf{F} , dV$$
3. Energy Conservation
$$\frac{d}{dt}\int_V e , dV + \int_S (e + p)\mathbf{u} \cdot \mathbf{n} , dS = \int_V Q , dV$$
where $e = \rho(c_v T + \frac{1}{2}|\mathbf{u}|^2)$ is total energy density.
4. Radiative Equilibrium
In equilibrium:
$$\kappa B = S$$
where $B$ is Planck function:
$$B_\lambda(T) = \frac{2hc^2}{\lambda^5} \frac{1}{e^{hc/(\lambda k_B T)} - 1}$$
5. Thermodynamic Consistency
First Law:
$$dE = TdS - PdV$$
Entropy Production:
$$\frac{dS}{dt} \geq 0$$
Enforced via entropic projection in Boltzmann solver.
Physical Constants
| Symbol | Value | Description |
|---|---|---|
| $g$ | 9.80665 m/s² | Gravitational acceleration |
| $\Omega$ | $7.2921 \times 10^{-5}$ rad/s | Earth rotation rate |
| $R_d$ | 287.058 J/(kg·K) | Dry air gas constant |
| $R_v$ | 461.495 J/(kg·K) | Water vapor gas constant |
| $c_p$ | 1005.0 J/(kg·K) | Specific heat (constant pressure) |
| $c_v$ | 718.0 J/(kg·K) | Specific heat (constant volume) |
| $L_v$ | $2.5 \times 10^6$ J/kg | Latent heat of vaporization |
| $k_B$ | $1.381 \times 10^{-23}$ J/K | Boltzmann constant |
| $h$ | $6.626 \times 10^{-34}$ J·s | Planck constant |
| $c$ | $2.998 \times 10^8$ m/s | Speed of light |
| $\kappa$ | 0.41 | von Kármán constant |
Computational Complexity
| Component | Complexity | Notes |
|---|---|---|
| FFT (3D) | $O(N^3 \log N)$ | $N$ points per dimension |
| FNO Layer | $O(N^3 \log N)$ | Dominated by FFT |
| Boltzmann Solver | $O(N_v^3 N_x^3)$ | $N_v$ velocity points, $N_x$ spatial points |
| Ray Marching | $O(N_r N_s)$ | $N_r$ rays, $N_s$ steps |
| Monte Carlo | $O(N_p N_{\text{scatter}})$ | $N_p$ photons, avg scattering events |
| Discrete Ordinates | $O(N_\Omega N_z)$ | $N_\Omega$ angles, $N_z$ vertical levels |
References
MDNO References
- Boltzmann equation: Chapman & Cowling, "The Mathematical Theory of Non-uniform Gases" (1970)
- Fourier Neural Operators: Li et al., "Fourier Neural Operator for Parametric PDEs" (2021)
- LES turbulence: Pope, "Turbulent Flows" (2000)
- Hamiltonian mechanics: Goldstein, "Classical Mechanics" (1980)
RTNO References
- Radiative transfer: Chandrasekhar, "Radiative Transfer" (1960)
- Mie scattering: Bohren & Huffman, "Absorption and Scattering of Light by Small Particles" (1983)
- Polarization: van de Hulst, "Light Scattering by Small Particles" (1957)
- Correlated-k: Goody et al., "The Correlated-k Method for Radiation Calculations" (1989)
STRATUS References
- Volume rendering: Max, "Optical Models for Direct Volume Rendering" (1995)
- Monte Carlo: Marchuk et al., "The Monte Carlo Methods in Atmospheric Optics" (1980)
- PINN: Raissi et al., "Physics-Informed Neural Networks" (2019)
- Discrete ordinates: Stamnes et al., "Numerically stable algorithm for discrete-ordinate-method" (1988)
Notation Summary
| Symbol | Description |
|---|---|
| $f(\mathbf{x}, \mathbf{v}, t)$ | Particle distribution function |
| $\mathbf{u}, \mathbf{v}, \mathbf{w}$ | Velocity vector components |
| $T, p, \rho$ | Temperature, pressure, density |
| $I_\lambda$ | Spectral radiance |
| $\kappa, \sigma$ | Extinction/scattering coefficients |
| $\omega$ | Single scattering albedo |
| $\tau$ | Optical depth |
| $P(\Theta)$ | Phase function |
| $\mathbf{I}$ | Stokes vector (I, Q, U, V) |
| $\mathbf{M}$ | Mueller matrix |
| $Y_l^m$ | Spherical harmonics |
| $\mathcal{F}, \mathcal{F}^{-1}$ | Fourier transform, inverse |
Document Version: 1.0
Last Updated: 04.11.2025
Models: MDNO v5.3, RTNO v4.3, STRATUS