htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

commit b063c15a54b6586cf520999e7b908924b588a6fd
parent 09b8fad5ea39a3fc20b3feb170f8b7f1e77b960f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  5 Mar 2021 16:10:15 +0100

Begin the implementation of the combustion draw map function

Diffstat:
Mcmake/combustion/CMakeLists.txt | 2++
Msrc/atmosphere/htrdr_atmosphere_draw_map.c | 2+-
Msrc/combustion/htrdr_combustion.c | 6++++--
Msrc/combustion/htrdr_combustion_c.h | 19+++++++++++++++++--
Asrc/combustion/htrdr_combustion_compute_radiance_sw.c | 55+++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/combustion/htrdr_combustion_draw_map.c | 186+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/core/htrdr_draw_map.c | 2+-
Msrc/core/htrdr_draw_map.h | 2--
8 files changed, 266 insertions(+), 8 deletions(-)

diff --git a/cmake/combustion/CMakeLists.txt b/cmake/combustion/CMakeLists.txt @@ -49,6 +49,8 @@ include_directories( set(HTRDR_COMBUSTION_FILES_SRC htrdr_combustion.c htrdr_combustion_args.c + htrdr_combustion_draw_map.c + htrdr_combustion_compute_radiance_sw.c htrdr_combustion_main.c) set(HTRDR_COMBUSTION_FILES_INC diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c @@ -395,7 +395,7 @@ draw_pixel_xwave time_sub(&t0, time_current(&t1), &t0); usec = (double)time_val(&t0, TIME_NSEC) * 0.001; - /* Update the pixel accumulator of the current channel */ + /* Update the pixel accumulator */ radiance.sum_weights += weight; radiance.sum_weights_sqr += weight*weight; radiance.nweights += 1; diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c @@ -166,7 +166,7 @@ setup_buffer /* Create the image buffer only on the master process; the image parts * rendered by the others processes are gathered onto the master process */ - if(!htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; + if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; res = htrdr_buffer_create(cmd->htrdr, &cmd->buf_layout, &cmd->buf); if(res != RES_OK) goto error; @@ -300,6 +300,7 @@ htrdr_combustion_create htrdr_ref_get(htrdr); cmd->htrdr = htrdr; + cmd->spp = args->image.spp; cmd->dump_volumetric_acceleration_structure = args->dump_volumetric_acceleration_structure; @@ -351,7 +352,8 @@ htrdr_combustion_run(struct htrdr_combustion* cmd) res = dump_volumetric_acceleration_structure(cmd); if(res != RES_OK) goto error; } else { - /* TODO */ + res = combustion_draw_map(cmd); + if(res != RES_OK) goto error; } exit: diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h @@ -32,6 +32,7 @@ struct htrdr_camera; struct htrdr_geometry; struct htrdr_materials; struct htrdr_rectangle; +struct ssp_rng; struct combustion_pixel { struct htrdr_estimate radiance; /* In W/m^2/sr */ @@ -53,13 +54,13 @@ struct htrdr_combustion { struct htrdr_rectangle* laser; /* Laser surface emission */ double wavelength; /* Wavelength of the laser in nanometer */ - struct htrdr_args_image image; /* Image definition */ struct htrdr_buffer_layout buf_layout; struct htrdr_buffer* buf; /* NULL on non master processes */ + size_t spp; /* #samples per pixel */ FILE* output; /* Output stream */ struct str output_name; /* Name of the output stream */ - int dump_volumetric_acceleration_structure; + int dump_volumetric_acceleration_structure; ref_T ref; struct htrdr* htrdr; @@ -70,4 +71,18 @@ combustion_get_pixel_format (const struct htrdr_combustion* cmd, struct htrdr_pixel_format* fmt); +extern LOCAL_SYM res_T +combustion_draw_map + (struct htrdr_combustion* cmd); + +/* Return the shortwave radiance in W/m^2/sr/m */ +extern LOCAL_SYM double +combustion_compute_radiance_sw + (struct htrdr_combustion* cmd, + const size_t ithread, + struct ssp_rng* rng, + const double pos_in[3], + const double dir_in[3], + const double wlen); /* In nanometer */ + #endif /* HTRDR_COMBUSTION_C_H */ diff --git a/src/combustion/htrdr_combustion_compute_radiance_sw.c b/src/combustion/htrdr_combustion_compute_radiance_sw.c @@ -0,0 +1,55 @@ +/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com) + * Copyright (C) 2018, 2019, 2021 CNRS + * Copyright (C) 2018, 2019, 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/>. */ + +#include "htrdr_combustion_c.h" + +#include <rsys/double2.h> +#include <rsys/double3.h> + +/******************************************************************************* + * Local functions + ******************************************************************************/ +extern LOCAL_SYM double +combustion_compute_radiance_sw + (struct htrdr_combustion* cmd, + const size_t ithread, + struct ssp_rng* rng, + const double pos_in[3], + const double dir_in[3], + const double wlen) /* In nanometer */ +{ + double pos[3]; + double dir[3]; + double range[2]; + double ksi; /* Throughput */ + double weight; + ASSERT(cmd && rng && pos_in && dir_in); + (void)cmd, (void)ithread, (void)rng, (void)pos_in, (void)dir_in, (void)wlen; + (void)ksi; + + d3_set(pos, pos_in); + d3_set(dir, dir_in); + d2(range, 0, FLT_MAX); + + ksi = 1; + weight = 0; + + /* TODO Radiative random walk */ + + return weight; +} + diff --git a/src/combustion/htrdr_combustion_draw_map.c b/src/combustion/htrdr_combustion_draw_map.c @@ -0,0 +1,186 @@ +/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com) + * Copyright (C) 2018, 2019, 2021 CNRS + * Copyright (C) 2018, 2019, 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/>. */ + +#include "combustion/htrdr_combustion_c.h" + +#include "core/htrdr_accum.h" +#include "core/htrdr_camera.h" +#include "core/htrdr_draw_map.h" +#include "core/htrdr_log.h" + +#include <star/ssp.h> + +#include <rsys/clock_time.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +draw_pixel + (struct htrdr* htrdr, + const struct htrdr_draw_pixel_args* args, + void* data) +{ + struct htrdr_accum radiance = HTRDR_ACCUM_NULL; + struct htrdr_accum time = HTRDR_ACCUM_NULL; + struct htrdr_combustion* cmd = NULL; + struct combustion_pixel* pixel = NULL; + size_t isamp; + ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); + (void)htrdr; + + cmd = args->context; + pixel = data; + + FOR_EACH(isamp, 0, args->spp) { + struct time t0, t1; + double pix_samp[2]; + double ray_org[3]; + double ray_dir[3]; + double weight; + double usec; + + /* Begin the registration of the time spent in the realisation */ + time_current(&t0); + + /* Sample a position into the pixel, in the normalized image plane */ + pix_samp[0] = (double)args->pixel_coord[0] + ssp_rng_canonical(args->rng); + pix_samp[1] = (double)args->pixel_coord[1] + ssp_rng_canonical(args->rng); + pix_samp[0] *= args->pixel_normalized_size[0]; + pix_samp[1] *= args->pixel_normalized_size[1]; + + /* Generate a ray starting from the pinhole camera and passing through the + * pixel sample */ + htrdr_camera_ray(cmd->camera, pix_samp, ray_org, ray_dir); + + /* Backward trace the path */ + weight = combustion_compute_radiance_sw(cmd, args->ithread, args->rng, + ray_org, ray_dir, cmd->wavelength); + + /* End the registration of the per realisation time */ + time_sub(&t0, time_current(&t1), &t0); + usec = (double)time_val(&t0, TIME_NSEC) * 0.001; + + /* Update the pixel accumulator */ + radiance.sum_weights += weight; + radiance.sum_weights_sqr += weight*weight; + radiance.nweights += 1; + + /* Update the pixel accumulator of per realisation time */ + time.sum_weights += usec; + time.sum_weights_sqr += usec*usec; + time.nweights += 1; + } + + /* Compute the estimation of the pixel radiance */ + htrdr_accum_get_estimation(&radiance, &pixel->radiance); + + /* Save the per realisation time */ + pixel->time = time; +} + +static void +dump_pixel + (const struct combustion_pixel* pix, + struct htrdr_accum* time_acc, /* May be NULL */ + FILE* stream) +{ + struct htrdr_estimate time = HTRDR_ESTIMATE_NULL; + ASSERT(pix && stream && pix->time.nweights); + + htrdr_accum_get_estimation(&pix->time, &time); + + fprintf(stream, "%g %g 0 0 0 0 %g %g\n", + pix->radiance.E, pix->radiance.SE, time.E,time.SE); + + if(time_acc) { + time_acc->sum_weights += pix->time.sum_weights; + time_acc->sum_weights_sqr += pix->time.sum_weights_sqr; + time_acc->nweights += pix->time.nweights; + } +} + +static res_T +dump_buffer + (struct htrdr_combustion* cmd, + struct htrdr_buffer* buf, + struct htrdr_accum* time_acc, /* May be NULL */ + FILE* stream) +{ + struct htrdr_pixel_format pixfmt; + struct htrdr_buffer_layout layout; + size_t x, y; + ASSERT(cmd && buf && stream); + + combustion_get_pixel_format(cmd, &pixfmt); + htrdr_buffer_get_layout(buf, &layout); + CHK(pixfmt.size == layout.elmt_size); + + fprintf(stream, "%lu %lu\n", layout.width, layout.height); + + if(time_acc) *time_acc = HTRDR_ACCUM_NULL; + + FOR_EACH(y, 0, layout.height) { + FOR_EACH(x, 0, layout.width) { + const struct combustion_pixel* pix = htrdr_buffer_at(buf, x, y); + ASSERT(IS_ALIGNED(pix, pixfmt.alignment)); + dump_pixel(pix, time_acc, stream); + } + fprintf(stream, "\n"); + } + return RES_OK; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +combustion_draw_map(struct htrdr_combustion* cmd) +{ + struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL; + struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL; + struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL; + res_T res = RES_OK; + ASSERT(cmd); + + /* Setup the draw map input arguments */ + args.draw_pixel = draw_pixel; + args.buffer_layout = cmd->buf_layout; + args.spp = cmd->spp; + args.context = cmd; + + res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf); + if(res != RES_OK) goto error; + + /* No more to do on non master processes */ + if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; + + /* Write buffer to output */ + res = dump_buffer(cmd, cmd->buf, &path_time_acc, cmd->output); + if(res != RES_OK) goto error; + + htrdr_accum_get_estimation(&path_time_acc, &path_time); + htrdr_log(cmd->htrdr, + "Time per radiative path (in micro seconds): %g +/- %g\n", + path_time.E, path_time.SE); + +exit: + return res; +error: + goto exit; +} + diff --git a/src/core/htrdr_draw_map.c b/src/core/htrdr_draw_map.c @@ -635,7 +635,7 @@ draw_map * is initialised in order to ensure that an unique and predictable set of * random numbers is used for the current tile. */ SSP(rng_proxy_create2 - (&htrdr->lifo_allocators[ithread], + (htrdr_get_thread_allocator(htrdr, (size_t)ithread), &ssp_rng_threefry, RNG_SEQUENCE_SIZE * (size_t)mcode, /* Offset */ RNG_SEQUENCE_SIZE, /* Size */ diff --git a/src/core/htrdr_draw_map.h b/src/core/htrdr_draw_map.h @@ -52,7 +52,6 @@ typedef void struct htrdr_draw_map_args { htrdr_draw_pixel_T draw_pixel; struct htrdr_buffer_layout buffer_layout; - struct ssp_rng* rng; size_t spp; /* Samples per pixel */ void* context; /* User defined data */ }; @@ -60,7 +59,6 @@ struct htrdr_draw_map_args { #define HTRDR_DRAW_MAP_ARGS_NULL__ { \ NULL, /* Draw pixel functor */ \ HTRDR_BUFFER_LAYOUT_NULL__, /* Layout of the destination buffer */ \ - NULL, /* Random Number Generator */ \ 0, /* #Samples per pixel */ \ NULL /* User defined data */ \ }