Soaring-Sim: A simulation of gliders soaring in a turbulent wind model

This repository contains a C++ code for simulating gliders that are soaring in turbulent wind field. The gliders measure the wind as they fly, and then control how they move in order to gain height. This code has been used to study how gliders—such as birds, or mechanical gliders—can navigate in the atmosphere to most efficiently stay aloft. This code is associated with the following scientific publication:

• Danyun He, Gautam Reddy, and Chris H. Rycroft, Energy-positive soaring using transient turbulent fluctuations, arXiv:2304.05983 (2024).

In particular, this publication shows that for realistic parameters, a glider in an idealized model of isotropic turbulence can continually gain height and extract energy.

Background

Soaring birds harvest energy by strategically gliding through atmospheric flows. For soaring birds like herring gulls and albatrosses, gliding uses oxygen at a rate 30% lower than flapping [1,2], which is important to make long-distance migration feasible [3]. In the absence of wind a glider will continuously lose energy due to drag, and will sink at a constant rate. Therefore, to compensate for drag, gliders will (A) actively localize at updrafts or (B) align themselves to take advantage of rapid gusts and wind shear. Birds take advantage of both of these mechanisms, which are known as thermal soaring and dynamic soaring respectively. Thermal soaring relies on ascending currents that are created by convection in the atmospheric boundary layer. Dynamic soaring allows for energy-neutral flight over oceans, via the bird moving through a stable shear layer generated by ocean waves.

Both of these scenarios rely on the formation of relatively stable convective plumes or wind shear. Previous studies have examined navigational strategies that can exploit these flow structures to gain energy. As an example, Reddy et al. used reinforcement learning to learn a glider soaring strategy in a Rayleigh–Bénard flow [4]. This flow has ascending thermals that the gliders are able to locate to gain energy continuously. Later work demonstrated this strategy in the field [5].

Atmospheric flows also contain short-lived turbulent eddies across a range of timescales. In this study, we asked whether it is possible for a glider to continuously gain energy using only these transient turbulent fluctuations, in a situation without any large scale flow features. To do this, we built an idealized model of an isotropic turbulent flow based on work by Fung et al. [6]. The three-dimensional wind field is modeled as a stochastically evolving set of Fourier modes whose statistics match Kolmogorov's theory of turbulence.

We then built an aerodynamic model for a glider moving through this wind field. We considered two different scenarios for how the glider interprets its surroundings: (A) it has complete information of the current state of the wind, or (B) it remembers the wind that it has encountered along its flight trajectory and estimates the state of the wind using Gaussian process regression (GPR).

To navigate through the flow, the glider can control its bank angle at regular intervals. The glider's planning strategy uses the Monte Carlo tree search (MCTS) to assess the energy gain of possible bank angle changes into the future, based on its interpretation of the wind field. The glider chooses the bank angle change to maximize its energy gain, and advances forward for an interval before executing another planning step.

Individual glider trajectories exhibit large variations. Therefore the code is designed to run an ensemble of gliders to extract reliable statistics. A typical run involves considering a number of different turbulent wind fields in sequence. Each turbulent flow field contains a number of gliders that are initialized at random locations within the flow field and operate independently. The code can collect and save a wide range of data and statistics about each run.

Overview of the repository

Directory	Contents
src	Library source and main programs, with a description of the C++ classes
tests	Test and analysis programs, and how to run them
config	Build configuration templates
sims	Sample simulation configuration files

Compiling the code

The code is written in C++ and uses the OpenMP library for multithreading. It has been tested on Linux, MacOS, and Windows via the Windows Subsystem for Linux (WSL). The code depends on three open source software libraries:

• The code makes use of the FFTW library for efficient computations of the fast Fourier transform.

• The code uses the GNU Scientific Library (GSL) for high-quality Gaussian random number generation.

• The Gaussian process regression routines involve numerical linear algebra. This is handled using custom routines based on the Woodbury formula. The code also makes use of LAPACK for diagnostic purposes, but this dependency can be removed. LAPACK is often installed by default on many new operating systems.

All three of these libraries are available via software package management systems, such as those available by default in Linux and in WSL. On the Mac, they are available via MacPorts and Homebrew.

To compile the code, it is necessary to create a common configuration file called config.mk in the top-level directory. Several templates are provided in the config directory. To use, copy one of the templates into the top-level directory. From the soaring-sim directory, on a Linux computer, type

cp config/config.mk.linux config.mk

On a Mac using GCC 15 installed via MacPorts, type

cp config/config.mk.mac_mp config.mk

On a Mac using GCC installed via Homebrew, type

cp config/config.mk.mac_hb config.mk

After this, the code can be compiled by typing

make

This will build several executables such as soar and unpack.

Running an example simulation

The main program for running the glider soaring simulation is called soar. It requires a text file that contains the simulation configuration. Several samples are provided in the sims directory. A good place to begin is the lfd64small.cfg file for simulating full information gliders in a turbulent wind field with a 64 s relaxation timescale.

The simulation can be run by typing

OMP_NUM_THREADS=<n> ./soar sims/lfd64small.cfg

where <n> is replaced with the number of OpenMP threads to use. The example may take approximately ten minutes to run, depending on hardware and the number of threads used. The following sections describe the example configuration file, and output data that is created from it.

The example simulation configuration file

The configuration file is divided into sections. The file consists of keywords followed by one or more parameters. A pound symbol in front of any line signifies that it is a comment. The first section sets up the physical scales, specified in SI units, which are used to scale all of the simulation variables:

# Physical scales in SI units
l_phys              10
g_phys              9.80665
nu_phys             1.48e-5

This sets the simulation length scale to be 10 m and the gravitational acceleration to be 9.80665 m/s², which is the mean value at the Earth's surface. From these two constants, a velocity scale and time scale can be derived as described in the paper. The velocity scale and time scale can also be set using the v_phys and t_phys keywords. Exactly two scales need to be set, from which the other can be derived.

The line for nu_phys contains the kinematic viscosity of the air in units of m²/s. It has no effect on the simulation, but is included for completeness since it can be used to compute the effective Reynolds number.

The next section contains several global parameters:

# Simulation setup
num_trials          4
gliders_per_trial   8
duration_phys       200
random_seed         100

num_trials sets the total number of random wind fields that will be simulated. gliders_per_trial sets the number of gliders that are simulated within each wind field. The product of these two numbers sets the total number of gliders to be considered, from which statistical averages will be calculated.

The total simulation duration can be set in simulation units using the duration keyword. Alternatively, a _phys suffix can be appended to specify the simulation duration in SI units. The random_seed value is used as a base for initializing the GSL random number generators.

The next section initializes the glider parameters, along with details of numerical integration:

# Lift and drag constants
c_L                 1
c_D                 0.0666666666666666666667
ctl_duration_phys   1
bank_angle_step     10
bank_angle_range    -4 4
bank_angle_control  -1 1
integration_type    improv_e
glider_ts_pad       0.5

c_L and c_D set the dimensionless lift and drag coefficients that are described in the paper. Here, the lift-to-drag ratio is set to 15, which is a typical value for soaring birds. ctl_duration sets the duration of time between the glider controlling its bank angle.

bank_angle_step sets the increment in bank angles for the glider in degrees. bank_angle_range sets the number of discrete bank angle steps that the glider can achieve. Here, the bank angle step is set to 10°, and the full step range is -40°, -30°, …, 40°. bank_angle_control sets the allowable number of step changes at each control time point. Here, the -1 to 1 range corresponds to three possible controls of -10°, 0°, and +10°.

integration_type sets the numerical integration method used for the glider ODE model. It can be set to either euler for the forward Euler method, or improv_e for the improved Euler method. Both integration methods are explicit, and therefore have a restriction on the numerical timestep. The selection of the timestep and the usage of glider_ts_pad is described in the next section.

The next section sets up the turbulent wind field:

# Wind field dimensions
wf_modes            64
box_size            50
wind_rms_phys       1.73205080756888
wind_C_phys         64
wind_ts_pad         0.2

wf_modes contains the number of Fourier modes in each dimension. box_size sets the side length of the simulation box, specified in simulation units.

Each component of the wind field in the three coordinates evolves independently. In this example, the root-mean-squared (RMS) of each velocity component is set to 1 m/s. wind_rms_phys sets the RMS wind velocity magnitude in physical units, which is √3 × 1 m/s = 1.732 m/s in this case. wind_C_phys sets the timescale for the temporal correlations in physical units. If the _phys suffix is dropped from either of these keywords, the values will be interpreted as simulation units instead.

During the simulation, the Fourier modes of the wind field are updated according to a stochastic differential equation (SDE) following an Ornstein–Uhlenbeck (OU) process. This numerical integration has a timestep restriction. The computation of the timestep and the usage of wind_ts_pad is described in the following section. wind_C_phys and wind_C both accept the frozen keyword if the wind field is frozen in time.

The next section sets up the wind prediction model:

# Wind prediction model
wind_model          full_lanczos2
hermite_segments    20

The wind_model keyword chooses the method that the glider uses to predict the wind. The full_cubic option corresponds to a glider with complete information about wind. This is achieved by taking the fast Fourier transform of the wind field, and then performing Lanczos2 interpolation [7] to evaluate the wind at an arbitrary location. Two alternative are available:

• The full_linear option uses trilinear interpolation, which is faster to evaluate but is slightly less accurate.
• The full_cubic [<a>] option uses tricubic interpolation based on Keys' formula [8]. This formula has a free parameter a that controls the shape of the interpolation kernel. The program defaults to a value of a = -0.84, which was determined empirically to give close results to the exact wind field. This parameter can be overridden.

If the wind field is time-varying, then complete-information models must also predict the wind up to a duration T into the future, by assuming that the Fourier modes will undergo mean reversion. To do this, the model divides the interval T into a number of segments specified by hermite_segments. It computes the expected wind field and its temporal derivative at each segment boundary, and uses cubic Hermite interpolation in time to evaluate the expected wind at any point in the future.

An alternative wind prediction model is to use Gaussian process regression (GPR), so that the glider only measures the wind that it has encountered during flight. A typical configuration of this is as follows:

# Wind prediction model (GPR)
#gpr_memory          50
#gpr_samp_per_ctl    2
#gpr_k_rt_param      160 160 1.25

gpr_memory sets the number of wind measurements that are retained in memory. gpr_samp_per_ctl sets the number of times that the wind is sampled per control duration. During these computations, it is necessary to evaluate the wind correlation function K(r,t) very frequently. Prior to starting the glider simulation, the code evaluates the kernel on a grid in r and t. This allows the kernel to be rapidly approximated using bilinear interpolation. The first two arguments to gpr_k_rt_param set the interpolation grid size in r and t, respectively. The code computes the maximum values of r and t that will be encountered in the simulation to decide how large the interpolation grid should be. These bounds are padded by the multiplicative factor given in the third argument.

The next section of the configuration sets up the planning strategy:

# Planning
planning            mcts
n_mcts              1000
mcts_depth          12
mcts_ex_fac         1

The mcts option enables the Monte Carlo tree search for planning. n_mcts sets the total number of MCTS trials. mcts_depth sets the depth of the MCTS tree, which corresponds to the number of actions to make. mcts_ex_fac controls the exploration factor, with higher values corresponding to a higher probability of choosing less favorable branches of the search tree. The planning keyword will also accept a random option for a glider that chooses among its available moves at each step with equal probability.

The final section of the configuration controls the simulation output:

# Output types
output_per_ctl      4
snapshots           glider_mb
summary             climb_stats

output_per_ctl sets the number of output snapshots to perform between each control duration. For example, if this parameter is set to 4, and the control duration is 1 s, the simulation output will be saved every 0.25 s. The snapshots keyword is followed by a list of the types of output snapshots to save, with four different options:

• glider_xyz outputs the glider positions in a text file.
• glider_mb outputs a minimal binary file that contains the glider positions that can be unpacked later. Binary output is more space-efficient than text.
• glider_lb outputs a larger binary file of glider information, including positions and wind information.
• wind_modes outputs the full three-dimensional grid of Fourier mode coefficients for the wind field.

The summary keyword is followed by a list of summary statistics to output. climb_stats outputs a file with statistics about the following four quantities:

1. glider height,
2. glider height rate-of-change,
3. glider energy,
4. glider energy rate-of-change.

For each quantity, the minimum, maximum, mean, and standard deviation across the entire data set of gliders are stored. The primary contribution to the energy comes from the potential energy of the glider. Therefore quantities 1 & 3 are similar, and quantities 2 & 4 are similar. Quantities 3 & 4 differ slightly, because they also contain contributions from kinetic energy.

Output data

For a given configuration file .cfg, all output data will be stored within a directory called .odr, where odr stands for “output dirctory”. When the example simulation is run, the data is therefore stored in lfd64small.odr.

The snapshots glider_mb option in the test simulation creates binary files glider_mb. where N varies from 0 to 3 for the four different wind fields that were simulated. Each file contains the all of the glider trajectories in a binary format. A utility called unpack is provided, which can output a particular glider trajectory as plain text, with each line listing four numbers for the time and (x,y,z) position of the glider, all given in simulation units. The utility has the syntax

./unpack [-w] [-z] <output_dir> <trial> <glider> <output_file>

Hence, to output the glider 0 in wind field 2 for the example simulation, the utility can be run with

./unpack lfd64small.cfg 2 0 gtraj.txt

The three-dimensional trajectory can be plotted using the freeware plotting program Gnuplot with the command

splot 'gtraj.txt' u 2:3:4 w l

The height of the glider as a function of time can be plotted using the command

splot 'gtraj.txt' u 1:4 w l

If the -w option is provided to the unpack utility, then the wind components at the glider's position are also outputted; this requires that the larger binary files were stored using the snapshots glider_lb option. If the -z option is provided, then the glider's initial height is subtracted out; this option can be useful to compare the gained height for many gliders.

If the output_file option is specified as -, then the utility outputs the data to standard output. This makes it possible to pipe data directly into Gnuplot without storing a separate file:

splot '<./unpack lfd64small.cfg 2 0 -' u 2:3:4 w l

All eight gliders in wind field 2 can be visualized with the command

splot for[n=0:7] sprintf('<./unpack lfd64small.cfg 2 %d -',n) \
      u 2:3:4 w l t sprintf("Glider %d",n)

Additional configuration keywords

The input files can contain a number of keywords for enabling further functionality or performing diagnostic tests:

• box_size_xyz accepts three arguments for the wind field size in the three coordinate directions, specified in simulation units. It can be used to simulate gliders in a non-cubic box.

• wf_modes_xyz accepts three arguments for the number of Fourier modes in the wind field in the three coordinate directions.

• For a full information glider in a non-frozen wind field, bypass_mr switches off the mean reversion in the glider's wind model. The glider assumes that the current wind will continue as is, and it will not capture how on average the wind will mean revert to zero. This is a useful diagnostic check of the simulation, since it shows that in rapidly time-evolving wind fields, mean-reversion substantially enhances the glider's performance.

• When the simulation terminates, the summary checksums option makes the code check that the kernel inverses that are efficiently computed using the Woodbury formula do indeed match a full LAPACK calculation. The code computes the Frobenius norm between the Woodbury-formula-based inverse and the LAPACK inverse, and scales it by the matrix size to give a measure of error per matrix entry.

• gpr_full_compute forces the Gaussian process regression class to bypass the Woodbury formula and instead always do a full matrix inverse via LAPACK.

• The keywords gpr_rs, gpr_ts can be used to set the number of gridpoints in space and time, respectively, for the precomputed table of kernel values K(r,t) in Gaussian process regression. gpr_dist_pad sets the padding factor to apply to the base estimate of the range of r values needed in the table. The gpr_k_r_param keyword requires two arguments to set the number of gridpoints in space, and the padding factor. It is intended to be used for frozen wind field simulation where the K is a function of r only.

• In extremely rare instances it is possible for two rows in the kernel matrix to be near-linearly-dependent. This happens most often in a frozen wind field when a glider circles back to a position that it had visited before, leading to two measurements at near-identical positions, and two rows with near-identical kernel entries. This can cause a loss of accuracy in the Woodbury formula. gpr_ker_tol sets a numerical tolerance for an experimental option that identifies these cases and skips over the additional measurement, to avoid the near-singular behavior.

• The path_interval option makes the simulation output binary files containing all of the MCTS paths for all of the gliders. This can be useful for diagnostic purposes. The option requires a single argument that sets the frequency at which the paths are stored. Note that this option can result in a large amount of data being generated.

• The wind_cor option computes the correlation coefficient between the real wind field and the one that GPR predicts, at various times into the future. It is based on GPR predictions along all of the MCTS paths that the gliders consider during planning. These statistics are only valid in a frozen wind field, where it is known what the real wind field will be in the future. The option requires one argument that sets a time (in simulation units) at which to start collecting data. The variant wind_cor_phys sets the time in physical units instead.

• The mcts_info option outputs statistical information about the MCTS. It outputs the average number of nodes in the search tree at different depths. The option requires one argument that sets a time in simulation units at which to start collecting data. The variant mcts_info_phys sets the time in physical units instead.

• wind_alpha sets the prefactor in the wind energy spectrum. It is an alternative to specifying the RMS wind velocity magnitude.

• summary en_components outputs information about the different components to the glider's energy gain.

• wind_slice outputs a two-dimensional cross-section of the wind field. It accepts three options: (1) the slice direction as x, y, or z; (2) the wind component as wx, wy, or wz; and (3) the grid index along the slice direction at which to output the field. The field is stored in a binary format that can be read by the freeware plotting program Gnuplot. See the utils-gp for details on the binary format.

Timestep selection

The code automatically computes the timestep to use based on two analytical calculations of timestep restrictions. The first restriction comes from the stochastic integration of the Fourier modes in the turbulence model. Each mode follows an Ornstein–Uhlenbeck (OU) process as described in the paper. The OU process features a mean-reversion term that places a restriction on the timestep in order for the integration to be stable. The tightest restriction comes from the mode with the smallest wavelength. Define this restriction to be Δt₁.

The glider model described in the paper also has a timestep restriction for stability. Here, the code computes a restriction Δt₂ based on several simplifying assumptions of the glider model:

1. the glider is assumed to be flying straight with no banking,
2. forward Euler integration is used,
3. there is no coupling to the wind field.

Even with these assumptions Δt₂ remains a good estimate of the timestep restriction that will differ from the true one by a constant factor.

After this, a proposed simulation timestep is chosen as Δt_p = min(p₁ Δt₁, p₂ Δt₂), where p₁ and p₂ are the padding factors of wind_ts_pad and glider_ts_pad described in the previous section. In the case when a frozen wind field is used, the first timestep is ignored and Δt_p =p₂ Δt₂.

Finally, the timestep Δt_p is adjusted downwards slightly to Δt_a so that an integer number of timesteps will exactly cover the integration intervals between output snapshots, and MCTS & GPR computations.

Code structure, additional utilities, and further tests

The code is structured around several C++ classes provided in the src directory, which are described in src/README.md. This directory contains the main program soar as well as three utilities unpack, get_paths, and climb_rate.

The tests directory contains several additional test programs, including some that were used to generate data for the paper. The programs are documented in tests/README.md.

Bibliography

1. R. V. Baudinette and K. Schmidt-Nielsen, Energy cost of gliding flight in herring gulls, Nature 248, 83–84 (1974). doi:10.1038/248083b0

2. Kentaro Q. Sakamoto, Akinori Takahashi, Takashi Iwata, Takashi Yamamoto, Maki Yamamoto, and Philip N. Trathan, Heart rate and estimated energy expenditure of flapping and gliding in black-browed albatrosses, 216, 3175–3182 (2013). doi:10.1242/jeb.079905

3. Vance A. Tucker, Metabolism during flight in the laughing gull, Larus atricilla, Am. J. Physiol. 222, 237–245 (1972).

4. Gautam Reddy, Antonio Celani, Terrence Sejnowski, and Massimo Vergassola, Learning to soar in turbulent environments, Proc. Natl. Acad. Sci. 113, E4877–E4884 (2016). doi:10.1073/pnas.1606075113

5. Gautam Reddy, Jerome Wong-Ng, Antonio Celani, Terrence Sejnowski, and Massimo Vergassola, Glider soaring via reinforcement learning in the field, Nature 562, 236–239 (2018). doi:10.1038/s41586-018-0533-0

6. Jimmy C. H. Fung, Julian C. R. Hunt, Nadeem Malik, and R. J. Perkins, Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes, J. Fluid Mech. 236, 281–318 (1992). doi:10.1017/S0022112092001423

7. Ken Turkowski, Filters for Common Resampling Tasks, in Graphics Gems (ed. Andrew S. Glassner), Academic Press, 1990.

8. Robert G. Keys, Cubic convolution interpolation for digital image processing, IEEE Trans. Acoust., Speech, Signal Process. 29, 1153–1160 (1981). doi:10.1109/TASSP.1981.1163711