stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

sdis_solve_camera.c (22497B)


      1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "sdis.h"
     17 #include "sdis_c.h"
     18 #include "sdis_camera.h"
     19 #include "sdis_device_c.h"
     20 #include "sdis_estimator_buffer_c.h"
     21 #include "sdis_log.h"
     22 #include "sdis_medium_c.h"
     23 #include "sdis_realisation.h"
     24 #include "sdis_scene_c.h"
     25 #include "sdis_tile.h"
     26 #ifdef SDIS_ENABLE_MPI
     27   #include "sdis_mpi.h"
     28 #endif
     29 
     30 #include <rsys/clock_time.h>
     31 #include <rsys/cstr.h>
     32 #include <rsys/list.h>
     33 #include <rsys/morton.h>
     34 
     35 #include <omp.h>
     36 
     37 /*******************************************************************************
     38  * Helper function
     39  ******************************************************************************/
     40 static res_T
     41 check_solve_camera_args(const struct sdis_solve_camera_args* args)
     42 {
     43   if(!args) return RES_BAD_ARG;
     44   if(!args->cam) return RES_BAD_ARG;
     45 
     46   /* Check the image resolution */
     47   if(!args->image_definition[0] || !args->image_definition[1]) {
     48     return RES_BAD_ARG;
     49   }
     50 
     51   /* Check the number of samples per pixel */
     52   if(!args->spp) {
     53     return RES_BAD_ARG;
     54   }
     55 
     56   /* Check the time range */
     57   if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
     58     return RES_BAD_ARG;
     59   }
     60   if(args->time_range[1] > DBL_MAX
     61   && args->time_range[0] != args->time_range[1]) {
     62     return RES_BAD_ARG;
     63   }
     64 
     65   /* Check the picard order */
     66   if(args->picard_order < 1) {
     67     return RES_BAD_ARG;
     68   }
     69 
     70   /* Check RNG type */
     71   if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) {
     72     return RES_BAD_ARG;
     73   }
     74 
     75   /* Check the diffusion algorithm */
     76   if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
     77     return RES_BAD_ARG;
     78   }
     79 
     80   return RES_OK;
     81 }
     82 
     83 static res_T
     84 solve_pixel
     85   (struct sdis_scene* scn,
     86    struct ssp_rng* rng,
     87    const unsigned enc_id,
     88    const struct sdis_camera* cam,
     89    const double time_range[2], /* Observation time */
     90    const size_t ipix[2], /* Pixel coordinate in the image space */
     91    const size_t nrealisations,
     92    const int register_paths, /* Combination of enum sdis_heat_path_flag */
     93    const double pix_sz[2], /* Pixel size in the normalized image plane */
     94    const size_t picard_order,
     95    const enum sdis_diffusion_algorithm diff_algo,
     96    struct sdis_estimator* estimator,
     97    struct pixel* pixel)
     98 {
     99   struct sdis_heat_path* pheat_path = NULL;
    100   size_t irealisation;
    101   res_T res = RES_OK;
    102   ASSERT(scn && rng && cam && ipix && nrealisations);
    103   ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
    104   ASSERT(pixel && time_range);
    105 
    106   pixel->acc_temp = ACCUM_NULL;
    107   pixel->acc_time = ACCUM_NULL;
    108 
    109   FOR_EACH(irealisation, 0, nrealisations) {
    110     struct ray_realisation_args realis_args = RAY_REALISATION_ARGS_NULL;
    111     struct time t0, t1;
    112     double samp[2]; /* Pixel sample */
    113     double ray_pos[3];
    114     double ray_dir[3];
    115     double w = 0;
    116     struct sdis_heat_path heat_path;
    117     double time;
    118     res_T res_simul = RES_OK;
    119 
    120     /* Begin time registration */
    121     time_current(&t0);
    122 
    123     time = sample_time(rng, time_range);
    124     if(register_paths) {
    125       heat_path_init(scn->dev->allocator, &heat_path);
    126       pheat_path = &heat_path;
    127     }
    128 
    129     /* Generate a sample into the pixel to estimate */
    130     samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
    131     samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
    132 
    133     /* Generate a ray starting from the camera position and passing through
    134      * pixel sample */
    135     camera_ray(cam, samp, ray_pos, ray_dir);
    136 
    137     /* Launch the realisation */
    138     realis_args.rng = rng;
    139     realis_args.enc_id = enc_id;
    140     realis_args.time = time;
    141     realis_args.picard_order = picard_order;
    142     realis_args.heat_path = pheat_path;
    143     realis_args.irealisation = (size_t)irealisation;
    144     realis_args.diff_algo = diff_algo;
    145     d3_set(realis_args.position, ray_pos);
    146     d3_set(realis_args.direction, ray_dir);
    147     res_simul = ray_realisation_3d(scn, &realis_args, &w);
    148 
    149     /* Handle fatal error */
    150     if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
    151       res = res_simul;
    152       goto error;
    153     }
    154 
    155     if(pheat_path) {
    156       pheat_path->status = res_simul == RES_OK
    157         ? SDIS_HEAT_PATH_SUCCESS
    158         : SDIS_HEAT_PATH_FAILURE;
    159 
    160       /* Check if the path must be saved regarding the register_paths mask */
    161       if(!(register_paths & (int)pheat_path->status)) {
    162         heat_path_release(pheat_path);
    163         pheat_path = NULL;
    164       } else { /* Register the sampled path */
    165         ASSERT(estimator);
    166         res = estimator_add_and_release_heat_path(estimator, pheat_path);
    167         if(res != RES_OK) goto error;
    168         pheat_path = NULL;
    169       }
    170     }
    171 
    172     /* Stop time registration */
    173     time_sub(&t0, time_current(&t1), &t0);
    174 
    175     if(res_simul == RES_OK) {
    176       /* Update pixel accumulators */
    177       const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
    178       pixel->acc_temp.sum += w;
    179       pixel->acc_temp.sum2 += w*w;
    180       pixel->acc_temp.count += 1;
    181       pixel->acc_time.sum += usec;
    182       pixel->acc_time.sum2 += usec*usec;
    183       pixel->acc_time.count += 1;
    184     }
    185   }
    186 
    187 exit:
    188   if(pheat_path) heat_path_release(pheat_path);
    189   return res;
    190 error:
    191   goto exit;
    192 }
    193 
    194 static res_T
    195 solve_tile
    196   (struct sdis_scene* scn,
    197    struct ssp_rng* rng,
    198    const unsigned enc_id,
    199    const struct sdis_camera* cam,
    200    const double time_range[2],
    201    const size_t tile_org[2], /* Origin of the tile in pixel space */
    202    const size_t tile_size[2], /* #pixels in the tile in X and Y */
    203    const size_t spp, /* #samples per pixel */
    204    const int register_paths, /* Combination of enum sdis_heat_path_flag */
    205    const double pix_sz[2], /* Pixel size in the normalized image plane */
    206    const size_t picard_order,
    207    const enum sdis_diffusion_algorithm diff_algo,
    208    struct sdis_estimator_buffer* buf,
    209    struct tile* tile)
    210 {
    211   size_t mcode; /* Morton code of the tile pixel */
    212   size_t npixels;
    213   res_T res = RES_OK;
    214   ASSERT(scn && rng && cam && spp);
    215   ASSERT(tile_size && tile_size[0] && tile_size[1]);
    216   ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range);
    217 
    218   /* Adjust the #pixels to process them wrt a morton order */
    219   npixels = round_up_pow2(MMAX(tile_size[0], tile_size[1]));
    220   npixels *= npixels;
    221 
    222   FOR_EACH(mcode, 0, npixels) {
    223     struct pixel* pixel = NULL;
    224     struct sdis_estimator* estimator = NULL;
    225     uint16_t ipix_tile[2];
    226     size_t ipix_image[2];
    227 
    228     ipix_tile[0] = morton2D_decode_u16((uint32_t)(mcode>>0));
    229     if(ipix_tile[0] >= tile_size[0]) continue;
    230     ipix_tile[1] = morton2D_decode_u16((uint32_t)(mcode>>1));
    231     if(ipix_tile[1] >= tile_size[1]) continue;
    232 
    233     pixel = tile_at(tile, ipix_tile[0], ipix_tile[1]);
    234 
    235     /* Compute the pixel coordinates in image space */
    236     ipix_image[0] = ipix_tile[0] + tile_org[0];
    237     ipix_image[1] = ipix_tile[1] + tile_org[1];
    238 
    239     if(register_paths != SDIS_HEAT_PATH_NONE) {
    240       ASSERT(buf);
    241       estimator = estimator_buffer_grab(buf, ipix_image[0], ipix_image[1]);
    242     }
    243     res = solve_pixel
    244       (scn, rng, enc_id, cam, time_range, ipix_image, spp, register_paths,
    245        pix_sz, picard_order, diff_algo, estimator, pixel);
    246     if(res != RES_OK) goto error;
    247   }
    248 
    249 exit:
    250   return res;
    251 error:
    252   goto exit;
    253 }
    254 
    255 static INLINE res_T
    256 setup_estimator_from_pixel
    257   (struct sdis_estimator* estimator,
    258    const size_t spp, /* #samples per pixel */
    259    struct pixel* pixel)
    260 {
    261   ASSERT(estimator && spp && pixel);
    262   ASSERT(pixel->acc_temp.count == pixel->acc_time.count);
    263   estimator_setup_realisations_count
    264     (estimator, spp, pixel->acc_temp.count);
    265   estimator_setup_temperature
    266     (estimator, pixel->acc_temp.sum, pixel->acc_temp.sum2);
    267   estimator_setup_realisation_time
    268     (estimator, pixel->acc_time.sum, pixel->acc_time.sum2);
    269   return RES_OK;
    270 }
    271 
    272 static res_T
    273 write_tile
    274   (struct sdis_estimator_buffer* buf,
    275    const size_t spp, /* #samples per pixel */
    276    struct tile* tile)
    277 {
    278   res_T res = RES_OK;
    279   size_t tile_org[2];
    280   size_t buf_sz[2];
    281   size_t tile_sz[2];
    282   uint16_t x, y;
    283   ASSERT(buf && spp && tile);
    284 
    285   SDIS(estimator_buffer_get_definition(buf, buf_sz));
    286 
    287   tile_org[0] = (size_t)(tile->data.x * TILE_SIZE);
    288   tile_org[1] = (size_t)(tile->data.y * TILE_SIZE);
    289   tile_sz[0] = MMIN(TILE_SIZE, buf_sz[0] - tile_org[0]);
    290   tile_sz[1] = MMIN(TILE_SIZE, buf_sz[1] - tile_org[1]);
    291 
    292   FOR_EACH(y, 0, tile_sz[1]) {
    293     const size_t pix_y = tile_org[1] + y;
    294     FOR_EACH(x, 0, tile_sz[0]) {
    295       const size_t pix_x = tile_org[0] + x;
    296       struct sdis_estimator* estimator = NULL;
    297       struct pixel* pixel = NULL;
    298 
    299       estimator = estimator_buffer_grab(buf, pix_x, pix_y);
    300       pixel = tile_at(tile, x, y);
    301 
    302       res = setup_estimator_from_pixel(estimator, spp, pixel);
    303       if(res != RES_OK) goto error;
    304     }
    305   }
    306 
    307 exit:
    308   return res;
    309 error:
    310   goto exit;
    311 }
    312 
    313 static res_T
    314 write_list_of_tiles
    315   (struct sdis_estimator_buffer* buf,
    316    const size_t spp, /* #samples per pixel */
    317    struct list_node* tiles) /* Tiles to write */
    318 {
    319   struct list_node* node = NULL;
    320   res_T res = RES_OK;
    321   ASSERT(buf && spp && tiles);
    322 
    323   LIST_FOR_EACH(node, tiles) {
    324     struct tile* tile = CONTAINER_OF(node, struct tile, node);
    325     res = write_tile(buf, spp, tile);
    326     if(res != RES_OK) goto error;
    327   }
    328 
    329 exit:
    330   return res;
    331 error:
    332   goto exit;
    333 }
    334 
    335 #ifndef SDIS_ENABLE_MPI
    336 static INLINE res_T
    337 gather_tiles
    338   (struct sdis_device* dev,
    339    struct sdis_estimator_buffer* buf, /* NULL on non master processes */
    340    const size_t spp, /* #realisations per pixel */
    341    const size_t ntiles, /* Overall #tiles that must be written into `buf' */
    342    struct list_node* tiles) /* List of tiles of the current process */
    343 {
    344   (void)dev, (void)ntiles;
    345   return write_list_of_tiles(buf, spp, tiles);
    346 }
    347 #else
    348 /* Gather the tiles and write them into the estimator buffer. The master process
    349  * write its own tiles into the estimator buffer. Then, it gathers the tiles
    350  * send by non master processes and write them too into the estimator buffer */
    351 static res_T
    352 gather_tiles
    353   (struct sdis_device* dev,
    354    struct sdis_estimator_buffer* buf, /* NULL on non master processes */
    355    const size_t spp, /* #realisations per pixel */
    356    const size_t ntiles, /* Overall #tiles that must be written into `buf' */
    357    struct list_node* tiles) /* List of tiles of the current process */
    358 {
    359   struct tile* tile_temp = NULL;
    360   struct list_node* node = NULL;
    361   res_T res = RES_OK;
    362   ASSERT(dev && spp && tiles);
    363 
    364   if(!dev->use_mpi) {
    365     res = write_list_of_tiles(buf, spp, tiles);
    366     if(res != RES_OK) goto error;
    367     goto exit; /* No more to do */
    368   }
    369 
    370   /* Non master process */
    371   if(dev->mpi_rank != 0) {
    372 
    373     /* Send to the master process the list of tiles solved by the current
    374      * process */
    375     LIST_FOR_EACH(node, tiles) {
    376       struct tile* tile = CONTAINER_OF(node, struct tile, node);
    377       mutex_lock(dev->mpi_mutex);
    378       MPI(Send(&tile->data, sizeof(tile->data), MPI_CHAR, 0, MPI_SDIS_MSG_TILE,
    379         MPI_COMM_WORLD));
    380       mutex_unlock(dev->mpi_mutex);
    381     }
    382 
    383   /* Master process */
    384   } else {
    385     size_t itile;
    386     size_t ntiles_master;
    387 
    388     /* Write into the buffer the tiles solved by the master process itself */
    389     res = write_list_of_tiles(buf, spp, tiles);
    390     if(res != RES_OK) goto error;
    391 
    392     /* Create a temporary tile to store the tile sent by the non master
    393      * processes */
    394     res = tile_create(dev->allocator, &tile_temp);
    395     if(res != RES_OK) {
    396       log_err(dev,
    397         "Could not allocate the tile to temporary store the tiles sent by "
    398         "the non master processes -- %s", res_to_cstr(res));
    399       goto error;
    400     }
    401 
    402     /* Count the number of tiles rendered onto the master process */
    403     ntiles_master = 0;
    404     LIST_FOR_EACH(node, tiles) ++ntiles_master;
    405     ASSERT(ntiles_master <= ntiles);
    406 
    407     /* Receive the remaining tiles sent by the non master processes */
    408     FOR_EACH(itile, ntiles_master, ntiles) {
    409       MPI_Request req;
    410 
    411       /* Asynchronously receive a tile */
    412       mutex_lock(dev->mpi_mutex);
    413       MPI(Irecv(&tile_temp->data, sizeof(tile_temp->data), MPI_CHAR,
    414         MPI_ANY_SOURCE, MPI_SDIS_MSG_TILE, MPI_COMM_WORLD, &req));
    415       mutex_unlock(dev->mpi_mutex);
    416       mpi_waiting_for_request(dev, &req);
    417 
    418       /* Write the tile into the estimator buffer */
    419       res = write_tile(buf, spp, tile_temp);
    420       if(res != RES_OK) goto error;
    421     }
    422   }
    423 
    424 exit:
    425   if(tile_temp) tile_ref_put(tile_temp);
    426   return res;
    427 error:
    428   goto exit;
    429 }
    430 #endif
    431 
    432 /* Setup the accumulators of the whole estimator buffer */
    433 static res_T
    434 finalize_estimator_buffer
    435   (struct sdis_estimator_buffer* buf,
    436    struct ssp_rng_proxy* rng_proxy,
    437    const size_t spp) /* #samples per pixel */
    438 {
    439   struct accum acc_temp = ACCUM_NULL;
    440   struct accum acc_time = ACCUM_NULL;
    441   size_t definition[2];
    442   size_t x, y;
    443   size_t nrealisations = 0;
    444   size_t nsuccesses = 0;
    445   res_T res = RES_OK;
    446   ASSERT(buf && rng_proxy && spp);
    447 
    448   SDIS(estimator_buffer_get_definition(buf, definition));
    449 
    450   FOR_EACH(y, 0, definition[1]) {
    451     FOR_EACH(x, 0, definition[0]) {
    452       const struct sdis_estimator* estimator;
    453       SDIS(estimator_buffer_at(buf, x, y, &estimator));
    454       acc_temp.sum += estimator->temperature.sum;
    455       acc_temp.sum2 += estimator->temperature.sum2;
    456       acc_temp.count += estimator->temperature.count;
    457       acc_time.sum += estimator->realisation_time.sum;
    458       acc_time.sum2 += estimator->realisation_time.sum2;
    459       acc_time.count += estimator->realisation_time.count;
    460       nsuccesses += estimator->nrealisations;
    461     }
    462   }
    463 
    464   nrealisations = definition[0]*definition[1]*spp;
    465   ASSERT(acc_temp.count == acc_time.count);
    466   ASSERT(acc_temp.count == nsuccesses);
    467   estimator_buffer_setup_realisations_count(buf, nrealisations, nsuccesses);
    468   estimator_buffer_setup_temperature(buf, acc_temp.sum, acc_temp.sum2);
    469   estimator_buffer_setup_realisation_time(buf, acc_time.sum, acc_time.sum2);
    470   res = estimator_buffer_save_rng_state(buf, rng_proxy);
    471   if(res != RES_OK) goto error;
    472 
    473 exit:
    474   return res;
    475 error:
    476   goto exit;
    477 }
    478 
    479 static void
    480 free_rendered_tiles(struct list_node* tiles)
    481 {
    482   struct list_node* node;
    483   struct list_node* tmp;
    484   ASSERT(tiles);
    485   LIST_FOR_EACH_SAFE(node, tmp, tiles) {
    486     struct tile* tile = CONTAINER_OF(node, struct tile, node);
    487     list_del(node);
    488     tile_ref_put(tile);
    489   }
    490 }
    491 
    492 /*******************************************************************************
    493  * Exported function
    494  ******************************************************************************/
    495 res_T
    496 sdis_solve_camera
    497   (struct sdis_scene* scn,
    498    const struct sdis_solve_camera_args* args,
    499    struct sdis_estimator_buffer** out_buf)
    500 {
    501   /* Time registration */
    502   struct time time0, time1;
    503   char buffer[128]; /* Temporary buffer used to store formated time */
    504 
    505   /* Stardis variables */
    506   struct sdis_estimator_buffer* buf = NULL;
    507 
    508   /* Random number generators */
    509   struct ssp_rng_proxy* rng_proxy = NULL;
    510   struct ssp_rng** per_thread_rng = NULL;
    511 
    512   /* Enclosure & medium in which the probe lies */
    513   unsigned enc_id = ENCLOSURE_ID_NULL;
    514 
    515   /* Miscellaneous */
    516   size_t ntiles_x, ntiles_y, ntiles, ntiles_adjusted;
    517   size_t ntiles_proc; /* #tiles for the current proc */
    518   struct list_node tiles; /* List of tiles rendered by the process */
    519   double pix_sz[2]; /* Size of a pixel in the normalized image plane */
    520   int64_t mcode; /* Morton code of a tile */
    521   int64_t mcode_1st; /* morton code of the 1st tile computed by the process */
    522   int64_t mcode_incr; /* Increment toward the next morton code */
    523   int32_t* progress = NULL; /* Per process progress bar */
    524   int pcent_progress = 1; /* Percentage requiring progress update */
    525   int register_paths = SDIS_HEAT_PATH_NONE;
    526   int is_master_process = 1;
    527   ATOMIC nsolved_tiles = 0;
    528   ATOMIC res = RES_OK;
    529 
    530   list_init(&tiles);
    531 
    532   if(!scn || !out_buf) { res = RES_BAD_ARG; goto error; }
    533   res = check_solve_camera_args(args);
    534   if(res != RES_OK) goto error;
    535 
    536   if(scene_is_2d(scn)) {
    537     log_err(scn->dev, "%s: 2D scenes are not supported.\n", FUNC_NAME);
    538     goto error;
    539   }
    540 
    541   /* Retrieve the medium in which the submitted position lies */
    542   res = scene_get_enclosure_id(scn, args->cam->position, &enc_id);
    543   if(res != RES_OK) goto error;
    544 
    545   /* Create the per thread RNGs */
    546   res = create_per_thread_rng
    547     (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng);
    548   if(res != RES_OK) goto error;
    549 
    550   /* Allocate the per process progress status */
    551   res = alloc_process_progress(scn->dev, &progress);
    552   if(res != RES_OK) goto error;
    553 
    554   /* Compute the overall number of tiles */
    555   ntiles_x = (args->image_definition[0] + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
    556   ntiles_y = (args->image_definition[1] + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
    557   ntiles = ntiles_x * ntiles_y;
    558   ntiles_adjusted = round_up_pow2(MMAX(ntiles_x, ntiles_y));
    559   ntiles_adjusted *= ntiles_adjusted;
    560 
    561 #ifdef SDIS_ENABLE_MPI
    562   is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0;
    563   if(scn->dev->use_mpi) {
    564     mcode_1st = scn->dev->mpi_rank;
    565     mcode_incr = scn->dev->mpi_nprocs;
    566 
    567     /* Compute the #tiles of the current proc */
    568     ntiles_proc = ntiles / (size_t)scn->dev->mpi_nprocs;
    569     if(ntiles % (size_t)scn->dev->mpi_nprocs > (size_t)scn->dev->mpi_rank) {
    570       ++ntiles_proc;
    571     }
    572   } else
    573 #endif
    574   {
    575     is_master_process = 1;
    576     mcode_1st = 0;
    577     mcode_incr = 1;
    578     ntiles_proc = ntiles;
    579   }
    580 
    581   /* Update the progress bar every percent if escape sequences are allowed in
    582    * log messages or only every 10 percent when only plain text is allowed.
    583    * This reduces the number of lines of plain text printed */
    584   pcent_progress = scn->dev->no_escape_sequence ? 10 : 1;
    585 
    586   /* Compute the normalized pixel size */
    587   pix_sz[0] = 1.0 / (double)args->image_definition[0];
    588   pix_sz[1] = 1.0 / (double)args->image_definition[1];
    589 
    590   /* Create the global estimator on the master process only */
    591   if(is_master_process) {
    592     res = estimator_buffer_create
    593       (scn->dev, args->image_definition[0], args->image_definition[1], &buf);
    594     if(res != RES_OK) goto error;
    595   }
    596 
    597   /* Synchronise the processes */
    598   process_barrier(scn->dev);
    599 
    600   #define PROGRESS_MSG "Rendering: "
    601   print_progress(scn->dev, progress, PROGRESS_MSG);
    602 
    603   /* Begin time registration of the computation */
    604   time_current(&time0);
    605 
    606   /* Here we go! Launch the Monte Carlo estimation */
    607   omp_set_num_threads((int)scn->dev->nthreads);
    608   register_paths = is_master_process
    609     ? args->register_paths : SDIS_HEAT_PATH_NONE;
    610   #pragma omp parallel for schedule(static, 1/*chunk size*/)
    611   for(mcode = mcode_1st; mcode < (int64_t)ntiles_adjusted; mcode+=mcode_incr) {
    612     size_t tile_org[2] = {0, 0};
    613     size_t tile_sz[2] = {0, 0};
    614     struct tile* tile = NULL;
    615     const int ithread = omp_get_thread_num();
    616     struct ssp_rng* rng = per_thread_rng[ithread];
    617     size_t n;
    618     int pcent;
    619     res_T res_local = RES_OK;
    620 
    621     if(ATOMIC_GET(&res) != RES_OK) continue;
    622 
    623     tile_org[0] = morton2D_decode_u16((uint32_t)(mcode>>0));
    624     if(tile_org[0] >= ntiles_x) continue; /* Discard tile */
    625     tile_org[1] = morton2D_decode_u16((uint32_t)(mcode>>1));
    626     if(tile_org[1] >= ntiles_y) continue; /* Discard tile */
    627 
    628     res_local = tile_create(scn->dev->allocator, &tile);
    629     if(tile == NULL) {
    630       log_err(scn->dev, "%s: error allocating the tile (%lu, %lu) -- %s\n",
    631         FUNC_NAME,
    632         (unsigned long)tile_org[0],
    633         (unsigned long)tile_org[1],
    634         res_to_cstr(res_local));
    635       ATOMIC_SET(&res, res_local);
    636       continue;
    637     }
    638 
    639     /* Register the tile */
    640     #pragma omp critical
    641     list_add_tail(&tiles, &tile->node);
    642 
    643     /* Setup the tile coordinates */
    644     tile->data.x = (uint16_t)tile_org[0];
    645     tile->data.y = (uint16_t)tile_org[1];
    646 
    647     /* Setup the tile coordinates in the image plane */
    648     tile_org[0] *= TILE_SIZE;
    649     tile_org[1] *= TILE_SIZE;
    650     tile_sz[0] = MMIN(TILE_SIZE, args->image_definition[0] - tile_org[0]);
    651     tile_sz[1] = MMIN(TILE_SIZE, args->image_definition[1] - tile_org[1]);
    652 
    653     /* Draw the tile */
    654     res_local = solve_tile
    655       (scn, rng, enc_id, args->cam, args->time_range, tile_org, tile_sz,
    656        args->spp, register_paths, pix_sz, args->picard_order, args->diff_algo,
    657        buf, tile);
    658     if(res_local != RES_OK) {
    659       ATOMIC_SET(&res, res_local);
    660       continue;
    661     }
    662 
    663     /* Update progress */
    664     n = (size_t)ATOMIC_INCR(&nsolved_tiles);
    665     pcent = (int)((double)n*100.0 / (double)ntiles_proc + 0.5/*round*/);
    666     #pragma omp critical
    667     if(pcent/pcent_progress > progress[0]/pcent_progress) {
    668       progress[0] = pcent;
    669       print_progress_update(scn->dev, progress, PROGRESS_MSG);
    670     }
    671   }
    672 
    673   /* Synchronise the processes */
    674   process_barrier(scn->dev);
    675 
    676   res = gather_res_T(scn->dev, (res_T)res);
    677   if(res != RES_OK) goto error;
    678 
    679   print_progress_completion(scn->dev, progress, PROGRESS_MSG);
    680   #undef PROGRESS_MSG
    681 
    682   /* Report computation time */
    683   time_sub(&time0, time_current(&time1), &time0);
    684   time_dump(&time0, TIME_ALL, NULL, buffer, sizeof(buffer));
    685   log_info(scn->dev, "Image rendered in %s.\n", buffer);
    686 
    687   /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of
    688    * the master process is greater than the RNG proxy state of all other
    689    * processes */
    690   res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy);
    691   if(res != RES_OK) goto error;
    692 
    693   time_current(&time0);
    694 
    695   res = gather_tiles(scn->dev, buf, args->spp, ntiles, &tiles);
    696   if(res != RES_OK) goto error;
    697 
    698   time_sub(&time0, time_current(&time1), &time0);
    699   time_dump(&time0, TIME_ALL, NULL, buffer, sizeof(buffer));
    700   log_info(scn->dev, "Image tiles gathered in %s.\n", buffer);
    701 
    702   if(is_master_process) {
    703     res = finalize_estimator_buffer(buf, rng_proxy, args->spp);
    704     if(res != RES_OK) goto error;
    705   }
    706 
    707 exit:
    708   free_rendered_tiles(&tiles);
    709   if(per_thread_rng)release_per_thread_rng(scn->dev, per_thread_rng);
    710   if(progress) free_process_progress(scn->dev, progress);
    711   if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
    712   if(out_buf) *out_buf = buf;
    713   return (res_T)res;
    714 error:
    715   if(buf) { SDIS(estimator_buffer_ref_put(buf)); buf = NULL; }
    716   goto exit;
    717 }