Spherical Geometry Problems

Overview

Dedalus provides three bases for spherical coordinate problems:

SphereBasis: For problems on the surface of a sphere (θ, φ)
ShellBasis: For problems in a spherical shell (φ, θ, r ∈ [R_i, R_o])
BallBasis: For problems inside a full sphere/ball (φ, θ, r ∈ [0, R])

These geometries are essential for geophysical, astrophysical, and planetary science applications.

Shallow Water on Sphere

File: examples/ivp_sphere_shallow_water/shallow_water.py

Description

Simulates the viscous shallow water equations on a rotating sphere, implementing the classic test case of a barotropically unstable mid-latitude jet from Galewsky et al. (2004). This example demonstrates:

Solving an LBVP for balanced initial conditions
Time-evolving an IVP on the sphere
Including rotation (Coriolis force)
Using dimensional units

Physical Setup

Geometry: Full sphere (Earth)
Resolution: (N_φ, N_θ) = (256, 128)
Initial condition: Zonal jet at mid-latitudes
Parameters:
- Radius: 6.37×10⁶ m
- Rotation: Ω = 7.29×10⁻⁵ s⁻¹
- Gravity: g = 9.81 m/s²
- Mean height: H = 10⁴ m

Setting Up the Geometry

# Bases
coords = d3.S2Coordinates('phi', 'theta')
dist = d3.Distributor(coords, dtype=dtype)
basis = d3.SphereBasis(coords, (Nphi, Ntheta), radius=R, dealias=dealias, dtype=dtype)

# Fields
u = dist.VectorField(coords, name='u', bases=basis)
h = dist.Field(name='h', bases=basis)

Initial Conditions: Balanced Jet

First, set up a zonal jet:

# Grid setup (phi = longitude, theta = colatitude)
phi, theta = dist.local_grids(basis)
lat = np.pi / 2 - theta  # Convert to latitude

# Jet profile
umax = 80 * meter / second
lat0 = np.pi / 7
lat1 = np.pi / 2 - lat0
en = np.exp(-4 / (lat1 - lat0)**2)
jet = (lat0 <= lat) * (lat <= lat1)
u_jet = umax / en * np.exp(1 / (lat[jet] - lat0) / (lat[jet] - lat1))
u['g'][0][jet] = u_jet

Then solve an LBVP for the balanced height field:

# Coriolis operator
zcross = lambda A: d3.MulCosine(d3.skew(A))

# Solve for balanced initial height
c = dist.Field(name='c')
problem = d3.LBVP([h, c], namespace=locals())
problem.add_equation("g*lap(h) + c = - div(u@grad(u) + 2*Omega*zcross(u))")
problem.add_equation("ave(h) = 0")
solver = problem.build_solver()
solver.solve()

Finally, add a localized perturbation:

hpert = 120 * meter
h['g'] += hpert * np.cos(lat) * np.exp(-(phi/alpha)**2) * np.exp(-((lat2-lat)/beta)**2)

Governing Equations

# Shallow water IVP with rotation and hyperdiffusion
problem = d3.IVP([u, h], namespace=locals())
problem.add_equation("dt(u) + nu*lap(lap(u)) + g*grad(h) + 2*Omega*zcross(u) = - u@grad(u)")
problem.add_equation("dt(h) + nu*lap(lap(h)) + H*div(u) = - div(h*u)")

Running the Example

mpiexec -n 4 python3 shallow_water.py
mpiexec -n 4 python3 plot_sphere.py snapshots/*.h5

Expected Results

Runtime: ~5 cpu-minutes
Physics: Mid-latitude jet becomes unstable, forming Rossby wave patterns
Output: Height and vorticity fields on the sphere

Shell Convection

File: examples/ivp_shell_convection/shell_convection.py

Description

Simulates Boussinesq convection in a spherical shell—the geometry relevant to planetary mantles and stellar convection zones. This example demonstrates the full 3D spherical shell geometry with radial boundary conditions.

Physical Setup

Geometry: Spherical shell with R_i = 14, R_o = 15 (thin shell)
Resolution: (N_φ, N_θ, N_r) = (192, 96, 6)
Boundary conditions: Fixed temperatures, no-slip walls
Parameters: Ra = 3500, Pr = 1
Non-dimensionalization: Shell thickness and freefall time

Setting Up the Shell

coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype, mesh=mesh)
shell = d3.ShellBasis(coords, shape=(Nphi, Ntheta, Nr), radii=(Ri, Ro), 
                      dealias=dealias, dtype=dtype)
sphere = shell.outer_surface  # For boundary conditions

Radial Direction

Define radial unit vector and position vector:

phi, theta, r = dist.local_grids(shell)

# Radial unit vector
er = dist.VectorField(coords, bases=shell.radial_basis)
er['g'][2] = 1

# Position vector
rvec = dist.VectorField(coords, bases=shell.radial_basis)
rvec['g'][2] = r

First-Order Formulation

lift_basis = shell.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
grad_u = d3.grad(u) + rvec*lift(tau_u1)
grad_b = d3.grad(b) + rvec*lift(tau_b1)

problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*er + lift(tau_u2) = - u@grad(u)")

Boundary Conditions

problem.add_equation("b(r=Ri) = 1")  # Hot inner boundary
problem.add_equation("u(r=Ri) = 0")  # No-slip
problem.add_equation("b(r=Ro) = 0")  # Cold outer boundary
problem.add_equation("u(r=Ro) = 0")  # No-slip
problem.add_equation("integ(p) = 0")  # Pressure gauge

Initial Conditions

# Random perturbation
b.fill_random('g', seed=42, distribution='normal', scale=1e-3)
b['g'] *= (r - Ri) * (Ro - r)  # Damp at boundaries

# Add conductive background
b['g'] += (Ri - Ri*Ro/r) / (Ri - Ro)

Running the Example

mpiexec -n 4 python3 shell_convection.py
mpiexec -n 4 python3 plot_shell.py snapshots/*.h5

Expected Results

Runtime: ~20 cpu-minutes
Physics: Convective plumes form and rise/fall in the shell
Output: Temperature slices at various radii and longitudes

Ball Internally Heated Convection

File: examples/ivp_ball_internally_heated_convection/internally_heated_convection.py

Description

Simulates internally-heated Boussinesq convection inside a full sphere (ball). The gravity is proportional to radius (constant density). Features stress-free boundary conditions and demonstrates restart capability.

Physical Setup

Geometry: Full ball with radius = 1
Resolution: (N_φ, N_θ, N_r) = (128, 64, 96)
Boundary conditions: Stress-free, constant flux
Parameters: Ra = 10⁶, Pr = 1
Heating: Volumetric internal heating (T_source = 6)
Equilibrium: Conductive state T(r) = 1 - r²

Setting Up the Ball

coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype, mesh=mesh)
ball = d3.BallBasis(coords, shape=(Nphi, Ntheta, Nr), radius=1, 
                    dealias=dealias, dtype=dtype)
sphere = ball.surface

Stress-Free Boundary Conditions

# Stress tensor
strain_rate = d3.grad(u) + d3.trans(d3.grad(u))
shear_stress = d3.angular(d3.radial(strain_rate(r=1), index=1))

problem.add_equation("shear_stress = 0")  # Stress-free
problem.add_equation("radial(u(r=1)) = 0")  # No penetration
problem.add_equation("radial(grad(T)(r=1)) = -2")  # Constant flux

Internal Heating

# Temperature equation with source term
problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T) = - u@grad(T) + kappa*T_source")

Restart Capability

import sys
restart = (len(sys.argv) > 1 and sys.argv[1] == '--restart')

if not restart:
    # Initial conditions
    T.fill_random('g', seed=42, distribution='normal', scale=0.01)
    T.low_pass_filter(scales=0.5)
    T['g'] += 1 - r**2
    file_handler_mode = 'overwrite'
else:
    # Restart from checkpoint
    write, initial_timestep = solver.load_state('checkpoints/checkpoints_s20.h5')
    file_handler_mode = 'append'

Running the Example

# Initial run
mpiexec -n 4 python3 internally_heated_convection.py

# Restart for additional time
mpiexec -n 4 python3 internally_heated_convection.py --restart

# Plot results
mpiexec -n 4 python3 plot_ball.py slices/*.h5

Expected Results

Runtime: ~30 cpu-minutes per segment
Physics: Convective plumes rise from interior to boundary
Output: Temperature slices at various radii and longitudes

Lane-Emden Equation

File: examples/nlbvp_ball_lane_emden/lane_emden.py

Description

Solves the Lane-Emden equation—a nonlinear boundary value problem describing polytropic stellar structure. This example demonstrates:

Nonlinear boundary value problems (NLBVP)
Newton iteration
Spherically symmetric problems (1D in ball)
Convergence to reference solutions

Mathematical Problem

The Lane-Emden equation in rescaled form:

∇²f + f^n = 0
f(r=1) = 0

where n is the polytropic index (n = 3 in the example).

Setting Up the NLBVP

# Bases (spherically symmetric: 1 mode in phi and theta)
coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype)
ball = d3.BallBasis(coords, (1, 1, Nr), radius=1, dtype=dtype, dealias=dealias)

# Fields
f = dist.Field(name='f', bases=ball)
tau = dist.Field(name='tau', bases=ball.surface)

# Problem
lift = lambda A: d3.Lift(A, ball, -1)
problem = d3.NLBVP([f, tau], namespace=locals())
problem.add_equation("lap(f) + lift(tau) = - f**n")
problem.add_equation("f(r=1) = 0")

Initial Guess

phi, theta, r = dist.local_grids(ball)
R0 = 5  # Initial guess for eigenvalue
f['g'] = R0**(2/(n-1)) * (1 - r**2)**2

Newton Iteration

solver = problem.build_solver(ncc_cutoff=ncc_cutoff)
pert_norm = np.inf
tolerance = 1e-10

while pert_norm > tolerance:
    solver.newton_iteration()
    pert_norm = sum(pert.allreduce_data_norm('c', 2) for pert in solver.perturbations)
    logger.info(f'Perturbation norm: {pert_norm:.3e}')
    
    # Extract eigenvalue R from f(r=0)
    f0 = f(r=0).evaluate().allgather_data('g')[0,0,0]
    Ri = f0**((n-1)/2)
    logger.info(f'R iterate: {Ri}')

Running the Example

python3 lane_emden.py

Expected Results

Runtime: A few seconds
Convergence: Typically ~10-12 Newton iterations
Output:
- Convergence information printed to console
- Plot showing Newton iteration progress lane_emden.png
- For n = 3: R ≈ 6.8968 (matches reference to high precision)

Reference Solutions

The example includes reference values from Boyd (2011) for validation:

R_ref = {0.0: np.sqrt(6),
         1.0: np.pi,
         2.0: 4.3528745959461246769735700,
         3.0: 6.896848619376960375454528,
         4.0: 14.971546348838095097611066}

Spherical Coordinate Tips

Surface Operators

# Extract surface values
field_at_surface = field(r=R).evaluate()

# Surface integrals
surface_average = d3.Average(field(r=R))

Radial Operations

# Radial derivative
df_dr = d3.radial(d3.grad(f))

# Radial component of vector
u_r = d3.radial(u)

# Angular components
u_angular = d3.angular(u)

Spherically Symmetric Problems

For 1D (spherically symmetric) problems, use shape (1, 1, Nr):

ball = d3.BallBasis(coords, (1, 1, Nr), radius=1, dtype=dtype)

Operators Specific to Spherical Surfaces

# Coriolis-like operator on sphere
zcross = lambda A: d3.MulCosine(d3.skew(A))

# Averaging on sphere
c = d3.Average(h)  # For constraint equations

Get Started

Core Concepts

Advanced Topics

Examples

​Overview

​Shallow Water on Sphere

​Description

​Physical Setup

​Setting Up the Geometry

​Initial Conditions: Balanced Jet

​Governing Equations

​Running the Example

​Expected Results

​Shell Convection

​Description

​Physical Setup

​Setting Up the Shell

​Radial Direction

​First-Order Formulation

​Boundary Conditions

​Initial Conditions

​Running the Example

​Expected Results

​Ball Internally Heated Convection

​Description

​Physical Setup

​Setting Up the Ball

​Stress-Free Boundary Conditions

​Internal Heating

​Restart Capability

​Running the Example

​Expected Results

​Lane-Emden Equation

​Description

​Mathematical Problem

​Setting Up the NLBVP

​Initial Guess

​Newton Iteration

​Running the Example

​Expected Results

​Reference Solutions

​Spherical Coordinate Tips

​Surface Operators

​Radial Operations

​Spherically Symmetric Problems

​Operators Specific to Spherical Surfaces

Build docs developers (and LLMs) love

Overview

Shallow Water on Sphere

Description

Physical Setup

Setting Up the Geometry

Initial Conditions: Balanced Jet

Governing Equations

Running the Example

Expected Results

Shell Convection

Description

Physical Setup

Setting Up the Shell

Radial Direction

First-Order Formulation

Boundary Conditions

Initial Conditions

Running the Example

Expected Results

Ball Internally Heated Convection

Description

Physical Setup

Setting Up the Ball

Stress-Free Boundary Conditions

Internal Heating

Restart Capability

Running the Example

Expected Results

Lane-Emden Equation

Description

Mathematical Problem

Setting Up the NLBVP

Initial Guess

Newton Iteration

Running the Example

Expected Results

Reference Solutions

Spherical Coordinate Tips

Surface Operators

Radial Operations

Spherically Symmetric Problems

Operators Specific to Spherical Surfaces