API Reference

HomogenousMaterialProperties{T}

Thermal properties with uniform rock properties (no stratification).

Simpler alternative to StratifiedMaterialProperties when rock is homogeneous.

StratifiedMaterialProperties{N,T}

Thermal properties for stratified rock layers and borehole materials.

Supports N horizontal rock layers with varying properties by depth. Includes material properties for water, steel, insulation, and cement backfill.

Fields

• k_rock_layers: thermal conductivities for each rock layer [W/(m·K)]
• rho_c_rock_layers: volumetric heat capacities for each layer [J/(m³·K)]
• layer_depths: depth boundaries for each layer [m]
• k_water, rho_c_water: water thermal properties
• k_steel, rho_c_steel: steel pipe thermal properties
• k_insulating, rho_c_insulating: insulation thermal properties
• k_backfill, rho_c_backfill: cement backfill thermal properties

Borehole{T}(xc, yc, h, r_inner, t_inner, r_outer, t_outer, r_backfill, ṁ, insulation_depth)

U-shaped coaxial borehole heat exchanger geometry.

Fluid flows down the inner pipe, turns around at depth h, and returns up the outer annulus. Precomputes squared radii and fluid velocities for performance.

Fields

• xc, yc: horizontal center coordinates [m]
• h: borehole depth [m]
• r_inner: inner pipe radius [m]
• t_inner: inner pipe wall thickness [m]
• r_outer: outer pipe inner radius [m]
• t_outer: outer pipe wall thickness [m]
• r_backfill: cement backfill outer radius [m]
• ṁ: mass flow rate [kg/s]
• insulation_depth: depth to which inner pipe is insulated [m]

ConstantInlet{T}(T_const)

Inlet model with constant temperature.

HeatExchangerInlet{T}(ΔT)

Inlet model applying temperature change ΔT [K] to outlet temperature.

Models a heat exchanger where inlet is T_outlet - ΔT. Typically computed as ΔT = Q / (ṁ * c_water) where Q [W] is heat extraction rate, ṁ [kg/s] is mass flow rate, and c_water [J/(kg·K)] is specific heat.

CustomInlet{F}(func)

Inlet model with custom function func(bh_idx, T_outlet_values, t).

Allows arbitrary time-varying or borehole-specific inlet temperatures.

create_adaptive_grid_1d(; xmin, xmax, dx_fine, growth_factor, dx_max, boreholes, backend, direction, Float_used=Float64)

Generate a non-uniform 1D grid with refined spacing around one or multiple boreholes.

Creates a 1D grid along the specified direction with fine resolution near boreholes and geometrically increasing spacing away from them, up to a maximum cell size. Works for both single borehole and multi-borehole configurations.

Arguments

• xmin, xmax: Domain boundaries [m]
• dx_fine: Fine grid spacing near boreholes [m]
• growth_factor: Geometric growth rate for coarsening (e.g., 1.3)
• dx_max: Maximum grid spacing [m]
• boreholes: Single borehole or collection of borehole objects with positions and radii
• backend: Computation backend (CPU/GPU) for array adaptation
• direction: Grid direction, either :x or :y
• Float_used=Float64: Floating point type for grid values

Grid generation strategy

For single or multiple boreholes:

1. Creates fine uniform grids around each borehole (spacing ≈ dx_fine)
2. Extends from domain boundaries toward boreholes with geometric growth
3. For multiple boreholes: fills gaps between boreholes symmetrically with progressive coarsening
4. Ensures all grid points are unique and sorted

Returns

Backend-adapted 1D grid array of the specified floating point type.

create_uniform_gridz_with_borehole_depths(; zmin, zmax, dz, boreholes, backend)

Create vertical grid including all borehole depths.

Ensures each borehole depth h is a grid point for accurate advection interpolation.

compute_domain(boreholes; buffer_x=100, buffer_y=100, buffer_z=200)

Compute simulation domain bounds from borehole positions.

Returns (xmin, xmax, ymin, ymax, zmin, zmax) encompassing all boreholes with specified buffers [m].

initial_condition_thermal_gradient(backend, gridx, gridy, gridz; T_surface, gradient, Float_used=Float64)
initial_condition_thermal_gradient(backend, Float_used, gridx, gridy, gridz; T_surface, gradient)

Create initial temperature field with linear thermal gradient.

Returns 3D array with T(z) = T_surface + gradient * z where T_surface is surface temperature [°C] and gradient is thermal gradient [K/m or °C/m].

Float_used defaults to Float64. The legacy 5-argument form (with Float_used as the second positional argument) is kept for backward compatibility.

create_cache(; backend, gridx, gridy, gridz, materials, boreholes, inlet_model, precompute_materials=true)

Create simulation cache.

By default, precompute_materials=:auto chooses the material evaluation mode based on the number of boreholes:

• one borehole: use :on_the_fly
• multiple boreholes: use :precomputed

This reflects the observed performance behavior: on-the-fly material lookup is slightly faster for a single well, while precomputed material arrays are significantly faster for well arrays.

You can override the automatic choice manually:

• precompute_materials=true or :precomputed
• precompute_materials=false or :on_the_fly

Returns named tuple containing grids, materials, index lists, outlet temperature arrays, eigenvalue estimates, material accessor, and precomputed Val types for kernel dispatch.

rhs_diffusion_z!(dϕ, ϕ, cache, t)

Compute right-hand side for vertical (z-direction) heat diffusion.

This is the RHS function passed to OrdinaryDiffEq.jl and is integrated using an explicit stabilized ROCK method. Vertical diffusion is treated separately from horizontal (x,y) diffusion because the ADI scheme is only unconditionally stable in 2D. Since Δz is much larger than the smallest Δx and Δy (due to the fine grid resolution needed near the borehole), the explicit ROCK method is sufficient for the z-direction without imposing prohibitive time step restrictions.

Horizontal diffusion via the ADI method and advection are handled separately using the callback functionality (ADI_and_ADV_callback!)—a workaround to implement operator splitting within the OrdinaryDiffEq.jl framework.

eigen_estimator(integrator)

Eigenvalue estimator function for adaptive time stepping in ROCK methods.

get_simulation_callback(; saveat, print_every_n=1000,
                          checkpoint_dir="", checkpoint_id="latest",
                          checkpoint_every_n=0, Float_used_to_save=Float32)

Create the required callback set for the simulation.

The callback combines up to four components:

1. ADI + Advection: Horizontal diffusion and fluid advection (runs every timestep)
2. Checkpointing (optional): Periodically saves (u, t) to disk for fault-tolerant restarts
3. Progress printing: Prints simulation progress every print_every_n steps
4. Solution saving: Saves temperature field at times specified by saveat. When checkpoint_dir is provided, snapshots are also written to disk so that data persists across crash/restart cycles.

The checkpoint callback is placed directly after the ADI + advection callback in the CallbackSet ordering. This ensures that a checkpoint always represents a fully consistent state (both ROCK2 vertical diffusion and ADI horizontal diffusion + advection have been applied), so restarting from a checkpoint produces identical results to an uninterrupted run.

Arguments

• saveat: Times at which to save the solution (e.g., range(0, 3600, 10) or [0.0, 3600.0])
• print_every_n=1000: Print progress every N accepted timesteps
• checkpoint_dir="": Directory for checkpoint and snapshot files. Empty string disables both checkpointing and persistent snapshots (solutions are only kept in memory).
• checkpoint_id="latest": Unique identifier for checkpoint/snapshot files (useful when multiple simulations share the same directory). Files are named checkpoint_{id}.jld2 (restart) and snapshot_{id}_NNNN.jld2 (data).
• checkpoint_every_n=0: Save a restart checkpoint every N accepted timesteps. Set to 0 to disable restart checkpointing (snapshots at saveat times are still written if checkpoint_dir is provided).
• Float_used_to_save=Float32: Floating point type for saved solution snapshots. Restart checkpoints always use full Float64 precision.

Returns

• callback: Combined CallbackSet to pass to solve(..., callback=callback)
• saved_values: SavedValues object. After solve, contains only snapshots from the current run. Call reload_snapshots! to populate with the full history from disk.

Example

callback, saved_values = get_simulation_callback(
    saveat=saveat,
    print_every_n=100_000,
    checkpoint_dir="output/",
    checkpoint_id="my_simulation",
    checkpoint_every_n=500_000
)
solve(prob, ROCK2(eigen_est=eigen_estimator), callback=callback, dt=80.0, adaptive=false)
 
# Assemble full history from all snapshot files (across all crash/restart cycles)
reload_snapshots!(saved_values, "output/", "my_simulation")

prepare_restart(T0, tspan, saveat; checkpoint_dir, checkpoint_id="latest", backend=CPU())

Check for an existing checkpoint and prepare the simulation for a fresh start or a restart.

If a checkpoint file exists at checkpoint_dir/checkpoint_{checkpoint_id}.jld2, the saved temperature field and time are loaded. The initial condition is replaced, tspan is adjusted to start from the checkpoint time, and saveat is filtered to only include future save times.

If no checkpoint file exists, all arguments are returned unchanged (identity operation).

This function pairs with get_simulation_callback: use the same checkpoint_dir and checkpoint_id for both so that checkpoints, snapshots, and restart logic all align.

Arguments

• T0: Fresh initial condition (e.g., from initial_condition_thermal_gradient)
• tspan: Full time span (t_start, t_end)
• saveat: Save times (any iterable of times)
• checkpoint_dir: Directory where checkpoint files are stored
• checkpoint_id="latest": Unique identifier matching the one used in get_simulation_callback
• backend=CPU(): Computation backend to adapt the loaded array to

Returns

(T0, tspan, saveat) – either unchanged (no checkpoint) or updated for restart.

Example

T0_fresh = initial_condition_thermal_gradient(backend, Float64, gridx, gridy, gridz;
    T_surface=10.0, gradient=0.035)
tspan_full = (0.0, 3600.0 * 24 * 365 * 20)
saveat_full = range(tspan_full..., 21)
 
T0, tspan, saveat = prepare_restart(
    T0_fresh, tspan_full, saveat_full;
    checkpoint_dir="output/",
    checkpoint_id="my_sim",
    backend=CUDABackend()
)
 
prob = ODEProblem(rhs_diffusion_z!, T0, tspan, cache)
callback, saved_values = get_simulation_callback(saveat=saveat,
    checkpoint_dir="output/", checkpoint_id="my_sim", checkpoint_every_n=500_000)
solve(prob, ROCK2(...), callback=callback, ...)
 
# Assemble full history from disk
reload_snapshots!(saved_values, "output/", "my_sim")

reload_snapshots!(saved_values, checkpoint_dir, checkpoint_id; Float_used_to_save=Float32)

Load all snapshot files from disk into saved_values, sorted by time.

This replaces the contents of saved_values.t and saved_values.saveval with the full history from all snapshot files matching the given checkpoint_id. Call this after solve to assemble the complete history across all crash/restart cycles.

The ODE solver clears saved_values at the start of each solve() call, so without calling this function, saved_values only contains snapshots from the most recent run.

Arguments

• saved_values: The SavedValues object returned by get_simulation_callback
• checkpoint_dir: Directory containing snapshot files
• checkpoint_id: Unique identifier matching the one used in get_simulation_callback
• Float_used_to_save=Float32: Floating point type matching the one used in get_simulation_callback

Example

callback, saved_values = get_simulation_callback(
    saveat=saveat, checkpoint_dir="output/", checkpoint_id="my_sim", checkpoint_every_n=500_000)
solve(prob, ROCK2(...), callback=callback, ...)
 
# After solve, saved_values only has this run's data.
# Reload to get the full history:
reload_snapshots!(saved_values, "output/", "my_sim")
 
# Now saved_values.t and saved_values.saveval contain ALL snapshots, sorted by time.
T_final = saved_values.saveval[end]

pkg_dir()

Return package root directory.

examples_dir()

Return path to package examples directory.

data_dir()

Return path to package data directory.

plot_grid(gridx, gridy; boreholes=nothing, xlims=nothing, ylims=nothing, ...)

Visualize the computational grid in the x-y plane.

Plots grid lines and optionally overlays borehole geometry (inner pipe, outer pipe, backfill boundaries). Useful for verifying grid resolution near boreholes.

Arguments

• gridx: Grid coordinates in x-direction
• gridy: Grid coordinates in y-direction
• boreholes=nothing: Tuple of Borehole objects to overlay
• xlims=nothing: x-axis limits (auto-determined if not provided)
• ylims=nothing: y-axis limits (auto-determined if not provided)
• size=(800, 800): Figure size in pixels
• annotate=true: Add text labels to borehole components
• legend=:bottomleft: Legend position

Returns

A Plots.jl plot object.

data_li(i) -> x_data, T_data

Load Li et al. validation data (Figure 5).

i=1,2,3,4 for 500m, 1000m, 1500m, 2000m depths respectively. From: Heat extraction model and characteristics of coaxial deep borehole heat exchanger

data_hu(years::Int) -> depth, temperature

Load Hu et al. validation data for specified simulation years (Figure 7).

From: Numerical modeling of a coaxial borehole heat exchanger (doi:10.1016/j.renene.2019.09.141) Returns depth [m] and temperature [°C].

data_brown_single_well_b() -> (z_beier, T_beier), (z_brown, T_brown)

Load Brown et al. single well validation data (Figure 4b).

From: Investigating scalability of deep borehole heat exchangers (doi:10.1016/j.renene.2021.01.036) Returns depth [m] and temperature [°C] profiles for Beier and Brown models.

data_brown_single_well_d(i) -> x_data, T_data

Load Brown et al. single well validation data (Figure 4d).

i=1,2,3 for 300m, 600m, 920m depths respectively. From: doi:10.1016/j.renene.2021.01.036

data_brown_array(spacing, i) -> x_data, T_data

Load Brown et al. array validation data (Figure 6).

spacing=20,30,40,50 for array spacing in meters. i=1,2,3 for 300m, 600m, 920m depths respectively. From: doi:10.1016/j.renene.2021.01.036