rnatm

rnatm_octree.c (52145B)

---

 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique
  * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace
  * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris
  * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2022, 2023, 2025 Observatoire de Paris
  * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne
  * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin
  * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #define _POSIX_C_SOURCE 200112L /* lround */
 
 #include "rnatm_c.h"
 #include "rnatm_log.h"
 #include "rnatm_voxel_partition.h"
 #include "rnatm_octrees_storage.h"
 
 #include <star/sars.h>
 #include <star/sck.h>
 #include <star/suvm.h>
 #include <star/svx.h>
 
 #include <rsys/clock_time.h>
 #include <rsys/condition.h>
 #include <rsys/cstr.h>
 #include <rsys/double3.h>
 #include <rsys/float4.h>
 #include <rsys/math.h>
 #include <rsys/morton.h>
 #include <rsys/mutex.h>
 #include <rsys/rsys.h>
 
 #include <math.h> /* lround */
 #include <omp.h>
 
 #define VOXELIZE_MSG "voxelize atmosphere: %3d%%\r"
 
 /* Structure used to synchronise the voxelization and the build threads */
 struct build_sync {
   struct mutex* mutex;
   struct cond* cond;
   size_t ibatch; /* Index of the batch currently consumed by the build thread */
 };
 #define BUILD_SYNC_NULL__ {NULL, NULL, 0}
 static const struct build_sync BUILD_SYNC_NULL = BUILD_SYNC_NULL__;
 
 struct build_octree_context {
   struct pool* pool;
   struct partition* part; /* Current partition */
   size_t iitem; /* Index of the item to handle in the partition */
 
   /* Optical thickness threshold criteria for merging voxels */
   double tau_threshold;
 };
 #define BUILD_OCTREE_CONTEXT_NULL__ {NULL, NULL, 0, 0}
 static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
   BUILD_OCTREE_CONTEXT_NULL__;
 
 struct radcoefs {
   float ka_min, ka_max;
   float ks_min, ks_max;
 };
 #define RADCOEFS_NULL__ {                                                      \
   FLT_MAX, -FLT_MAX,                                                           \
   FLT_MAX, -FLT_MAX,                                                           \
   FLT_MAX, -FLT_MAX                                                            \
 }
 
 /* Generate the dynamic array of radiative coefficient lists */
 #define DARRAY_NAME radcoef_list
 #define DARRAY_DATA struct radcoefs*
 #include <rsys/dynamic_array.h>
 
 /* Generate the dynmaic array of partitions */
 #define DARRAY_NAME partition
 #define DARRAY_DATA struct partition*
 #include <rsys/dynamic_array.h>
 
 struct voxelize_batch_args {
   /* Working data structure */
   struct darray_size_t_list* per_thread_tetra_list;
   struct darray_radcoef_list* per_thread_radcoef_list;
   struct darray_partition* per_thread_dummy_partition;
   struct pool* pool;
 
   size_t part_def; /* Partition definition */
   size_t nparts[3]; /* #partitions required to cover the entire grid */
   size_t nparts_adjusted; /* #partitions allowing their indexing by morton id */
   size_t nparts_overall; /* Total number of partitions to voxelize */
   size_t nparts_voxelized; /* #partitions already voxelized */
 
   /* AABB of the atmosphere */
   double atm_low[3];
   double atm_upp[3];
 
   double vxsz[3]; /* Size of a voxel */
 
   struct accel_struct* accel_structs;
   size_t batch_size; /* #acceleration structures */
 };
 #define VOXELIZE_BATCH_ARGS_NULL__ {                                           \
   NULL, /* Per thread tetra list */                                            \
   NULL, /* Per thread radiative coefs */                                       \
   NULL, /* Per thread dummy partition */                                       \
   NULL, /* Partition pool */                                                   \
                                                                                \
   0, /* Partition definition */                                                \
   {0,0,0}, /* #partitions to cover the entire grid */                          \
   0, /* #partitions adjusted wrt morton indexing */                            \
   0, /* Total number of partitions to voxelize */                              \
   0, /* #partitions already voxelized */                                       \
                                                                                \
   /* AABB of the atmosphere */                                                 \
   { DBL_MAX, DBL_MAX, DBL_MAX},                                                \
   {-DBL_MAX,-DBL_MAX,-DBL_MAX},                                                \
                                                                                \
   {0,0,0}, /* Size of a voxel */                                               \
                                                                                \
   NULL, /* Spectral items */                                                   \
   0 /* Batch size */                                                           \
 }
 static const struct voxelize_batch_args VOXELIZE_BATCH_ARGS_NULL =
   VOXELIZE_BATCH_ARGS_NULL__;
 
 struct voxelize_args {
   struct darray_size_t* tetra_ids;
   struct radcoefs* per_item_radcoefs;
 
   struct partition* part;
   struct partition* part_dummy;
   size_t part_def; /* Partition definition */
 
   /* AABB of the partition */
   double part_low[3];
   double part_upp[3];
 
   double vxsz[3]; /* Size of a voxel */
 
   const struct accel_struct* accel_structs;
   size_t batch_size; /* #spectral items */
 
 };
 #define VOXELIZE_ARGS_NULL__ {                                                 \
   NULL, /* Tetra ids */                                                        \
   NULL, /* Radiative coefs */                                                  \
                                                                                \
   NULL, /* Partition */                                                        \
   NULL, /* Dummy partition */                                                  \
   0, /* Partition definition */                                                \
                                                                                \
   /* AABB of the partition */                                                  \
   { DBL_MAX, DBL_MAX, DBL_MAX},                                                \
   {-DBL_MAX,-DBL_MAX,-DBL_MAX},                                                \
                                                                                \
   {0,0,0}, /* Size of a voxel */                                               \
                                                                                \
   NULL, /* Spectral items */                                                   \
   0 /* Batch size */                                                           \
 }
 static const struct voxelize_args VOXELIZE_ARGS_NULL = VOXELIZE_ARGS_NULL__;
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static INLINE res_T
 check_voxelize_batch_args
   (struct rnatm* atm,
    const struct voxelize_batch_args* args)
 {
   size_t i;
   if(!args
   || !args->per_thread_tetra_list
   || !args->per_thread_radcoef_list
   || !args->pool
   || !IS_POW2(args->part_def)
   || !args->nparts[0]
   || !args->nparts[1]
   || !args->nparts[2]
   || !IS_POW2(args->nparts_adjusted)
   || args->nparts_adjusted < args->nparts[0]*args->nparts[1]*args->nparts[2]
   || args->atm_low[0] >= args->atm_upp[0]
   || args->atm_low[1] >= args->atm_upp[1]
   || args->atm_low[2] >= args->atm_upp[2]
   || !args->vxsz[0]
   || !args->vxsz[1]
   || !args->vxsz[2]
   || !args->accel_structs
   || !args->batch_size) {
     return RES_BAD_ARG;
   }
   if(atm->nthreads != darray_size_t_list_size_get(args->per_thread_tetra_list)
   || atm->nthreads != darray_radcoef_list_size_get(args->per_thread_radcoef_list)) {
     return RES_BAD_ARG;
   }
   FOR_EACH(i, 0, atm->nthreads) {
     if(!darray_radcoef_list_cdata_get(args->per_thread_radcoef_list)[i]) {
       return RES_BAD_ARG;
     }
   }
   return RES_OK;
 }
 
 static INLINE res_T
 check_voxelize_args(const struct voxelize_args* args)
 {
   if(!args->tetra_ids
   || !args->per_item_radcoefs
   || args->part_low[0] >= args->part_upp[0]
   || args->part_low[1] >= args->part_upp[1]
   || args->part_low[2] >= args->part_upp[2]
   || !IS_POW2(args->part_def)
   || !args->vxsz[0]
   || !args->vxsz[1]
   || !args->vxsz[2]
   || !args->accel_structs
   || !args->batch_size) {
     return RES_BAD_ARG;
   }
   return RES_OK;
 }
 
 static INLINE res_T
 check_trace_ray_args
   (const struct rnatm* atm,
    const struct rnatm_trace_ray_args* args)
 {
   size_t n;
   if(!args) return RES_BAD_ARG;
 
   /* Invalid band index */
   n = darray_band_size_get(&atm->bands) - 1;
   if(args->iband < darray_band_cdata_get(&atm->bands)[0].index
   || args->iband > darray_band_cdata_get(&atm->bands)[n].index)
     return RES_BAD_ARG;
 
   return RES_OK;
 }
 
 static void
 build_sync_release(struct build_sync* sync)
 {
   if(sync->mutex) mutex_destroy(sync->mutex);
   if(sync->cond) cond_destroy(sync->cond);
 }
 
 static res_T
 build_sync_init(struct build_sync* sync)
 {
   res_T res = RES_OK;
   ASSERT(sync);
 
   memset(sync, 0, sizeof(*sync));
 
   sync->mutex = mutex_create();
   if(!sync->mutex) { res = RES_UNKNOWN_ERR; goto error; }
   sync->cond = cond_create();
   if(!sync->cond) { res = RES_UNKNOWN_ERR; goto error; }
   sync->ibatch = 0;
 
 exit:
   return res;
 error:
   build_sync_release(sync);
   goto exit;
 }
 
 static FINLINE unsigned
 round_pow2(const unsigned val)
 {
   const unsigned next_pow2 = (unsigned)round_up_pow2(val);
   if(val == 0 || next_pow2 - val <= next_pow2/4) {
     return next_pow2;
   } else {
     return next_pow2/2;
   }
 }
 
 /* Return the number of quadrature points for the range of bands overlapped by
  * the input spectral range */
 static INLINE size_t
 compute_nquad_pts(const struct rnatm* atm)
 {
   size_t nquad_pts = 0;
   size_t iband = 0;
   ASSERT(atm);
 
   FOR_EACH(iband, 0, darray_band_size_get(&atm->bands)) {
     nquad_pts += darray_band_cdata_get(&atm->bands)[iband].nquad_pts;
   }
   return nquad_pts;
 }
 
 static res_T
 setup_accel_structs(struct rnatm* atm)
 {
   struct accel_struct* accel_structs = NULL;
   size_t istruct;
   size_t naccel_structs;
   size_t i;
   res_T res = RES_OK;
   ASSERT(atm);
 
   naccel_structs = compute_nquad_pts(atm);
 
   res = darray_accel_struct_resize(&atm->accel_structs, naccel_structs);
   if(res != RES_OK) {
     log_err(atm,
       "%s: failed to allocate the list of acceleration structures -- %s",
       FUNC_NAME, res_to_cstr(res));
     goto error;
   }
   accel_structs = darray_accel_struct_data_get(&atm->accel_structs);
 
   istruct = 0;
   FOR_EACH(i, 0, darray_band_size_get(&atm->bands)) {
     const struct band* band = darray_band_cdata_get(&atm->bands)+i;
     size_t iquad_pt;
 
     FOR_EACH(iquad_pt, 0, band->nquad_pts) {
       ASSERT(istruct < naccel_structs);
       accel_structs[istruct].iband = band->index;
       accel_structs[istruct].iquad_pt = iquad_pt;
       ++istruct;
     }
   }
 
 exit:
   return res;
 error:
   darray_accel_struct_clear(&atm->accel_structs);
   goto exit;
 }
 
 static INLINE void
 register_tetra
   (const struct suvm_primitive* prim,
    const double low[3],
    const double upp[3],
    void* context)
 {
   struct darray_size_t* tetra_ids = context;
   ASSERT(prim && low && upp && context);
   ASSERT(low[0] < upp[0]);
   ASSERT(low[1] < upp[1]);
   ASSERT(low[2] < upp[2]);
   (void)low, (void)upp;
   CHK(darray_size_t_push_back(tetra_ids, &prim->iprim) == RES_OK);
 }
 
 static res_T
 compute_grid_definition(struct rnatm* atm, const struct rnatm_create_args* args)
 {
   double low[3];
   double upp[3];
   double sz[3];
   double sz_max;
   double vxsz;
   unsigned def[3];
   int iaxis_max;
   int iaxis_remain[2];
   int i;
   ASSERT(atm && args);
 
   /* Recover the AABB from the atmosphere. Note that we have already made sure
    * when setting up the meshes of the gas and aerosols that the aerosols are
    * included in the gas. Therefore, the AABB of the gas is the same as the
    * AABB of the atmosphere */
   SUVM(volume_get_aabb(atm->gas.volume, low, upp));
   sz[0] = upp[0] - low[0];
   sz[1] = upp[1] - low[1];
   sz[2] = upp[2] - low[2];
 
   /* Find the axis along which the atmosphere's AABB is greatest */
   sz_max = -DBL_MAX;
   FOR_EACH(i, 0, 3) {
     if(sz[i] > sz_max) { iaxis_max = i; sz_max = sz[i]; }
   }
 
   /* Define the other axes */
   iaxis_remain[0] = (iaxis_max + 1) % 3;
   iaxis_remain[1] = (iaxis_max + 2) % 3;
 
   /* Fix the definition along the maximum axis and calculate the size of a
    * voxel along this axis */
   def[iaxis_max] = round_pow2(args->grid_definition_hint);
   vxsz = sz[iaxis_max] / def[iaxis_max];
 
   /* Calculate the definitions along the remaining 2 axes. First calculates
    * them assuming the voxels are cubic, then rounds them to their nearest
    * power of two */
   def[iaxis_remain[0]] = (unsigned)lround(sz[iaxis_remain[0]]/vxsz);
   def[iaxis_remain[1]] = (unsigned)lround(sz[iaxis_remain[1]]/vxsz);
   def[iaxis_remain[0]] = round_pow2(def[iaxis_remain[0]]);
   def[iaxis_remain[1]] = round_pow2(def[iaxis_remain[1]]);
 
   /* Setup grid definition */
   atm->grid_definition[0] = def[0];
   atm->grid_definition[1] = def[1];
   atm->grid_definition[2] = def[2];
 
   return RES_OK;
 }
 
 /* Compute the radiative coefficients range of the tetrahedron */
 static INLINE void
 setup_tetra_radcoefs
   (struct rnatm* atm,
    const struct suvm_primitive* tetra,
    const size_t cpnt,
    const size_t iband,
    const size_t iquad,
    struct radcoefs* radcoefs)
 {
   float ka[4];
   float ks[4];
   ASSERT(atm && tetra && radcoefs);
   ASSERT(cpnt == RNATM_GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
 
   if(tetra_get_radcoef(atm, tetra, cpnt, iband, iquad, RNATM_RADCOEF_Ka, ka)
   && tetra_get_radcoef(atm, tetra, cpnt, iband, iquad, RNATM_RADCOEF_Ks, ks)) {
     /* Radiative coefficients are defined for this component at the considered
      * spectral band */
     radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
     radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
     radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
     radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
   } else {
     /* No valid radiative coefficients are set: ignore the tetrahedron */
     radcoefs->ka_min = FLT_MAX;
     radcoefs->ka_max =-FLT_MAX;
     radcoefs->ks_min = FLT_MAX;
     radcoefs->ks_max =-FLT_MAX;
   }
 }
 
 static INLINE void
 update_voxel
   (struct rnatm* atm,
    const struct radcoefs* radcoefs,
    struct partition* part,
    const uint64_t vx_mcode,
    const uint64_t vx_iitem)
 {
   float vx_ka_min, vx_ka_max;
   float vx_ks_min, vx_ks_max;
 
   float* vx = NULL;
   ASSERT(atm && radcoefs && part);
   (void)atm;
 
   vx = partition_get_voxel(part, vx_mcode, vx_iitem);
 
   /* Update the range of the radiative coefficients of the voxel */
   vx_ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
   vx_ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
   vx_ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
   vx_ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
   vx_ka_min = MMIN(vx_ka_min, radcoefs->ka_min);
   vx_ka_max = MMAX(vx_ka_max, radcoefs->ka_max);
   vx_ks_min = MMIN(vx_ks_min, radcoefs->ks_min);
   vx_ks_max = MMAX(vx_ks_max, radcoefs->ks_max);
   vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)] = vx_ka_min;
   vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)] = vx_ka_max;
   vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)] = vx_ks_min;
   vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] = vx_ks_max;
 }
 
 static res_T
 voxelize_tetra
   (struct rnatm* atm,
    const struct voxelize_args* args,
    struct partition* part, /* Partition to fill */
    const size_t cpnt)
 {
   const struct suvm_volume* volume = NULL;
   size_t i;
   ASSERT(atm && check_voxelize_args(args) == RES_OK);
   ASSERT(cpnt == RNATM_GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
 
   /* Fetch the component volumetric mesh */
   if(cpnt == RNATM_GAS) {
     volume = atm->gas.volume;
   } else {
     volume = darray_aerosol_cdata_get(&atm->aerosols)[cpnt].volume;
   }
 
   FOR_EACH(i, 0, darray_size_t_size_get(args->tetra_ids)) {
     struct suvm_primitive tetra;
     struct suvm_polyhedron poly;
     double poly_low[3];
     double poly_upp[3];
     float vx_low[3];
     float vx_upp[3];
     uint32_t ivx_low[3];
     uint32_t ivx_upp[3];
     uint32_t ivx[3];
     uint64_t mcode[3]; /* Cache of 3D morton code */
     size_t iitem;
     const size_t itetra = darray_size_t_cdata_get(args->tetra_ids)[i];
     enum suvm_intersection_type intersect;
 
     /* Recover the tetrahedron and setup its polyhedron */
     SUVM(volume_get_primitive(volume, itetra, &tetra));
     SUVM(primitive_setup_polyhedron(&tetra, &poly));
     ASSERT(poly.lower[0] <= args->part_upp[0]);
     ASSERT(poly.lower[1] <= args->part_upp[1]);
     ASSERT(poly.lower[2] <= args->part_upp[2]);
     ASSERT(poly.upper[0] >= args->part_low[0]);
     ASSERT(poly.upper[1] >= args->part_low[1]);
     ASSERT(poly.upper[2] >= args->part_low[2]);
 
     /* Clamp the AABB of the polyhedra to the partition bounds */
     poly_low[0] = MMAX(poly.lower[0], args->part_low[0]);
     poly_low[1] = MMAX(poly.lower[1], args->part_low[1]);
     poly_low[2] = MMAX(poly.lower[2], args->part_low[2]);
     poly_upp[0] = MMIN(poly.upper[0], args->part_upp[0]);
     poly_upp[1] = MMIN(poly.upper[1], args->part_upp[1]);
     poly_upp[2] = MMIN(poly.upper[2], args->part_upp[2]);
 
     /* Transform the AABB of the polyhedron in voxel space of the partition */
     ivx_low[0] = (uint32_t)((poly_low[0] - args->part_low[0]) / args->vxsz[0]);
     ivx_low[1] = (uint32_t)((poly_low[1] - args->part_low[1]) / args->vxsz[1]);
     ivx_low[2] = (uint32_t)((poly_low[2] - args->part_low[2]) / args->vxsz[2]);
     ivx_upp[0] = (uint32_t)ceil((poly_upp[0] - args->part_low[0]) / args->vxsz[0]);
     ivx_upp[1] = (uint32_t)ceil((poly_upp[1] - args->part_low[1]) / args->vxsz[1]);
     ivx_upp[2] = (uint32_t)ceil((poly_upp[2] - args->part_low[2]) / args->vxsz[2]);
     ASSERT(ivx_upp[0] <= args->part_def);
     ASSERT(ivx_upp[1] <= args->part_def);
     ASSERT(ivx_upp[2] <= args->part_def);
 
     /* Compute the range of the tetrahedron radiative coefficients */
     FOR_EACH(iitem, 0, args->batch_size) {
       const size_t iband = args->accel_structs[iitem].iband;
       const size_t iquad_pt = args->accel_structs[iitem].iquad_pt;
       struct radcoefs* radcoefs = args->per_item_radcoefs + iitem;
       setup_tetra_radcoefs(atm, &tetra, cpnt, iband, iquad_pt, radcoefs);
     }
 
     /* Iterate voxels intersected by the AABB of the polyedron */
     FOR_EACH(ivx[2], ivx_low[2], ivx_upp[2]) {
       vx_low[2] = (float)((double)ivx[2]*args->vxsz[2] + args->part_low[2]);
       vx_upp[2] = vx_low[2] + (float)args->vxsz[2];
       mcode[2] = morton3D_encode_u21(ivx[2]);
 
       FOR_EACH(ivx[1], ivx_low[1], ivx_upp[1]) {
         vx_low[1] = (float)((double)ivx[1]*args->vxsz[1] + args->part_low[1]);
         vx_upp[1] = vx_low[1] + (float)args->vxsz[1];
         mcode[1] = (morton3D_encode_u21(ivx[1]) << 1) | mcode[2];
 
         FOR_EACH(ivx[0], ivx_low[0], ivx_upp[0]) {
           vx_low[0] = (float)((double)ivx[0]*args->vxsz[0] + args->part_low[0]);
           vx_upp[0] = vx_low[0] + (float)args->vxsz[0];
           mcode[0] = (morton3D_encode_u21(ivx[0]) << 2) | mcode[1];
 
           intersect = suvm_polyhedron_intersect_aabb(&poly, vx_low, vx_upp);
           if(intersect == SUVM_INTERSECT_NONE) continue;
 
           /* Update the intersected voxel */
           FOR_EACH(iitem, 0, args->batch_size) {
             struct radcoefs* radcoefs = args->per_item_radcoefs + iitem;
             update_voxel(atm, radcoefs, part, mcode[0], iitem);
           }
         }
       }
     }
   }
 
   return RES_OK;
 }
 
 static res_T
 voxelize_partition(struct rnatm* atm, const struct voxelize_args* args)
 {
   STATIC_ASSERT((RNATM_GAS+1) == 0, Unexpected_constant_value);
 
   size_t naerosols = 0;
   size_t cpnt = 0;
   int part_is_empty = 1;
   res_T res = RES_OK;
   ASSERT(atm && check_voxelize_args(args) == RES_OK);
 
   partition_clear_voxels(args->part);
 
   /* Voxelize atmospheric components */
   naerosols = darray_aerosol_size_get(&atm->aerosols);
 
   /* CAUTION: in the following loop, it is assumed that the first component is
    * necessarily the gas which is therefore voxelized directly into the
    * partition. Aerosols are voxelized in a dummy partition before their
    * accumulation in the partition */
   cpnt = RNATM_GAS;
 
   do {
     struct suvm_volume* volume = NULL;
 
     /* Get the volumetric mesh of the component */
     if(cpnt == RNATM_GAS) {
       volume = atm->gas.volume;
     } else {
       volume = darray_aerosol_cdata_get(&atm->aerosols)[cpnt].volume;
     }
 
     /* Find the list of tetrahedra that overlap the partition */
     darray_size_t_clear(args->tetra_ids);
     res = suvm_volume_intersect_aabb
       (volume, args->part_low, args->part_upp, register_tetra, args->tetra_ids);
     if(res != RES_OK) goto error;
 
     /* The partition is not covered by any tetrahedron of the component */
     if(darray_size_t_size_get(args->tetra_ids) == 0) continue;
 
     /* The partition is covered by at least one tetrahedron and is therefore no
      * longer empty */
     part_is_empty = 0;
 
     /* Voxelize the gas directly into the partition */
     if(cpnt == RNATM_GAS) {
       res = voxelize_tetra(atm, args, args->part, cpnt);
       if(res != RES_OK) goto error;
 
     /* Voxelize each aerosol into the dummy partition before its accumulation */
     } else {
       partition_clear_voxels(args->part_dummy);
       res = voxelize_tetra(atm, args, args->part_dummy, cpnt);
       if(res != RES_OK) goto error;
       partition_accum(args->part, args->part_dummy);
     }
   } while(++cpnt < naerosols);
 
   /* The partition is not covered by any tetrahedron */
   if(part_is_empty) partition_empty(args->part);
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
 {
   int64_t i;
   int progress = 0;
   ATOMIC nparts_voxelized;
   ATOMIC res = RES_OK;
   ASSERT(atm && check_voxelize_batch_args(atm, args) == RES_OK);
 
   pool_reset(args->pool); /* Reset the next partition id to 0 */
 
   /* #partitions already voxelized */
   nparts_voxelized = (ATOMIC)args->nparts_voxelized;
 
   /* Iterates over the partitions of the grid according to their Morton order
    * and voxelizes the tetrahedrons that overlap them */
   omp_set_num_threads((int)atm->nthreads);
   #pragma omp parallel for schedule(static, 1/*chunk size*/)
   for(i = 0; i < (int64_t)args->nparts_adjusted; ++i) {
     struct voxelize_args voxelize_args = VOXELIZE_ARGS_NULL;
     struct partition* part_dummy = NULL;
     struct darray_size_t* tetra_ids = NULL;
     struct partition* part = NULL;
     struct radcoefs* per_item_radcoefs = NULL;
     double part_low[3];
     double part_upp[3];
     uint32_t part_ids[3];
     size_t n;
     const int ithread = omp_get_thread_num();
     int pcent;
     res_T res_local = RES_OK;
 
     if(ATOMIC_GET(&res) != RES_OK) continue;
 
     /* Recover the batch partitions */
     part = pool_next_partition(args->pool);
     if(!part) {
       ATOMIC_SET(&res, RES_UNKNOWN_ERR);
       continue;
     };
 
     /* Is the current partition out of bounds of the atmosphere grid */
     morton_xyz_decode_u21((uint64_t)partition_get_id(part), part_ids);
     if(part_ids[0] >= args->nparts[0]
     || part_ids[1] >= args->nparts[1]
     || part_ids[2] >= args->nparts[2]) {
       partition_free(part);
       continue;
     }
 
     /* Compute the AABB of the partition */
     part_low[0] = (double)(part_ids[0] * args->part_def) * args->vxsz[0];
     part_low[1] = (double)(part_ids[1] * args->part_def) * args->vxsz[1];
     part_low[2] = (double)(part_ids[2] * args->part_def) * args->vxsz[2];
     part_low[0] = part_low[0] + args->atm_low[0];
     part_low[1] = part_low[1] + args->atm_low[1];
     part_low[2] = part_low[2] + args->atm_low[2];
     part_upp[0] = part_low[0] + (double)args->part_def * args->vxsz[0];
     part_upp[1] = part_low[1] + (double)args->part_def * args->vxsz[1];
     part_upp[2] = part_low[2] + (double)args->part_def * args->vxsz[2];
 
     /* Retrieves the array where to store the indices of tetrahedra that
      * overlap the partition */
     tetra_ids = darray_size_t_list_data_get(args->per_thread_tetra_list) + ithread;
 
     /* Retrieve the spectral item data structure used to store radiative
      * coeficients */
     per_item_radcoefs = darray_radcoef_list_data_get
       (args->per_thread_radcoef_list)[ithread];
 
     /* Obtain the dummy partition needed to temporarily store aerosol
      * voxelization prior to accumulation in the current partition */
     part_dummy = darray_partition_data_get
       (args->per_thread_dummy_partition)[ithread];
 
     /* Voxelizes the partition and once done, commits */
     voxelize_args.tetra_ids = tetra_ids;
     voxelize_args.part_def = args->part_def;
     d3_set(voxelize_args.part_low, part_low);
     d3_set(voxelize_args.part_upp, part_upp);
     d3_set(voxelize_args.vxsz, args->vxsz);
     voxelize_args.part = part;
     voxelize_args.part_dummy = part_dummy;
     voxelize_args.per_item_radcoefs = per_item_radcoefs;
     voxelize_args.accel_structs = args->accel_structs;
     voxelize_args.batch_size = args->batch_size;
     res_local = voxelize_partition(atm, &voxelize_args);
     if(res_local == RES_OK) {
       partition_commit(part, args->batch_size);
     } else {
       partition_free(part);
       ATOMIC_SET(&res, res_local);
       continue;
     };
 
     /* Update progress bar */
     n = (size_t)ATOMIC_INCR(&nparts_voxelized);
     pcent = (int)((n * 100) / args->nparts_overall);
     #pragma omp critical
     if(pcent > progress) {
       progress = pcent;
       log_info(atm, VOXELIZE_MSG, pcent);
     }
   }
   if(res != RES_OK) goto error;
 
 exit:
   return (res_T)res;
 error:
   goto exit;
 }
 
 static res_T
 voxelize_atmosphere
   (struct rnatm* atm,
    struct pool* pool,
    struct build_sync* sync)
 {
   /* Batch of spectral items to process in a single voxelization */
   struct voxelize_batch_args batch_args = VOXELIZE_BATCH_ARGS_NULL;
   size_t batch_size = 0;
   size_t nbatches = 0;
   size_t ibatch = 0;
 
   /* Working data structures */
   struct darray_size_t_list per_thread_tetra_list;
   struct darray_radcoef_list per_thread_radcoef_list;
   struct darray_partition per_thread_dummy_partition;
 
   /* Voxelization parameters */
   double atm_low[3];
   double atm_upp[3];
   double vxsz[3];
   size_t nparts[3]; /* #partitions along the 3 axis */
   size_t nparts_adjusted; /* #partitions allowing their indexing by morton id */
   size_t nparts_overall; /* Total number of voxelized partitions */
   size_t part_def; /* Definition of a partition */
 
   /* Miscellaneous variables */
   struct accel_struct* accel_structs = NULL;
   size_t naccel_structs = 0;
   size_t i = 0;
   res_T res = RES_OK;
 
   ASSERT(atm && pool && sync);
 
   darray_size_t_list_init(atm->allocator, &per_thread_tetra_list);
   darray_radcoef_list_init(atm->allocator, &per_thread_radcoef_list);
   darray_partition_init(atm->allocator, &per_thread_dummy_partition);
 
   /* Allocate the per thread lists */
   res = darray_size_t_list_resize(&per_thread_tetra_list, atm->nthreads);
   if(res != RES_OK) goto error;
   res = darray_radcoef_list_resize(&per_thread_radcoef_list, atm->nthreads);
   if(res != RES_OK) goto error;
   res = darray_partition_resize(&per_thread_dummy_partition, atm->nthreads);
   if(res != RES_OK) goto error;
 
   /* Setup the per thread and per item working data structures */
   batch_size = pool_get_voxel_width(pool);
   FOR_EACH(i, 0, atm->nthreads) {
     struct radcoefs* per_item_k = NULL;
     per_item_k = MEM_CALLOC(atm->allocator, batch_size, sizeof(*per_item_k));
     if(!per_item_k) { res = RES_MEM_ERR; goto error; }
     darray_radcoef_list_data_get(&per_thread_radcoef_list)[i] = per_item_k;
   }
 
   /* Setup the per thread dummy partition */
   FOR_EACH(i, 0, atm->nthreads) {
     struct partition* part = pool_dummy_partition(pool);
     darray_partition_data_get(&per_thread_dummy_partition)[i] = part;
   }
 
   /* Recover the AABB atmosphere and compute the size of a voxel */
   SUVM(volume_get_aabb(atm->gas.volume, atm_low, atm_upp));
   vxsz[0] = (atm_upp[0] - atm_low[0]) / (double)atm->grid_definition[0];
   vxsz[1] = (atm_upp[1] - atm_low[1]) / (double)atm->grid_definition[1];
   vxsz[2] = (atm_upp[2] - atm_low[2]) / (double)atm->grid_definition[2];
 
   /* Number of partitions required to cover the entire atmosphere grid */
   part_def = pool_get_partition_definition(pool);
   nparts[0] = (atm->grid_definition[0] + (part_def-1)/*ceil*/) / part_def;
   nparts[1] = (atm->grid_definition[1] + (part_def-1)/*ceil*/) / part_def;
   nparts[2] = (atm->grid_definition[2] + (part_def-1)/*ceil*/) / part_def;
 
   /* Adjust the #partitions allowing their indexing by their morton code */
   nparts_adjusted = MMAX(nparts[0], MMAX(nparts[1], nparts[2]));
   nparts_adjusted = round_up_pow2(nparts_adjusted);
   nparts_adjusted = nparts_adjusted * nparts_adjusted * nparts_adjusted;
 
   /* Calculate the number of batches required to build the accelerating
    * structures for the spectral bands submitted; currently we are building one
    * octree per quadrature point of each band. And each octree requires a
    * complete voxelization of the atmospheric meshes which takes time. To
    * amortize this cost without prohibitive increase in memory space, one fills
    * up to N octrees from a single voxelization of the mesh. The number of
    * batches corresponds to the total number of voxelization required */
   accel_structs = darray_accel_struct_data_get(&atm->accel_structs);
   naccel_structs = darray_accel_struct_size_get(&atm->accel_structs);
   nbatches = (naccel_structs + (batch_size - 1)/*ceil*/) / batch_size;
 
   /* Calculate the total number of partitions to voxelize. Note that the same
    * partition can be voxelized several times, once per batch */
   nparts_overall = nparts[0] * nparts[1] * nparts[2] * nbatches;
 
   /* Print the size of a voxel */
   log_info(atm, "voxel size = {%g, %g, %g}\n", SPLIT3(vxsz));
 
   /* Setup regular batch arguments */
   batch_args.per_thread_tetra_list = &per_thread_tetra_list;
   batch_args.per_thread_radcoef_list = &per_thread_radcoef_list;
   batch_args.per_thread_dummy_partition = &per_thread_dummy_partition;
   batch_args.pool = pool;
   batch_args.part_def = part_def;
   batch_args.nparts[0] = nparts[0];
   batch_args.nparts[1] = nparts[1];
   batch_args.nparts[2] = nparts[2];
   batch_args.nparts_adjusted = nparts_adjusted;
   batch_args.nparts_overall = nparts_overall;
   batch_args.nparts_voxelized = 0;
   d3_set(batch_args.atm_low, atm_low);
   d3_set(batch_args.atm_upp, atm_upp);
   d3_set(batch_args.vxsz, vxsz);
 
   /* Print voxelization status */
   log_info(atm, VOXELIZE_MSG, 0);
 
   /* Voxelize the batches */
   FOR_EACH(ibatch, 0, nbatches) {
     size_t item_range[2];
     item_range[0] = ibatch * batch_size;
     item_range[1] = MMIN(item_range[0] + batch_size, naccel_structs);
 
     batch_args.accel_structs = accel_structs + item_range[0];
     batch_args.batch_size = item_range[1] - item_range[0];
 
     /* Wait for the building thread to finish consuming the previous batch */
     mutex_lock(sync->mutex);
     if(sync->ibatch != ibatch) {
       ASSERT(sync->ibatch == ibatch - 1);
       cond_wait(sync->cond, sync->mutex);
       /* An error occured in the building thread */
       if(sync->ibatch != ibatch) res = RES_BAD_ARG;
     }
     mutex_unlock(sync->mutex);
     if(res != RES_OK) goto error;
 
     /* Generate the voxels of the current batch */
     res = voxelize_batch(atm, &batch_args);
     if(res != RES_OK) goto error;
 
     batch_args.nparts_voxelized += nparts[0]*nparts[1]*nparts[2];
   }
 
   /* Print final status message */
   log_info(atm, VOXELIZE_MSG"\n", 100);
 
 exit:
   FOR_EACH(i, 0, darray_radcoef_list_size_get(&per_thread_radcoef_list)) {
     struct radcoefs* per_item_k = NULL;
     per_item_k = darray_radcoef_list_data_get(&per_thread_radcoef_list)[i];
     MEM_RM(atm->allocator,per_item_k);
   }
   FOR_EACH(i, 0, darray_partition_size_get(&per_thread_dummy_partition)) {
     struct partition* part = NULL;
     part = darray_partition_data_get(&per_thread_dummy_partition)[i];
     partition_free(part);
   }
   darray_size_t_list_release(&per_thread_tetra_list);
   darray_radcoef_list_release(&per_thread_radcoef_list);
   darray_partition_release(&per_thread_dummy_partition);
   return res;
 error:
   goto exit;
 }
 
 static void
 vx_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
 {
   struct build_octree_context* ctx = context;
   const float* vx = NULL;
   uint64_t ivx, ipart;
   int log2_part_def;
   ASSERT(xyz && dst && ctx);
   (void)xyz;
 
   /* Retrieve the partition's morton id and the voxel's morton id in the
    * partition from the morton code of the voxel in the whole octree. Note
    * that, like voxels, partitions are sorted in morton order. Therefore, the
    * partition's morton id is encoded in the most significant bits of the
    * voxel's morton code, while the voxel's morton id in the partition is
    * stored in its least significant bits */
   log2_part_def = log2i((int)pool_get_partition_definition(ctx->pool));
   ipart = (mcode >> (log2_part_def*3));
   ivx = (mcode & (BIT_U64(log2_part_def*3)-1));
 
   /* Recover the partition storing the voxel */
   if(ctx->part == NULL || partition_get_id(ctx->part) != ipart) {
     if(ctx->part) partition_free(ctx->part);
     ctx->part = pool_fetch_partition(ctx->pool, ipart);
 
     if(ctx->part == NULL) { /* An error occurs */
       memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
       return;
     }
   }
 
   vx = partition_cget_voxel(ctx->part, ivx, ctx->iitem);
   memcpy(dst, vx, NFLOATS_PER_VOXEL * sizeof(float));
 }
 
 static void
 vx_merge(void* dst, const void* vxs[], const size_t nvxs, void* ctx)
 {
   float ka_min = FLT_MAX, ka_max = -FLT_MAX;
   float ks_min = FLT_MAX, ks_max = -FLT_MAX;
   float* merged_vx = dst;
   size_t ivx;
   ASSERT(dst && vxs && nvxs);
   (void)ctx;
 
   FOR_EACH(ivx, 0, nvxs) {
     const float* vx = vxs[ivx];
     ka_min = MMIN(ka_min, vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)]);
     ka_max = MMAX(ka_max, vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)]);
     ks_min = MMIN(ks_min, vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)]);
     ks_max = MMAX(ks_max, vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)]);
   }
 
   merged_vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)] = ka_min;
   merged_vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)] = ka_max;
   merged_vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)] = ks_min;
   merged_vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] = ks_max;
 }
 
 static int
 vx_challenge_merge
   (const struct svx_voxel vxs[],
    const size_t nvxs,
    void* context)
 {
   const struct build_octree_context* ctx = context;
   double low[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
   double upp[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
   double tau;
   float kext_min = FLT_MAX, kext_max = -FLT_MAX;
   double sz_max;
   size_t ivx;
   ASSERT(vxs && nvxs && context);
 
   /* Compute the range of the extinction coefficients of the submitted voxels */
   FOR_EACH(ivx, 0, nvxs) {
     const float* vx = vxs[ivx].data;
     const float ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
     const float ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
     const float ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
     const float ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
     if(ka_min > ka_max) { ASSERT(ks_min > ks_max); continue; } /* Empty voxel */
     kext_min = MMIN(kext_min, ka_min + ks_min);
     kext_max = MMAX(kext_max, ka_max + ks_max);
   }
 
   /* Always merge empty voxels */
   if(kext_min > kext_max)  return 1;
 
   /* Compute the AABB of the submitted voxels */
   FOR_EACH(ivx, 0, nvxs) {
     low[0] = MMIN(vxs[ivx].lower[0], low[0]);
     low[1] = MMIN(vxs[ivx].lower[1], low[1]);
     low[2] = MMIN(vxs[ivx].lower[2], low[2]);
     upp[0] = MMAX(vxs[ivx].upper[0], upp[0]);
     upp[1] = MMAX(vxs[ivx].upper[1], upp[1]);
     upp[2] = MMAX(vxs[ivx].upper[2], upp[2]);
   }
 
   /* Compute the size of the largest dimension of the AABB */
   sz_max = upp[0] - low[0];
   sz_max = MMAX(upp[1] - low[1], sz_max);
   sz_max = MMAX(upp[2] - low[2], sz_max);
 
   /* Check if voxels can be merged by comparing the optical thickness along the
    * maximum size of the merged AABB against a user-defined threshold */
   tau = (kext_max - kext_min) * sz_max;
   return tau < ctx->tau_threshold;
 }
 
 static res_T
 init_octrees_storage(struct rnatm* atm, const struct rnatm_create_args* args)
 {
   struct octrees_storage_desc challenge = OCTREES_STORAGE_DESC_NULL;
   struct octrees_storage_desc reference = OCTREES_STORAGE_DESC_NULL;
   res_T res = RES_OK;
   ASSERT(atm && args);
 
   if(!args->octrees_storage) goto exit;
 
   /* Register the storage file */
   atm->octrees_storage = args->octrees_storage;
 
   res = octrees_storage_desc_init(atm, &challenge);
   if(res != RES_OK) goto error;
 
   /* Save the computed descriptor in the storage */
   if(!args->load_octrees_from_storage) {
     res = write_octrees_storage_desc(atm, &challenge, atm->octrees_storage);
     if(res != RES_OK) goto error;
 
   /* The octrees must be loaded from storage. Verify that the data and settings
    * used by the saved octrees are compatibles with the current parameters */
   } else {
 
     res = octrees_storage_desc_init_from_stream
       (atm, &reference, args->octrees_storage);
     if(res != RES_OK) goto error;
 
     res = check_octrees_storage_compatibility(&reference, &challenge);
     if(res != RES_OK) {
       log_err(atm,
         "unable to load octrees from the supplied storage. The saved "
         "octrees are built from different data\n");
       res = RES_BAD_ARG;
       goto error;
     }
   }
 
 exit:
   octrees_storage_desc_release(&challenge);
   octrees_storage_desc_release(&reference);
   return res;
 error:
   goto exit;
 }
 
 static res_T
 create_storage_toc(struct rnatm* atm, fpos_t* toc)
 {
   int err = 0;
   res_T res = RES_OK;
   ASSERT(atm && toc);
 
   if(!atm->octrees_storage) goto exit;
 
   err = fgetpos(atm->octrees_storage, toc);
   if(err != 0) {
     log_err(atm, "error retrieving toc position of octrees storage -- %s\n",
       strerror(errno));
     res = RES_IO_ERR;
     goto error;
   }
 
   res = reserve_octrees_storage_toc(atm, atm->octrees_storage);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   log_err(atm, "error reserving octree storage header -- %s\n",
     strerror(errno));
   goto exit;
 }
 
 static res_T
 offload_octrees
   (struct rnatm* atm,
    const size_t ioctrees[2]) /* Octrees to write. Limits are inclusive */
 {
   size_t ioctree = 0;
   res_T res = RES_OK;
   ASSERT(atm && ioctrees && ioctrees[0] <= ioctrees[1]);
   ASSERT(ioctrees[1] < darray_accel_struct_size_get(&atm->accel_structs));
 
   /* No storage: nothing to do */
   if(!atm->octrees_storage) goto exit;
 
   res = store_octrees(atm, ioctrees, atm->octrees_storage);
   if(res != RES_OK) goto error;
 
   FOR_EACH(ioctree, ioctrees[0], ioctrees[1]+1) {
     unload_octree(atm, ioctree); /* Free the octree memory */
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 finalize_storage(struct rnatm* atm, const fpos_t* toc)
 {
   int err = 0;
   res_T res = RES_OK;
   ASSERT(atm && toc);
 
   /* No storage: nothing to do */
   if(!atm->octrees_storage) goto exit;
 
   err = fsetpos(atm->octrees_storage, toc);
   if(err != 0) {
     log_err(atm, "error positioning on octree storage toc -- %s\n",
       strerror(errno));
     res = RES_IO_ERR;
     goto error;
   }
 
   /* Update the table of content */
   res = write_octrees_storage_toc(atm, atm->octrees_storage);
   if(res != RES_OK) goto exit;
 
   CHK(fflush(atm->octrees_storage) == 0);
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 build_octrees
   (struct rnatm* atm,
    const struct rnatm_create_args* args,
    struct pool* pool,
    struct build_sync* sync)
 {
   fpos_t storage_toc;
   struct accel_struct* accel_structs = NULL;
   double low[3], upp[3];
   size_t def[3];
   size_t istruct;
   size_t naccel_structs;
   size_t voxel_width;
   ATOMIC res = RES_OK;
   ASSERT(atm && args && pool);
 
   /* Recover the AABB from the atmosphere. Note that we have already made sure
    * when setting up the meshes of the gas and aerosols that the aerosols are
    * included in the gas. Therefore, the AABB of the gas is the same as the
    * AABB of the atmosphere */
   SUVM(volume_get_aabb(atm->gas.volume, low, upp));
 
   /* Retrieve grid definition */
   def[0] = (size_t)atm->grid_definition[0];
   def[1] = (size_t)atm->grid_definition[1];
   def[2] = (size_t)atm->grid_definition[2];
 
   accel_structs = darray_accel_struct_data_get(&atm->accel_structs);
   naccel_structs = darray_accel_struct_size_get(&atm->accel_structs);
   voxel_width = pool_get_voxel_width(pool);
 
   /* Setup the disk storage of the octrees */
   res = create_storage_toc(atm, &storage_toc);
   if(res != RES_OK) goto error;
 
   /* Build the octrees. Each thread consumes an element of the voxels generated
    * by the voxelization thread, each element corresponding to the voxel of an
    * octree to be constructed. By fixing the number of threads to the width of
    * the voxel, we therefore build `voxel_width' octrees in parallel from a
    * single voxelization of the atmospheric meshes */
   for(istruct = 0; istruct < naccel_structs; istruct += voxel_width) {
     const size_t batch_size = MMIN(voxel_width, naccel_structs - istruct);
     size_t ioctrees[2];
     omp_set_num_threads((int)batch_size);
 
     /* Note that we are using a parallel block rather than a parallel loop in
      * order to add an implicit barrier after a batch has been fully consumed.
      * This is necessary to prevent a thread from consuming voxels from the
      * previous batch */
     #pragma omp parallel
     {
       struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
       struct svx_voxel_desc vx_desc = SVX_VOXEL_DESC_NULL;
       struct svx_tree* octree = NULL;
       const int ithread = omp_get_thread_num();
       const size_t istruct_curr = (size_t)ithread + istruct;
       res_T res_local = RES_OK;
 
       /* Setup the build context */
       ctx.pool = pool;
       ctx.part = NULL;
       ctx.iitem = (size_t)ithread;
       ctx.tau_threshold = args->optical_thickness;
 
       /* Setup the voxel descriptor */
       vx_desc.get = vx_get;
       vx_desc.merge = vx_merge;
       vx_desc.challenge_merge = vx_challenge_merge;
       vx_desc.context = &ctx;
       vx_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
 
       res_local = svx_octree_create(atm->svx, low, upp, def, &vx_desc, &octree);
       if(ctx.part) partition_free(ctx.part);
       if(res_local != RES_OK) {
         ATOMIC_SET(&res, res_local);
       } else { /* Register the built octree */
         accel_structs[istruct_curr].octree = octree;
       }
     }
     if(res != RES_OK) goto error;
 
     /* Signal the voxelization thread to generate the next batch */
     mutex_lock(sync->mutex);
     sync->ibatch += 1;
     mutex_unlock(sync->mutex);
     cond_signal(sync->cond);
 
     /* Offload builded octrees */
     ioctrees[0] = istruct;
     ioctrees[1] = istruct + batch_size - 1;
     res = offload_octrees(atm, ioctrees);
     if(res != RES_OK) goto error;
   }
 
   /* Finalize the file where octrees are offloaded */
   res = finalize_storage(atm, &storage_toc);
   if(res != RES_OK) goto error;
 
 exit:
   return (res_T)res;
 error:
   /* Signal to the voxelization thread that there is no need to wait for the
    * build thread */
   cond_signal(sync->cond);
   darray_accel_struct_clear(&atm->accel_structs);
   goto exit;
 }
 
 static res_T
 create_pool(struct rnatm* atm, struct pool** out_pool, const size_t nbands)
 {
   /* Empirical constant that defines the number of non empty partitions to
    * pre-allocate per thread to avoid contension between the thread building the
    * octrees from the partitions and the threads that fill these partitions */
   const size_t NPARTITIONS_PER_THREAD = 64;
 
   /* Number of dummy partitions per thread. These partitions are used to
    * temporarily store the voxelization of aerosols that are then accumulated
    * in the voxelized partition. Actually, one such partition per voxelization
    * thread is required */
   size_t NPARTITIONS_PER_THREAD_DUMMY = 1;
 
   /* Empirical constant that defines the total number of partitions to
    * pre-allocate per thread. This constant is necessarily greater than or
    * equal to (NPARTITIONS_PER_THREAD + NPARTITIONS_PER_THREAD_DUMMY), the
    * difference representing the number of completely empty partitions. Such
    * partitions help reduce thread contention with a sparse grid without
    * significantly increasing memory usage */
   const size_t OVERALL_NPARTITIONS_PER_THREAD =
     NPARTITIONS_PER_THREAD * 64
   + NPARTITIONS_PER_THREAD_DUMMY;
 
   /* Number of items stored per voxel data. A width greater than 1 makes it
    * possible to store by partition the radiative coefficients of several
    * spectral data. The goal is to voxelize the volumetric mesh once for N
    * spectral data */
   const size_t VOXEL_WIDTH = MMIN(4, nbands);
 
   /* Definition of a partition on the 3 axis */
   const size_t PARTITION_DEFINITION = 8;
 
   struct pool_create_args pool_args = POOL_CREATE_ARGS_DEFAULT;
   struct pool* pool = NULL;
   res_T res = RES_OK;
   ASSERT(atm && out_pool);
 
   /* Create the partition pool */
   pool_args.npartitions = atm->nthreads * OVERALL_NPARTITIONS_PER_THREAD;
   pool_args.npreallocated_partitions =  atm->nthreads * NPARTITIONS_PER_THREAD;
   pool_args.partition_definition = PARTITION_DEFINITION;
   pool_args.voxel_width = VOXEL_WIDTH;
   pool_args.allocator = atm->allocator;
   res = pool_create(&pool_args, &pool);
   if(res != RES_OK) goto error;
 
 exit:
   *out_pool = pool;
   return res;
 error:
   log_err(atm, "failed to create the partition pool -- %s\n",
     res_to_cstr(res));
   if(pool) { pool_ref_put(pool); pool = NULL; }
   goto exit;
 }
 
 static res_T
 create_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
 {
   struct build_sync sync = BUILD_SYNC_NULL;
   struct pool* pool = NULL;
   size_t nbands = 0;
 
   ATOMIC res = RES_OK;
   ASSERT(atm);
 
   nbands = darray_band_size_get(&atm->bands);
 
   res = create_pool(atm, &pool, nbands);
   if(res != RES_OK) goto error;
   res = build_sync_init(&sync);
   if(res != RES_OK) goto error;
 
   log_info(atm,
     "partitionning of radiative properties "
     "(grid definition = %ux%ux%u; #octrees = %lu)\n",
     SPLIT3(atm->grid_definition),
     (unsigned long)darray_accel_struct_size_get(&atm->accel_structs));
 
   /* Enable nested threads to allow multi-threading when voxelizing tetrahedra
    * and building octrees */
   omp_set_nested(1);
 
   #pragma omp parallel sections num_threads(2)
   {
     #pragma omp section
     {
       const res_T res_local = voxelize_atmosphere(atm, pool, &sync);
       if(res_local != RES_OK) {
         log_err(atm, "error when voxelizing the atmosphere -- %s\n",
           res_to_cstr((res_T)res_local));
         pool_invalidate(pool);
         ATOMIC_SET(&res, res_local);
       }
     }
 
     #pragma omp section
     {
       const res_T res_local = build_octrees(atm, args, pool, &sync);
       if(res_local != RES_OK) {
         log_err(atm, "error building octrees -- %s\n",
           res_to_cstr((res_T)res_local));
         pool_invalidate(pool);
         ATOMIC_SET(&res, res_local);
       }
     }
   }
   if(res != RES_OK) goto error;
 
 exit:
   build_sync_release(&sync);
   if(pool) pool_ref_put(pool);
   return (res_T)res;
 error:
   goto exit;
 }
 
 static INLINE res_T
 load_octrees(struct rnatm* atm)
 {
   res_T res = RES_OK;
   ASSERT(atm && atm->octrees_storage);
 
   /* Only load the storage table of content. The octrees will be loaded on
    * demand */
   res = read_octrees_storage_toc(atm, atm->octrees_storage);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /*******************************************************************************
  * Exported functions
  ******************************************************************************/
 res_T
 rnatm_trace_ray
   (struct rnatm* atm,
    const struct rnatm_trace_ray_args* args,
    struct svx_hit* hit)
 {
   const struct accel_struct* accel_struct = NULL;
   size_t i;
   size_t iaccel;
   res_T res = RES_OK;
 
   if(!atm || !hit) { res = RES_BAD_ARG; goto error; }
   res = check_trace_ray_args(atm, args);
   if(res != RES_OK) goto error;
 
   /* Start calculating the id of the acceleration structure to consider */
   i = args->iband - darray_band_cdata_get(&atm->bands)[0].index;
   ASSERT(i < darray_band_size_get(&atm->bands));
   ASSERT(args->iband == darray_band_cdata_get(&atm->bands)[i].index);
 
   /* Invalid quadrature point index */
   if(args->iquad >= darray_band_cdata_get(&atm->bands)[i].nquad_pts)
     return RES_BAD_ARG;
 
   /*  Find the id of the acceleration structure */
   iaccel = darray_band_cdata_get(&atm->bands)[i].offset_accel_struct + args->iquad;
 
   /* Retrieve the acceleration structure */
   ASSERT(i < darray_accel_struct_size_get(&atm->accel_structs));
   accel_struct = darray_accel_struct_cdata_get(&atm->accel_structs) + iaccel;
   ASSERT(args->iband == accel_struct->iband);
   ASSERT(args->iquad == accel_struct->iquad_pt);
 
   res = make_sure_octree_is_loaded(atm, iaccel);
   if(res != RES_OK) goto error;
 
   *hit = SVX_HIT_NULL;
 
   res = svx_tree_trace_ray
     (accel_struct->octree,
      args->ray_org,
      args->ray_dir,
      args->ray_range,
      args->challenge,
      args->filter,
      args->context,
      hit);
   if(res != RES_OK) {
     log_err(atm,
       "%s: error tracing ray "
       "(origin = %g, %g, %g; direction = %g, %g, %g; range = %g, %g) -- %s\n",
       FUNC_NAME, SPLIT3(args->ray_org), SPLIT3(args->ray_dir),
       SPLIT2(args->ray_range), res_to_cstr(res));
     goto error;
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
 res_T
 setup_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
 {
   char buf[128];
   struct time t0, t1;
   size_t sz;
   res_T res = RES_OK;
   ASSERT(atm && args);
 
   /* Start time recording */
   time_current(&t0);
 
   /* Create the Star-VoXel device */
   res = svx_device_create
     (atm->logger, &atm->svx_allocator, atm->verbose, &atm->svx);
   if(res != RES_OK) goto error;
 
   /* Setup miscellaneous parameters */
   atm->optical_thickness = args->optical_thickness;
 
   res = compute_grid_definition(atm, args);
   if(res != RES_OK) goto error;
 
   /* Initialize storage *after* calculating the grid definition since it is
    * saved in the storage file */
   res = init_octrees_storage(atm, args);
   if(res != RES_OK) goto error;
 
   res = setup_accel_structs(atm);
   if(res != RES_OK) goto error;
 
   if(args->load_octrees_from_storage) {
     res = load_octrees(atm);
     if(res != RES_OK) goto error;
   } else {
     res = create_octrees(atm, args);
     if(res != RES_OK) goto error;
   }
 
   /* Log elapsed time */
   time_sub(&t0, time_current(&t1), &t0);
   time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
   log_info(atm, "acceleration structures built in %s\n", buf);
 
   /* Log memory space used by Star-VX */
   sz = MEM_ALLOCATED_SIZE(&atm->svx_allocator);
   size_to_cstr(sz, SIZE_ALL, NULL, buf, sizeof(buf));
   log_info(atm, "Star-VoXel memory space: %s\n", buf);
 
 exit:
   return res;
 error:
   goto exit;
 }