GitHub

Getting Started

This guide walks through a complete simulation of a single deep borehole heat exchanger.

Overview

A typical GeothermalWells.jl simulation involves:

  1. Defining material properties (rock, water, steel, etc.)
  2. Specifying borehole geometry
  3. Creating computational grids
  4. Setting initial conditions
  5. Choosing an inlet model
  6. Running the simulation

Complete Example

using GeothermalWells
using OrdinaryDiffEqStabilizedRK: ODEProblem, solve, ROCK2
using KernelAbstractions: CPU

# Choose backend: CPU() for testing, CUDABackend() for NVIDIA GPU
backend = CPU()

Material Properties

Define thermal properties for all materials in the simulation:

materials = HomogenousMaterialProperties(
    2.88,       # k_rock - rock thermal conductivity [W/(m·K)]
    2.17e6,     # rho_c_rock - rock volumetric heat capacity [J/(m³·K)]
    0.6,        # k_water [W/(m·K)]
    4.18e6,     # rho_c_water [J/(m³·K)]
    44.5,       # k_steel [W/(m·K)]
    3.73e6,     # rho_c_steel [J/(m³·K)]
    0.26,       # k_insulating [W/(m·K)]
    1.96e6,     # rho_c_insulating [J/(m³·K)]
    1.0,        # k_backfill [W/(m·K)]
    1.0         # rho_c_backfill [J/(m³·K)]
)

For layered rock formations, use StratifiedMaterialProperties instead.

Borehole Geometry

Define the coaxial borehole heat exchanger geometry:

borehole = Borehole(
    0.0,             # xc - x-coordinate of center [m]
    0.0,             # yc - y-coordinate of center [m]
    2000.0,          # h - borehole depth [m]
    0.0381,          # r_inner - inner pipe radius [m]
    0.01,            # t_inner - inner pipe wall thickness [m]
    0.0889,          # r_outer - outer pipe inner radius [m]
    0.01,            # t_outer - outer pipe wall thickness [m]
    0.0989,          # r_backfill - backfill outer radius [m]
    10.0,            # ṁ - mass flow rate [kg/s]
    500.0            # insulation_depth - depth of inner pipe insulation [m]
)

boreholes = (borehole,)  # tuple of all boreholes

Grid Setup

Create adaptive grids that are fine near the borehole and coarse far away:

# Domain boundaries
xmin, xmax = -100.0, 100.0
ymin, ymax = -100.0, 100.0
zmin, zmax = 0.0, 2200.0

# Grid parameters
dx_fine = 0.0025      # fine spacing near borehole [m]
growth_factor = 1.3   # geometric growth rate
dx_max = 10.0         # maximum spacing far from borehole [m]
dz = 100.0            # vertical spacing [m]

gridx = create_adaptive_grid_1d(
    xmin=xmin, xmax=xmax,
    dx_fine=dx_fine, growth_factor=growth_factor, dx_max=dx_max,
    boreholes=boreholes, backend=backend, direction=:x
)

gridy = create_adaptive_grid_1d(
    xmin=ymin, xmax=ymax,
    dx_fine=dx_fine, growth_factor=growth_factor, dx_max=dx_max,
    boreholes=boreholes, backend=backend, direction=:y
)

gridz = create_uniform_gridz_with_borehole_depths(
    zmin=zmin, zmax=zmax, dz=dz,
    boreholes=boreholes, backend=backend
)

println("Grid size: $(length(gridx)) x $(length(gridy)) x $(length(gridz))")

Initial Condition

Set the initial temperature field with a geothermal gradient:

T0 = initial_condition_thermal_gradient(
    backend, gridx, gridy, gridz;
    T_surface=10.0,    # surface temperature [°C]
    gradient=0.03      # thermal gradient [K/m]
)

Inlet Model

Choose how the inlet temperature is determined:

# Option 1: Constant inlet temperature
inlet_model = ConstantInlet(20.0)  # 20°C

# Option 2: Heat exchanger (inlet = outlet - ΔT)
# inlet_model = HeatExchangerInlet(5.0)  # ΔT = 5K

Create Cache and Solve

cache = create_cache(
    backend=backend,
    gridx=gridx,
    gridy=gridy,
    gridz=gridz,
    materials=materials,
    boreholes=boreholes,
    inlet_model=inlet_model
)

# Time span: simulate for 1 day
tspan = (0.0, 3600.0 * 24.0)

prob = ODEProblem(rhs_diffusion_z!, T0, tspan, cache)

# Callback handles ADI (horizontal diffusion) and advection, plus saves solutions
# This is required - without it, only vertical diffusion is computed
callback, saved_values = get_simulation_callback(
    saveat=range(tspan..., 10),  # times to save solutions
    print_every_n=100            # print progress every N steps
)

# Solve with ROCK2 (stabilized explicit Runge-Kutta)
Δt = 80.0  # time step [s]

solve(
    prob,
    ROCK2(max_stages=100, eigen_est=eigen_estimator),
    save_everystep=false,
    callback=callback,
    adaptive=false,
    dt=Δt,
    maxiters=Int(1e10)
)

Checkpointing and Restart

For long-running simulations that may be interrupted, pass checkpoint_dir and checkpoint_id to enable fault-tolerant checkpointing. Use prepare_restart at the top and reload_snapshots! at the bottom, and the same script can be resubmitted without changes:

checkpoint_dir = joinpath(@__DIR__, "simulation_data")
checkpoint_id = splitext(basename(@__FILE__))[1]

T0_fresh = initial_condition_thermal_gradient(backend, gridx, gridy, gridz;
    T_surface=10.0, gradient=0.03)
tspan_full = (0.0, 3600.0 * 24 * 365 * 20)
saveat_full = range(tspan_full..., 21)

# Loads checkpoint if available, otherwise passes through unchanged
T0, tspan, saveat = prepare_restart(T0_fresh, tspan_full, saveat_full;
    checkpoint_dir=checkpoint_dir, checkpoint_id=checkpoint_id, backend=backend)

prob = ODEProblem(rhs_diffusion_z!, T0, tspan, cache)

callback, saved_values = get_simulation_callback(
    saveat=saveat,
    checkpoint_dir=checkpoint_dir,
    checkpoint_id=checkpoint_id,
    checkpoint_every_n=500_000
)

solve(prob, ROCK2(max_stages=100, eigen_est=eigen_estimator),
    save_everystep=false, callback=callback, adaptive=false, dt=80.0, maxiters=Int(1e10))

# Assemble full snapshot history from disk (across all restart cycles)
reload_snapshots!(saved_values, checkpoint_dir, checkpoint_id)

Use distinct checkpoint_id values when multiple simulations share the same directory. See the API Reference for details on all keyword arguments.

GPU Acceleration

For GPU acceleration, simply change the backend:

using CUDA: CUDABackend

backend = CUDABackend()

# All other code remains the same!

Well Arrays

To simulate multiple boreholes, create a tuple of Borehole objects:

# 3x3 array with 30m spacing
boreholes = tuple(
    (Borehole(
        xc, yc,          # center coordinates
        2000.0,          # depth
        # ... other parameters
    ) for xc in [-30.0, 0.0, 30.0], yc in [-30.0, 0.0, 30.0])...
)

Next Steps

  • See the API Reference for detailed documentation
  • Check the examples/ folder:
    • single_borehole_basic.jl - Single borehole simulation
    • well_array_basic.jl - Multi-well array simulation
  • See the reproducibility repository for advanced examples reproducing results from literature (Li et al. (2021), Hu et al. (2020), Brown et al. (2023))