Soaring-Sim:

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 turbulunce 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 oxgyen 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 a idealized model of an isotropic turbulent flow based on work by Fung et al. [6]. The three-dimensional wind field is modeled as a combination of stochastically evolving set of Fourier modes whose statistics match Kolmogorov's theory of turbulence.

We then built a aerodynamical 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 encounted along its flight trajectory an 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.

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 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 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 parent directory. Several templates are provided in the config directory. To use, copy one of the templates into the parent directory. From the soaring-sim directory, on a Linux computer, type

cp config/config.mk.linux ../config.mk

On a Mac using GCC 14 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.

Configuring and running a 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 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 use 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.80655 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 stastical 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 set the number of discrete bank angle steps that the glider can achieve. Here, the bank angle 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&deg, 0°, and +10&deg.

integration_type set 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.

wind_rms sets the root-mean-squared (RMS) wind velocity in simulation units, and wind_C_phys sets the timescale for the temporal correlations. 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 is 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 also accepts 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_cubic
hermite_segments    20

The wind_model keyword chooses 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 tricubic interpolation to evaluate the wind at an arbitrary location. The alternative full_linear option uses trilinear interpolation, which is faster to evaluate but is slightly less accurate.

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 is has encountered during flight. A typical configration 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 probabilities of choosing less favorable branches of the search tree.

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 the glider positions in a binary file that can be unpacked later. Binary output is more space-efficient than text.
• glider_fb outputs full 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 quantity:

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 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.

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 Orstein–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_p so that an integer number of timesteps will exactly cover the integration intervals between output snapshots, and MCTS & GPR computations.

Code structure

The code is structured around several C++ classes:

• glider – This class represents the state of a single glider, including its position, velocity, and bank angle. The class contains routines for integrating the glider state via the forward Euler method or improved Euler method.

• glider_model - This class contains all of the glider model parameters. It also calculates several trigonometrical tables of bank angles that are used in the differential equations for the glider. Only one instance of this class is needed, and can be shared amongst many glider instances.

• turb_fluid – This class represents a model of three-dimensional turbuluent wind field made up of stochastically evolving set of a Fourier modes. It contains routines for stepping the modes forward in time, as well as evaluating the wind at one or several locations.

• turb_fluid_grid – This class is derived from turb_fluid. It contains additional functionality to evaluate the wind on a full three-dimensional grid using the FFTW library. It also contains routines for rapidly evaluating the wind at any location using either trilinear or tricubic interpolation.

• turb_fluid_grid_mr - This class is derived from turb_fluid_grid. It contains additional functionality to estimate the wind field into the future assuming that the Fourier modes undergo reversion to the mean, following the OU process.

• kernel_func, kernel_rt, & kernel_r – These classes are used to evaluate the kernel function used in the Gaussian process regression (GPR). The kernel function is derived from the correlations in the wind field. Since the kernel function needs to be called often, the class creates a lookup table of values,

• fileinfo – This class parses the text configuration files and reads in all of the required parameters. It also performs calculations for initializing the simulation, such as calculating the timestep.

• cli_stats & mti_stats – Small structures for computing the mean, standard deviation, minimum, and maximum of a group of numbers. cli_stats is used for collecting statistics on climb rates, and mti_stats is used for collecting statistics on the MCTS.

• soaring_sim – This is the main class to performing the soaring simulations, and makes use of all the other classes as components. This class also contains output routines.

The code is commented in the style of Doxygen, where each class and function has a special comment block beginning with /** that describes it.

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