Full tutorial

This tutorial covers every major BLDFM capability, organized by what you want to accomplish. It assumes you have completed the Quickstart tutorial and have BLDFM installed.

Prerequisites

$ pip install -e ".[dev,plotting]"

All code below assumes you are in the repository root and have an active Python session. Small grid sizes are used throughout for fast execution.

import bldfm
bldfm.initialize()

Configuration

BLDFM configurations define the computational domain, tower positions, meteorological forcing, and solver settings. You can load configs from YAML files or build them programmatically in Python.

Loading a YAML config

BLDFM ships with several example configs in examples/configs/. load_config() returns a BLDFMConfig dataclass:

config = bldfm.load_config("examples/configs/footprint.yaml")

print(f"Domain: {config.domain.nx}x{config.domain.ny}, nz={config.domain.nz}")
print(f"Tower: {config.towers[0].name}, z_m={config.towers[0].z_m}")
print(f"Solver: closure={config.solver.closure}, footprint={config.solver.footprint}")

Building a config in Python

For programmatic workflows – parameter sweeps, scripted experiments – use parse_config_dict() with a plain Python dictionary that mirrors the YAML structure. When ref_lat/ref_lon are set, each tower’s local (x, y) coordinates are computed automatically from its (lat, lon) via an equirectangular projection:

from bldfm import parse_config_dict

config = parse_config_dict({
    "domain": {
        "nx": 128, "ny": 64, "xmax": 500.0, "ymax": 250.0, "nz": 16,
        "ref_lat": 50.95, "ref_lon": 11.586,
    },
    "towers": [{"name": "tower_A", "lat": 50.9505, "lon": 11.5865, "z_m": 10.0}],
    "met": {"ustar": 0.4, "mol": -100.0, "wind_speed": 5.0, "wind_dir": 270.0},
    "solver": {"closure": "MOST", "footprint": True},
})

print(f"Tower local coords: x={config.towers[0].x:.1f} m, y={config.towers[0].y:.1f} m")

Config schema reference

A BLDFMConfig contains six sub-configs:

DomainConfig – computational grid geometry.

  • nx, ny: grid points in x (cross-wind) and y (along-wind)

  • nz: number of vertical levels

  • xmax, ymax: domain extent in metres

  • modes: Fourier modes [kx, ky]; higher values improve accuracy

  • halo: zero-padding width (metres) to reduce spectral leakage

  • ref_lat, ref_lon: reference origin for lat/lon to local coordinate conversion

TowerConfig – measurement tower location.

  • name, lat, lon, z_m: tower identity, coordinates, and height

  • x, y: local coordinates (auto-computed from ref_lat/ref_lon)

MetConfig – meteorological forcing. Use scalars for a single timestep or lists of equal length for a timeseries.

  • ustar: friction velocity (m/s). Provide either ustar or z0, not both. When z0 (roughness length) is given, the model derives ustar internally via the PBL closure.

  • mol: Monin-Obukhov length (m). Negative = unstable (daytime convection), positive = stable (nighttime), very large (~1e9) = neutral.

  • wind_speed, wind_dir: horizontal wind speed (m/s) and direction (degrees; 0 = N, 90 = E, 180 = S, 270 = W).

  • timestamps: optional list of timestamp labels for each timestep.

SolverConfig – numerical solver settings.

  • closure: PBL closure scheme – MOST, MOSTM, CONSTANT, or OAAHOC.

  • footprint: true computes a flux footprint (Green’s function at measurement height); false computes a concentration field from a surface source.

  • precision: single or double floating-point arithmetic.

  • analytic: if true, uses the Kormann-Meixner analytical footprint model for comparison.

OutputConfig – output format and directory.

ParallelConfig – parallelism settings.

  • num_threads: BLAS/numba threads per worker.

  • max_workers: number of parallel worker processes.

  • use_cache: enable disk-based solver result caching.

Running simulations

BLDFM provides four functions at progressive levels of complexity, from a single solve to fully parallel multi-tower timeseries.

Single solve

run_bldfm_single() is the fundamental building block. It takes a config and a tower, runs one timestep, and returns a result dictionary:

from bldfm import run_bldfm_single

result = run_bldfm_single(config, config.towers[0])

print(f"Keys: {sorted(result.keys())}")
print(f"Footprint shape: {result['flx'].shape}")  # (ny, nx)
print(f"Met used: wind_speed={result['params']['wind_speed']}")

The result dict contains:

  • grid: tuple of (X, Y, Z) coordinate arrays

  • flx: 2D flux footprint field, shape (ny, nx)

  • conc: 2D concentration field

  • tower_name, tower_xy: tower identity and local coordinates

  • timestamp: timestep label or index

  • params: dict of meteorological parameters used for this solve

Timeseries

For time-varying meteorology, provide lists in the met config. run_bldfm_timeseries() loops over all timesteps for a single tower and returns a list of result dicts. Use generate_synthetic_timeseries() to create reproducible test data:

from bldfm import run_bldfm_timeseries, generate_synthetic_timeseries

met = generate_synthetic_timeseries(n_timesteps=3, seed=42)
config_ts = parse_config_dict({
    "domain": {
        "nx": 128, "ny": 64, "xmax": 500.0, "ymax": 250.0, "nz": 16,
        "ref_lat": 50.95, "ref_lon": 11.586,
    },
    "towers": [{"name": "tower_A", "lat": 50.9505, "lon": 11.5865, "z_m": 10.0}],
    "met": met,
    "solver": {"closure": "MOST", "footprint": True},
})

results = run_bldfm_timeseries(config_ts, config_ts.towers[0])
print(f"Timesteps computed: {len(results)}")
for r in results:
    print(f"  t={r['timestamp']}: max flux = {r['flx'].max():.6f}")

Multi-tower

run_bldfm_multitower() runs all towers across all timesteps, returning a dict of {tower_name: [result, ...]}. Use generate_towers_grid() to create synthetic tower layouts ("grid", "transect", or "random"):

from bldfm import run_bldfm_multitower, generate_towers_grid

towers = generate_towers_grid(n_towers=2, z_m=10.0, layout="transect", seed=42)
config_mt = parse_config_dict({
    "domain": {
        "nx": 128, "ny": 64, "xmax": 500.0, "ymax": 250.0, "nz": 16,
        "ref_lat": towers[0]["lat"], "ref_lon": towers[0]["lon"],
    },
    "towers": towers,
    "met": met,
    "solver": {"closure": "MOST", "footprint": True},
})

all_results = run_bldfm_multitower(config_mt)
for name, res_list in all_results.items():
    print(f"{name}: {len(res_list)} timesteps")

Parallel execution

For large runs, run_bldfm_parallel() distributes work across CPU cores using ProcessPoolExecutor. Three strategies are available:

  • "towers": each worker handles one tower’s full timeseries

  • "time": distribute timesteps across workers (per tower)

  • "both": flatten all tower x timestep pairs across workers

Each subprocess automatically sets NUMBA_NUM_THREADS=1 to avoid CPU oversubscription. The return format is identical to run_bldfm_multitower():

from bldfm import run_bldfm_parallel

par_results = run_bldfm_parallel(config_mt, max_workers=2, parallel_over="towers")

Working with results

Solver results are plain Python dictionaries containing NumPy arrays. This section covers common post-processing tasks.

Caching

Green’s function caching

In footprint mode, the solver output (Green’s function) depends only on vertical profiles and domain geometry – not on the surface flux. The GreensFunctionCache stores results as SHA-256-keyed .npz files so that identical configurations skip the solver entirely:

from bldfm import GreensFunctionCache, run_bldfm_single

cache = GreensFunctionCache(cache_dir=".bldfm_cache")

result1 = run_bldfm_single(config, config.towers[0], cache=cache)  # cache miss
result2 = run_bldfm_single(config, config.towers[0], cache=cache)  # cache hit (fast)

cache.clear()  # clean up

Caching can also be enabled globally via parallel.use_cache: true in the YAML config.

Visualization

Plotting functions are provided by abltk.plotting. Core plots require only matplotlib; optional dependencies (contextily, windrose, plotly) are imported lazily and produce helpful messages when missing.

import matplotlib
matplotlib.use("Agg")          # remove for interactive display
import matplotlib.pyplot as plt

Footprint field

plot_footprint_field() renders a 2D pcolormesh with optional percentile contour overlays. It works for both footprint and concentration fields:

from abltk.plotting import plot_footprint_field

ax = plot_footprint_field(
    result["flx"], result["grid"],
    contour_pcts=[0.5, 0.8],
    title="Footprint with 50% and 80% contours",
)
plt.savefig("plots/tutorial_contours.png", dpi=150, bbox_inches="tight")
plt.close()

Map overlay (requires contextily)

plot_footprint_on_map() overlays footprint contours and tower markers on web map tiles. Set land_cover=True to use ESA WorldCover 2021 instead of street maps (requires owslib):

try:
    from abltk.plotting import plot_footprint_on_map

    ax = plot_footprint_on_map(
        result["flx"], result["grid"], config,
        tower=config.towers[0],
        contour_pcts=[0.5, 0.8],
    )
    plt.savefig("plots/tutorial_map.png", dpi=150, bbox_inches="tight")
    plt.close()
except ImportError:
    print("Install contextily for map overlays: pip install contextily")

Wind rose (requires windrose)

plot_wind_rose() creates a polar wind rose from wind speed and direction arrays:

try:
    from abltk.plotting import plot_wind_rose

    ax = plot_wind_rose(
        met["wind_speed"], met["wind_dir"],
        title="Synthetic wind rose",
    )
    plt.savefig("plots/tutorial_windrose.png", dpi=150, bbox_inches="tight")
    plt.close()
except ImportError:
    print("Install windrose: pip install windrose")

Footprint timeseries

plot_footprint_timeseries() tracks how footprint extent (area at given percentiles) changes across timesteps:

from abltk.plotting import plot_footprint_timeseries

ax = plot_footprint_timeseries(
    results, results[0]["grid"],
    pcts=[0.5, 0.8],
    title="Footprint area evolution",
)
plt.savefig("plots/tutorial_timeseries.png", dpi=150, bbox_inches="tight")
plt.close()

Interactive plot (requires plotly)

plot_footprint_interactive() creates a Plotly figure that can be saved as HTML for exploration in a browser:

try:
    from abltk.plotting import plot_footprint_interactive

    fig = plot_footprint_interactive(result["flx"], result["grid"])
    fig.write_html("plots/tutorial_interactive.html")
except ImportError:
    print("Install plotly: pip install plotly")

Diagnostic plots

Functions for inspecting solver behaviour and vertical structure (plot_convergence(), plot_vertical_profiles(), plot_vertical_slice()) are available in abltk.plotting. See the abl-tk documentation for full signatures and usage.

Low-level API

For full control over every step of the simulation, you can call the solver pipeline directly. The high-level run_bldfm_single() wraps three steps:

  1. compute_wind_fields(speed, direction) – decompose wind into (u, v) components.

  2. vertical_profiles(nz, z_m, wind, ustar) – compute vertical profiles of wind and eddy diffusivity. Returns (z, profiles) where profiles is a 5-tuple (u, v, Kx, Ky, Kz).

  3. steady_state_transport_solver(srf_flx, z, profiles, domain, levels, ...) – solve the advection-diffusion equation.

from bldfm import compute_wind_fields, ideal_source, steady_state_transport_solver
from bldfm.pbl_model import vertical_profiles

# Step 1: wind components
u, v = compute_wind_fields(wind_speed=5.0, wind_dir=270.0)

# Step 2: vertical profiles
z, profiles = vertical_profiles(n=32, meas_height=10.0, wind=(u, v), ustar=0.4)
u_prof, v_prof, Kx, Ky, Kz = profiles

# Step 3: surface flux (or use np.zeros for footprint mode)
srf_flx = ideal_source((256, 128), (1000.0, 500.0))

# Step 4: solve
grid, conc, flx = steady_state_transport_solver(
    srf_flx, z, profiles, domain=(1000.0, 500.0), levels=32,
)

See the examples/low_level/ directory for complete scripts: minimal_example.py, footprint_example.py, plot_profiles.py, and more. For source area analysis, see runs/low_level/source_area_example.py.

Command-line interface

The bldfm CLI wraps the full workflow and calls initialize() automatically:

# Validate a config without running the solver
$ bldfm run examples/configs/footprint.yaml --dry-run

# Run all towers and timesteps
$ bldfm run examples/configs/multitower.yaml

# Run and save footprint plots to plots/
$ bldfm run examples/configs/multitower.yaml --plot

Note

--plot generates quick-inspection plots with default settings (concentration + flux/footprint fields, or vertical slices for 3D runs). For publication-quality figures with full control over layout and styling, run the example scripts in examples/ directly. See Example Workflows for the complete list.

Further resources

  • Quick Reference: concise code snippets for common workflows.

  • Example Scripts: examples/ (config-driven), examples/low_level/ (direct API), and runs/manuscript/ (paper reproduction).

  • API Reference: full function signatures and docstrings.

  • Manuscript figures: python runs/manuscript/generate_all.py regenerates all paper figures.