stardis-solver

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

commit 8bcedc68395190ee034c13cbe2fa071d0c5ce5e7
parent 267793ef65cfbc627c56e5f4b986dc834be91c27
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 20 May 2019 16:14:40 +0200

Major update of the camera solver API

Replace the user defined buffer of accumulation by a buffer of
estimators. Remove the sdis_accum_buffer API.

Diffstat:
Mcmake/CMakeLists.txt | 3+--
Msrc/sdis.h | 83++++++++++++++++---------------------------------------------------------------
Dsrc/sdis_accum_buffer.c | 160-------------------------------------------------------------------------------
Msrc/sdis_device.c | 9---------
Msrc/sdis_device_c.h | 13-------------
Msrc/sdis_estimator.c | 19++++++++++++++-----
Asrc/sdis_estimator_buffer.c | 159+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sdis_estimator_buffer_c.h | 40++++++++++++++++++++++++++++++++++++++++
Msrc/sdis_estimator_c.h | 28+++++++++++++---------------
Msrc/sdis_solve.c | 117+++++++++++++++++++-------------------------------------------------------------
Msrc/test_sdis_solve_camera.c | 107++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
11 files changed, 333 insertions(+), 405 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -65,11 +65,11 @@ set(VERSION_PATCH 1) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SDIS_FILES_SRC - sdis_accum_buffer.c sdis_camera.c sdis_data.c sdis_device.c sdis_estimator.c + sdis_estimator_buffer.c sdis_green.c sdis_heat_path.c sdis_interface.c @@ -171,7 +171,6 @@ if(NOT NO_TEST) register_test(${_name} ${_name}) endfunction() - new_test(test_sdis_accum_buffer) new_test(test_sdis_camera) new_test(test_sdis_conducto_radiative) new_test(test_sdis_conducto_radiative_2d) diff --git a/src/sdis.h b/src/sdis.h @@ -56,11 +56,11 @@ struct senc_descriptor; * get or release a reference on the data, i.e. they increment or decrement the * reference counter, respectively. When this counter reaches 0, the object is * silently destroyed and cannot be used anymore. */ -struct sdis_accum_buffer; struct sdis_camera; struct sdis_data; struct sdis_device; struct sdis_estimator; +struct sdis_estimator_buffer; struct sdis_green_function; struct sdis_interface; struct sdis_medium; @@ -137,14 +137,6 @@ struct sdis_interface_fragment { static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL = SDIS_INTERFACE_FRAGMENT_NULL__; -/* Monte-Carlo accumulator */ -struct sdis_accum { - double sum_weights; /* Sum of Monte-Carlo weights */ - double sum_weights_sqr; /* Sum of Monte-Carlo square weights */ - size_t nweights; /* #accumulated weights */ - size_t nfailures; /* #failures */ -}; - /* Monte-Carlo estimation */ struct sdis_mc { double E; /* Expected value */ @@ -244,22 +236,6 @@ struct sdis_interface_shader { static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL = SDIS_INTERFACE_SHADER_NULL__; -struct sdis_accum_buffer_layout { - size_t width; - size_t height; -}; -#define SDIS_ACCUM_BUFFER_LAYOUT_NULL__ {0, 0} -static const struct sdis_accum_buffer_layout SDIS_ACCUM_BUFFER_LAYOUT_NULL = - SDIS_ACCUM_BUFFER_LAYOUT_NULL__; - -/* Functor use to write accumulations performed by sdis_solve_camera */ -typedef res_T -(*sdis_write_accums_T) - (void* context, /* User data */ - const size_t origin[2], /* Coordinates of the 1st accumulation in image plane */ - const size_t naccums[2], /* #accumulations in X and Y */ - const struct sdis_accum* accums); /* List of row ordered accumulations */ - /* Vertex of heat path v*/ struct sdis_heat_vertex { double P[3]; @@ -418,51 +394,27 @@ sdis_camera_look_at const double up[3]); /******************************************************************************* - * A buffer of accumulations is a 2D array whose each cell stores an - * Monte-Carlo accumulation, i.e. a sum of MC weights, the sum of their square - * and the overall number of summed weights (see struct sdis_accum) - ******************************************************************************/ + * An estimator buffer is 2D array of estimators +******************************************************************************/ SDIS_API res_T -sdis_accum_buffer_create - (struct sdis_device* dev, - const size_t width, /* #cells in X */ - const size_t height, /* #cells in Y */ - struct sdis_accum_buffer** buf); +sdis_estimator_buffer_ref_get + (struct sdis_estimator_buffer* buf); SDIS_API res_T -sdis_accum_buffer_ref_get - (struct sdis_accum_buffer* buf); +sdis_estimator_buffer_ref_put + (struct sdis_estimator_buffer* buf); SDIS_API res_T -sdis_accum_buffer_ref_put - (struct sdis_accum_buffer* buf); +sdis_estimator_buffer_get_definition + (const struct sdis_estimator_buffer* buf, + size_t definition[2]); SDIS_API res_T -sdis_accum_buffer_get_layout - (const struct sdis_accum_buffer* buf, - struct sdis_accum_buffer_layout* layout); - -/* Get a read only pointer toward the memory space of the accum buffer */ -SDIS_API res_T -sdis_accum_buffer_map - (const struct sdis_accum_buffer* buf, - const struct sdis_accum** accums); - -SDIS_API res_T -sdis_accum_buffer_unmap - (const struct sdis_accum_buffer* buf); - -/* Helper function that matches the `sdis_write_accums_T' functor type. One can - * send this function directly to the sdis_solve_camera function, to fill the - * accum buffer with the estimation of the radiative temperature that reaches - * each pixel of an image whose definition matches the definition of the accum - * buffer. */ -SDIS_API res_T -sdis_accum_buffer_write - (void* buf, /* User data */ - const size_t origin[2], /* Coordinates of the 1st accum in image plane */ - const size_t naccum[2], /* #accum in X and Y */ - const struct sdis_accum* accums); /* List of row ordered accum */ +sdis_estimator_buffer_at + (const struct sdis_estimator_buffer* buf, + const size_t x, + const size_t y, + const struct sdis_estimator** estimator); /******************************************************************************* * A medium encapsulates the properties of either a fluid or a solid. @@ -927,10 +879,7 @@ sdis_solve_camera const size_t height, /* Image definition in Y */ const size_t spp, /* #samples per pixel */ const int register_paths, /* Combination of enum sdis_heat_path_flag */ - sdis_write_accums_T writer, - void* writer_data, - /* Estimator of the whole image. May be NULL */ - struct sdis_estimator** estimator); + struct sdis_estimator_buffer** buf); SDIS_API res_T sdis_solve_medium diff --git a/src/sdis_accum_buffer.c b/src/sdis_accum_buffer.c @@ -1,160 +0,0 @@ -/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com) - * - * 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 "sdis.h" -#include "sdis_device_c.h" - -struct sdis_accum_buffer { - struct sdis_accum* accums; - size_t width; - size_t height; - - ref_T ref; - struct sdis_device* dev; -}; - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static void -accum_buffer_release(ref_T* ref) -{ - struct sdis_accum_buffer* buf; - struct sdis_device* dev; - ASSERT(ref); - buf = CONTAINER_OF(ref, struct sdis_accum_buffer, ref); - dev = buf->dev; - if(buf->accums) MEM_RM(dev->allocator, buf->accums); - MEM_RM(dev->allocator, buf); - SDIS(device_ref_put(dev)); -} - -/******************************************************************************* - * Export functions - ******************************************************************************/ -res_T -sdis_accum_buffer_create - (struct sdis_device* dev, - const size_t width, - const size_t height, - struct sdis_accum_buffer** out_buf) -{ - struct sdis_accum_buffer* buf = NULL; - size_t sz; - res_T res = RES_OK; - - if(!dev || !width || !height || !out_buf) { - res = RES_BAD_ARG; - goto error; - } - - buf = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_accum_buffer)); - if(!buf) { - res = RES_MEM_ERR; - goto error; - } - ref_init(&buf->ref); - SDIS(device_ref_get(dev)); - buf->dev = dev; - buf->width = width; - buf->height = height; - - sz = buf->width * buf->height * sizeof(struct sdis_accum); - buf->accums = MEM_ALLOC_ALIGNED(dev->allocator, sz, 16); - if(!buf->accums) { - log_err(dev, "%s: could not allocate a buffer accumulations of %lux%lu.\n", - FUNC_NAME, (unsigned long)width, (unsigned long)height); - res = RES_MEM_ERR; - goto error; - } - -exit: - if(out_buf) *out_buf = buf; - return res; -error: - if(buf) { - SDIS(accum_buffer_ref_put(buf)); - buf = NULL; - } - goto exit; -} - -res_T -sdis_accum_buffer_ref_get(struct sdis_accum_buffer* buf) -{ - if(!buf) return RES_BAD_ARG; - ref_get(&buf->ref); - return RES_OK; -} - -res_T -sdis_accum_buffer_ref_put(struct sdis_accum_buffer* buf) -{ - if(!buf) return RES_BAD_ARG; - ref_put(&buf->ref, accum_buffer_release); - return RES_OK; -} - -res_T -sdis_accum_buffer_get_layout - (const struct sdis_accum_buffer* buf, - struct sdis_accum_buffer_layout* layout) -{ - if(!buf || !layout) return RES_BAD_ARG; - layout->width = buf->width; - layout->height = buf->height; - return RES_OK; -} - -res_T -sdis_accum_buffer_map - (const struct sdis_accum_buffer* buf, const struct sdis_accum** accums) -{ - if(!buf || !accums) return RES_BAD_ARG; - *accums = buf->accums; - return RES_OK; -} - -res_T -sdis_accum_buffer_unmap(const struct sdis_accum_buffer* buf) -{ - if(!buf) return RES_BAD_ARG; - /* Do nothing */ - return RES_OK; -} - -res_T -sdis_accum_buffer_write - (void* buffer, - const size_t origin[2], - const size_t size[2], - const struct sdis_accum* accums) -{ - struct sdis_accum_buffer* buf = buffer; - const struct sdis_accum* src_row = accums; - size_t dst_iy; - - if(UNLIKELY(!buffer || !origin || !size || !accums)) return RES_BAD_ARG; - if(UNLIKELY((origin[0] + size[0]) > buf->width)) return RES_BAD_ARG; - if(UNLIKELY((origin[1] + size[1]) > buf->height)) return RES_BAD_ARG; - - FOR_EACH(dst_iy, origin[1], origin[1] + size[1]) { - const size_t dst_irow = dst_iy*buf->width + origin[0]; - memcpy(buf->accums + dst_irow, src_row, size[0] * sizeof(struct sdis_accum)); - src_row += size[0]; - } - return RES_OK; -} - diff --git a/src/sdis_device.c b/src/sdis_device.c @@ -54,7 +54,6 @@ device_release(ref_T* ref) ASSERT(flist_name_is_empty(&dev->media_names)); flist_name_release(&dev->interfaces_names); flist_name_release(&dev->media_names); - darray_tile_release(&dev->tiles); MEM_RM(dev->allocator, dev); } @@ -98,14 +97,6 @@ sdis_device_create ref_init(&dev->ref); flist_name_init(allocator, &dev->interfaces_names); flist_name_init(allocator, &dev->media_names); - darray_tile_init(allocator, &dev->tiles); - - res = darray_tile_resize(&dev->tiles, dev->nthreads); - if(res != RES_OK) { - log_err(dev, - "%s: could not allocate the per thread buffer of estimations.\n", FUNC_NAME); - goto error; - } res = s2d_device_create(log, allocator, 0, &dev->s2d); if(res != RES_OK) { diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h @@ -26,18 +26,6 @@ struct name { FITEM; }; #define FITEM_TYPE name #include <rsys/free_list.h> -#define DARRAY_NAME accum -#define DARRAY_DATA struct sdis_accum -#include <rsys/dynamic_array.h> - -#define DARRAY_NAME tile -#define DARRAY_DATA struct darray_accum -#define DARRAY_FUNCTOR_INIT darray_accum_init -#define DARRAY_FUNCTOR_RELEASE darray_accum_release -#define DARRAY_FUNCTOR_COPY darray_accum_copy -#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_accum_copy_and_release -#include <rsys/dynamic_array.h> - struct sdis_device { struct logger* logger; struct mem_allocator* allocator; @@ -46,7 +34,6 @@ struct sdis_device { struct flist_name interfaces_names; struct flist_name media_names; - struct darray_tile tiles; struct s2d_device* s2d; struct s3d_device* s3d; diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -82,12 +82,19 @@ sdis_estimator_get_failure_count return RES_OK; } +#define SETUP_MC(McName, Acc) { \ + (McName)->E = (Acc)->sum / (double)(Acc)->count; \ + (McName)->V = (Acc)->sum2 / (double)(Acc)->count - (McName)->E*(McName)->E; \ + (McName)->V = MMAX((McName)->V, 0); \ + (McName)->SE = sqrt((McName)->V / (double)(Acc)->count); \ +} (void)0 + res_T sdis_estimator_get_temperature (const struct sdis_estimator* estimator, struct sdis_mc* mc) { if(!estimator || !mc) return RES_BAD_ARG; - *mc = estimator->temperature; + SETUP_MC(mc, &estimator->temperature); return RES_OK; } @@ -96,7 +103,7 @@ sdis_estimator_get_realisation_time (const struct sdis_estimator* estimator, struct sdis_mc* mc) { if(!estimator || !mc) return RES_BAD_ARG; - *mc = estimator->realisation_time; + SETUP_MC(mc, &estimator->realisation_time); return RES_OK; } @@ -107,7 +114,7 @@ sdis_estimator_get_convective_flux if(!estimator || !flux ||estimator->type != SDIS_ESTIMATOR_FLUX) return RES_BAD_ARG; ASSERT(estimator->fluxes); - *flux = estimator->fluxes[FLUX_CONVECTIVE]; + SETUP_MC(flux, &estimator->fluxes[FLUX_CONVECTIVE]); return RES_OK; } @@ -118,7 +125,7 @@ sdis_estimator_get_radiative_flux if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX) return RES_BAD_ARG; ASSERT(estimator->fluxes); - *flux = estimator->fluxes[FLUX_RADIATIVE]; + SETUP_MC(flux, &estimator->fluxes[FLUX_RADIATIVE]); return RES_OK; } @@ -129,10 +136,12 @@ sdis_estimator_get_total_flux if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX) return RES_BAD_ARG; ASSERT(estimator->fluxes); - *flux = estimator->fluxes[FLUX_TOTAL]; + SETUP_MC(flux, &estimator->fluxes[FLUX_TOTAL]); return RES_OK; } +#undef SETUP_MC + res_T sdis_estimator_get_paths_count (const struct sdis_estimator* estimator, size_t* npaths) diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c @@ -0,0 +1,159 @@ +/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com) + * + * 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 "sdis.h" +#include "sdis_device_c.h" +#include "sdis_estimator_c.h" +#include "sdis_estimator_buffer_c.h" + +struct sdis_estimator_buffer { + struct sdis_estimator** estimators; /* Row major estimators list */ + size_t width; + size_t height; + + ref_T ref; + struct sdis_device* dev; +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +estimator_buffer_release(ref_T* ref) +{ + struct sdis_estimator_buffer* buf; + struct sdis_device* dev; + size_t i; + ASSERT(ref); + + buf = CONTAINER_OF(ref, struct sdis_estimator_buffer, ref); + dev = buf->dev; + + if(buf->estimators) { + FOR_EACH(i, 0, buf->width*buf->height) { + if(buf->estimators[i]) SDIS(estimator_ref_put(buf->estimators[i])); + } + MEM_RM(dev->allocator, buf->estimators); + } + + MEM_RM(dev->allocator, buf); + SDIS(device_ref_put(dev)); +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +sdis_estimator_buffer_ref_get(struct sdis_estimator_buffer* buf) +{ + if(!buf) return RES_BAD_ARG; + ref_get(&buf->ref); + return RES_OK; +} + +res_T +sdis_estimator_buffer_ref_put(struct sdis_estimator_buffer* buf) +{ + if(!buf) return RES_BAD_ARG; + ref_put(&buf->ref, estimator_buffer_release); + return RES_OK; +} + +res_T +sdis_estimator_buffer_get_definition + (const struct sdis_estimator_buffer* buf, size_t definition[2]) +{ + if(!buf || !definition) return RES_BAD_ARG; + definition[0] = buf->width; + definition[1] = buf->height; + return RES_OK; +} + +res_T +sdis_estimator_buffer_at + (const struct sdis_estimator_buffer* buf, + const size_t x, + const size_t y, + const struct sdis_estimator** estimator) +{ + if(!buf || x >= buf->width || y >= buf->height || !estimator) + return RES_BAD_ARG; + *estimator = estimator_buffer_grab(buf, x, y); + return RES_OK; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +estimator_buffer_create + (struct sdis_device* dev, + const size_t width, + const size_t height, + struct sdis_estimator_buffer** out_buf) +{ + struct sdis_estimator_buffer* buf = NULL; + size_t i; + res_T res = RES_OK; + + if(!dev || !width || !height || !out_buf) { + res = RES_BAD_ARG; + goto error; + } + + buf = MEM_CALLOC(dev->allocator, 1, sizeof(*buf)); + if(!buf) { res = RES_MEM_ERR; goto error; } + + ref_init(&buf->ref); + SDIS(device_ref_get(dev)); + buf->dev = dev; + buf->width = width; + buf->height = height; + + buf->estimators = MEM_CALLOC + (dev->allocator, width*height, sizeof(*buf->estimators)); + if(!buf->estimators) { + log_err(dev, "%s: could not allocate a buffer of estimators of %lux%lu.\n", + FUNC_NAME, (unsigned long)width, (unsigned long)height); + res = RES_MEM_ERR; + goto error; + } + + FOR_EACH(i, 0, width*height) { + res = estimator_create(dev, SDIS_ESTIMATOR_TEMPERATURE, buf->estimators+i); + if(res != RES_OK) goto error; + } + +exit: + if(out_buf) *out_buf = buf; + return res; +error: + if(buf) { + SDIS(estimator_buffer_ref_put(buf)); + buf = NULL; + } + goto exit; +} + +struct sdis_estimator* +estimator_buffer_grab + (const struct sdis_estimator_buffer* buf, + const size_t x, + const size_t y) +{ + ASSERT(x < buf->width && y < buf->height); + return buf->estimators[y*buf->width + x]; +} + diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h @@ -0,0 +1,40 @@ +/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com) + * + * 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/>. */ + +#ifndef SDIS_ESTIMATOR_BUFFER_C_H +#define SDIS_ESTIMATOR_BUFFER_C_H + +#include <rsys/rsys.h> + +/* Forward declaration */ +struct sdis_device; +struct sdis_estimator; +struct sdis_estimator_buffer; + +extern LOCAL_SYM res_T +estimator_buffer_create + (struct sdis_device* dev, + const size_t width, + const size_t height, + struct sdis_estimator_buffer** buf); + +extern LOCAL_SYM struct sdis_estimator* +estimator_buffer_grab + (const struct sdis_estimator_buffer* buf, + const size_t x, + const size_t y); + +#endif /* SDIS_ESTIMATOR_BUFFER_C_H */ + diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -17,6 +17,7 @@ #define SDIS_ESTIMATOR_C_H #include "sdis_heat_path.h" +#include "sdis_misc.h" #include <rsys/math.h> #include <rsys/ref_count.h> @@ -34,9 +35,9 @@ enum flux_name { }; struct sdis_estimator { - struct sdis_mc temperature; - struct sdis_mc realisation_time; - struct sdis_mc fluxes[FLUX_NAMES_COUNT__]; + struct accum temperature; + struct accum realisation_time; + struct accum fluxes[FLUX_NAMES_COUNT__]; size_t nrealisations; size_t nfailures; @@ -77,13 +78,6 @@ estimator_setup_realisations_count estimator->nfailures = nrealisations - nsuccesses; } -#define SETUP_MC(McName, Sum, Sum2, Count) { \ - (McName).E = (Sum) / (double)(Count); \ - (McName).V = (Sum2) / (double)(Count) - (McName).E*(McName).E; \ - (McName).V = MMAX((McName).V, 0); \ - (McName).SE = sqrt((McName).V / (double)(Count)); \ -} (void)0 - static INLINE void estimator_setup_temperature (struct sdis_estimator* estim, @@ -91,7 +85,9 @@ estimator_setup_temperature const double sum2) { ASSERT(estim && estim->nrealisations); - SETUP_MC(estim->temperature, sum, sum2, estim->nrealisations); + estim->temperature.sum = sum; + estim->temperature.sum2 = sum2; + estim->temperature.count = estim->nrealisations; } static INLINE void @@ -101,7 +97,9 @@ estimator_setup_realisation_time const double sum2) { ASSERT(estim && estim->nrealisations); - SETUP_MC(estim->realisation_time, sum, sum2, estim->nrealisations); + estim->realisation_time.sum = sum; + estim->realisation_time.sum2 = sum2; + estim->realisation_time.count = estim->nrealisations; } static INLINE void @@ -112,10 +110,10 @@ estimator_setup_flux const double sum2) { ASSERT(estim && (unsigned)name < FLUX_NAMES_COUNT__ && estim->nrealisations); - SETUP_MC(estim->fluxes[name], sum, sum2, estim->nrealisations); + estim->fluxes[name].sum = sum; + estim->fluxes[name].sum2 = sum2; + estim->fluxes[name].count = estim->nrealisations; } -#undef SETUP_MC - #endif /* SDIS_PROBE_ESTIMATOR_C_H */ diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -17,6 +17,7 @@ #include "sdis_camera.h" #include "sdis_device_c.h" #include "sdis_estimator_c.h" +#include "sdis_estimator_buffer_c.h" #include "sdis_interface_c.h" #include <star/ssp.h> @@ -74,19 +75,15 @@ solve_pixel const size_t nrealisations, const int register_paths, /* Combination of enum sdis_heat_path_flag */ const double pix_sz[2], /* Pixel size in the normalized image plane */ - struct accum* acc_temp, - struct accum* acc_time, - struct sdis_accum* accum, struct sdis_estimator* estimator) { - double sum_weights = 0; - double sum_weights_sqr = 0; - size_t N = 0; /* #realisations that do not fail */ + struct accum acc_temp = ACCUM_NULL; + struct accum acc_time = ACCUM_NULL; size_t irealisation; res_T res = RES_OK; ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0); ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0); - ASSERT(acc_time && acc_temp && accum && estimator && time_range); + ASSERT(estimator && time_range); FOR_EACH(irealisation, 0, nrealisations) { struct time t0, t1; @@ -146,21 +143,16 @@ solve_pixel if(res_simul == RES_OK) { /* Update global accumulators */ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; - acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count; - acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count; - - /* Update per pixel accumulator */ - sum_weights += w; - sum_weights_sqr += w*w; - ++N; + acc_temp.sum += w; acc_temp.sum2 += w*w; ++acc_temp.count; + acc_time.sum += usec; acc_time.sum2 += usec*usec; ++acc_time.count; } } - /* Setup per pixel accumulator */ - accum->sum_weights = sum_weights; - accum->sum_weights_sqr = sum_weights_sqr; - accum->nweights = N; - accum->nfailures = nrealisations - N; + /* Setup the pixel estimator */ + ASSERT(acc_temp.count == acc_time.count); + estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); + estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); + estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); exit: return res; @@ -183,16 +175,13 @@ solve_tile const size_t spp, /* #samples per pixel */ const int register_paths, /* Combination of enum sdis_heat_path_flag */ const double pix_sz[2], /* Pixel size in the normalized image plane */ - struct accum* acc_temp, - struct accum* acc_time, - struct sdis_accum* accums, - struct sdis_estimator* estimator) + struct sdis_estimator_buffer* buf) { size_t mcode; /* Morton code of the tile pixel */ size_t npixels; res_T res = RES_OK; - ASSERT(scn && rng && mdm && cam && spp && origin && accums && Tref >= 0); - ASSERT(size &&size[0] && size[1] && acc_temp && acc_time && estimator); + ASSERT(scn && rng && mdm && cam && spp && origin && Tref >= 0); + ASSERT(size &&size[0] && size[1] && buf); ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range); /* Adjust the #pixels to process them wrt a morton order */ @@ -201,19 +190,21 @@ solve_tile FOR_EACH(mcode, 0, npixels) { size_t ipix[2]; - struct sdis_accum* accum; + struct sdis_estimator* estimator; ipix[0] = morton2D_decode((uint32_t)(mcode>>0)); if(ipix[0] >= size[0]) continue; ipix[1] = morton2D_decode((uint32_t)(mcode>>1)); if(ipix[1] >= size[1]) continue; - accum = accums + ipix[1]*size[0] + ipix[0]; ipix[0] = ipix[0] + origin[0]; ipix[1] = ipix[1] + origin[1]; + /* Fetch the pixel estimator */ + estimator = estimator_buffer_grab(buf, ipix[0], ipix[1]); + res = solve_pixel(scn, rng, mdm, cam, time_range, fp_to_meter, Tarad, Tref, - ipix, spp, register_paths, pix_sz, acc_temp, acc_time, accum, estimator); + ipix, spp, register_paths, pix_sz, estimator); if(res != RES_OK) goto error; } @@ -418,18 +409,13 @@ sdis_solve_camera const size_t height, /* #pixels in Y */ const size_t spp, /* #samples per pixel */ const int register_paths, /* Combination of enum sdis_heat_path_flag */ - sdis_write_accums_T writer, - void* writer_data, - struct sdis_estimator** out_estimator) + struct sdis_estimator_buffer** out_buf) { #define TILE_SIZE 32 /* definition in X & Y of a tile */ STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2); - struct sdis_estimator* estimator = NULL; + struct sdis_estimator_buffer* buf = NULL; struct sdis_medium* medium = NULL; - struct darray_accum* tiles = NULL; - struct accum* acc_temps = NULL; - struct accum* acc_times = NULL; struct ssp_rng_proxy* rng_proxy = NULL; struct ssp_rng** rngs = NULL; size_t ntiles_x, ntiles_y, ntiles; @@ -438,8 +424,8 @@ sdis_solve_camera size_t i; ATOMIC res = RES_OK; - if(!scn || !cam || fp_to_meter <= 0 || Tref < 0 || !width - || !height || !spp || !writer || !out_estimator) { + if(!scn || !cam || fp_to_meter <= 0 || Tref < 0 || !width || !height || !spp + || !out_buf) { res = RES_BAD_ARG; goto error; } @@ -481,23 +467,6 @@ sdis_solve_camera if(res != RES_OK) goto error; } - /* Create the per thread accumulators */ - acc_temps = MEM_CALLOC - (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps)); - if(!acc_temps) { res = RES_MEM_ERR; goto error; } - acc_times = MEM_CALLOC - (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times)); - if(!acc_times) { res = RES_MEM_ERR; goto error; } - - /* Allocate per thread buffer of accumulations */ - tiles = darray_tile_data_get(&scn->dev->tiles); - ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads); - FOR_EACH(i, 0, scn->dev->nthreads) { - const size_t naccums = TILE_SIZE * TILE_SIZE; - res = darray_accum_resize(tiles+i, naccums); - if(res != RES_OK) goto error; - } - ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y)); @@ -506,8 +475,8 @@ sdis_solve_camera pix_sz[0] = 1.0 / (double)width; pix_sz[1] = 1.0 / (double)height; - /* Create the estimator */ - res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator); + /* Create the global estimator */ + res = estimator_buffer_create(scn->dev, width, height, &buf); if(res != RES_OK) goto error; omp_set_num_threads((int)scn->dev->nthreads); @@ -516,10 +485,7 @@ sdis_solve_camera size_t tile_org[2] = {0, 0}; size_t tile_sz[2] = {0, 0}; const int ithread = omp_get_thread_num(); - struct sdis_accum* accums = NULL; struct ssp_rng* rng = rngs[ithread]; - struct accum* acc_temp = &acc_temps[ithread]; - struct accum* acc_time = &acc_times[ithread]; res_T res_local = RES_OK; if(ATOMIC_GET(&res) != RES_OK) continue; @@ -535,40 +501,13 @@ sdis_solve_camera tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]); tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]); - /* Fetch the accumulations buffer */ - accums = darray_accum_data_get(tiles+ithread); - /* Draw the tile */ res_local = solve_tile(scn, rng, medium, cam, time_range, fp_to_meter, - Tarad, Tref, tile_org, tile_sz, spp, register_paths, pix_sz, acc_temp, - acc_time, accums, estimator); + Tarad, Tref, tile_org, tile_sz, spp, register_paths, pix_sz, buf); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } - - /* Write the accumulations */ - res_local = writer(writer_data, tile_org, tile_sz, accums); - if(res_local != RES_OK) { - ATOMIC_SET(&res, res_local); - continue; - } - } - - /* Setup the estimated temperature and per realisation time for the whole - * image */ - if(out_estimator) { - const size_t nrealisations = width*height*spp; - struct accum acc_temp; - struct accum acc_time; - - sum_accums(acc_temps, scn->dev->nthreads, &acc_temp); - sum_accums(acc_times, scn->dev->nthreads, &acc_time); - ASSERT(acc_temp.count == acc_time.count); - - estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); - estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); - estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); } exit: @@ -578,10 +517,8 @@ exit: } MEM_RM(scn->dev->allocator, rngs); } - if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps); - if(acc_times) MEM_RM(scn->dev->allocator, acc_times); if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); - if(out_estimator) *out_estimator = estimator; + if(out_buf) *out_buf = buf; return (res_T)res; error: goto exit; diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -433,39 +433,59 @@ geometry_add_shape } static void -dump_image(const struct sdis_accum_buffer* buf) +dump_image(const struct sdis_estimator_buffer* buf) { - struct sdis_accum_buffer_layout layout = SDIS_ACCUM_BUFFER_LAYOUT_NULL; struct image img; double* temps = NULL; - const struct sdis_accum* accums = NULL; double Tmax = -DBL_MAX; double Tmin = DBL_MAX; double norm; + size_t definition[2]; size_t i, ix, iy; CHK(buf != NULL); - OK(sdis_accum_buffer_get_layout(buf, &layout)); + OK(sdis_estimator_buffer_get_definition(buf, definition)); - temps = mem_alloc(layout.width*layout.height*sizeof(double)); + temps = mem_alloc(definition[0]*definition[1]*sizeof(double)); CHK(temps != NULL); - OK(sdis_accum_buffer_map(buf, &accums)); - /* Check the results validity */ - FOR_EACH(i, 0, layout.height * layout.width) { - CHK(accums[i].nweights + accums[i].nfailures == SPP); - CHK(accums[i].nfailures <= SPP/100); - CHK(accums[i].sum_weights >= 0); - CHK(accums[i].sum_weights_sqr >= 0); + FOR_EACH(iy, 0, definition[1]) { + FOR_EACH(ix, 0, definition[0]) { + const struct sdis_estimator* estimator; + struct sdis_mc T; + struct sdis_mc time; + size_t nreals; + size_t nfails; + + BA(sdis_estimator_buffer_at(NULL, ix, iy, &estimator)); + BA(sdis_estimator_buffer_at(buf, definition[0]+1, iy, &estimator)); + BA(sdis_estimator_buffer_at(buf, ix, definition[1]+1, &estimator)); + BA(sdis_estimator_buffer_at(buf, ix, iy, NULL)); + OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + CHK(nreals + nfails == SPP); + CHK(T.E > 0); + CHK(time.E > 0); + } } - + /* Compute the per pixel temperature */ - FOR_EACH(iy, 0, layout.height) { - const struct sdis_accum* row_accums = accums + iy * layout.width; - double* row = temps + iy * layout.width; - FOR_EACH(ix, 0, layout.width) { - row[ix] = row_accums[ix].sum_weights / (double)row_accums[ix].nweights; + FOR_EACH(iy, 0, definition[1]) { + double* row = temps + iy * definition[0]; + FOR_EACH(ix, 0, definition[0]) { + const struct sdis_estimator* estimator; + struct sdis_mc T; + + OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + OK(sdis_estimator_get_temperature(estimator, &T)); + + row[ix] = T.E; Tmax = MMAX(row[ix], Tmax); Tmin = MMIN(row[ix], Tmin); } @@ -481,10 +501,10 @@ dump_image(const struct sdis_accum_buffer* buf) OK(image_init(NULL, &img)); OK(image_setup(&img, IMG_WIDTH, IMG_HEIGHT, IMG_WIDTH*3, IMAGE_RGB8, NULL)); - FOR_EACH(iy, 0, layout.height) { - const double* src_row = temps + iy*layout.width; + FOR_EACH(iy, 0, definition[1]) { + const double* src_row = temps + iy*definition[0]; char* dst_row = img.pixels + iy*img.pitch; - FOR_EACH(ix, 0, layout.width) { + FOR_EACH(ix, 0, definition[0]) { unsigned char* pixels = (unsigned char*) (dst_row + ix * sizeof_image_format(img.format)); const unsigned char T = (unsigned char) @@ -511,10 +531,9 @@ main(int argc, char** argv) struct s3dut_mesh_data msh_data; struct sdis_mc T = SDIS_MC_NULL; struct sdis_mc time = SDIS_MC_NULL; - struct sdis_accum_buffer* buf = NULL; struct sdis_camera* cam = NULL; struct sdis_device* dev = NULL; - struct sdis_estimator* estimator = NULL; + struct sdis_estimator_buffer* buf = NULL; struct sdis_medium* solid = NULL; struct sdis_medium* fluid0 = NULL; struct sdis_medium* fluid1 = NULL; @@ -526,6 +545,7 @@ main(int argc, char** argv) struct interf interface_param = INTERF_NULL; size_t ntris, npos; size_t nreals, nfails; + size_t definition[2]; double pos[3]; double tgt[3]; double up[3]; @@ -604,47 +624,42 @@ main(int argc, char** argv) OK(sdis_camera_set_fov(cam, MDEG2RAD(70))); OK(sdis_camera_look_at(cam, pos, tgt, up)); - /* Create the accum buffer */ - OK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf)); - #if 0 dump_mesh(stdout, geom.positions, npos, geom.indices, ntris); exit(0); #endif BA(sdis_solve_camera(NULL, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, NULL, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, NULL, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 0, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 1, 300, -1, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, 0, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, 0, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - 0, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, NULL, buf, &estimator)); + 0, SDIS_HEAT_PATH_NONE, &buf)); BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, NULL)); + SPP, SDIS_HEAT_PATH_NONE, NULL)); trange[0] = -1; BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); trange[0] = 10; trange[1] = 1; BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); trange[0] = trange[1] = INF; /* Launch the simulation */ OK(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, sdis_accum_buffer_write, buf, &estimator)); + SPP, SDIS_HEAT_PATH_NONE, &buf)); - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + /*OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); @@ -654,19 +669,23 @@ main(int argc, char** argv) fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE); fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE); fprintf(stderr, "#failures = %lu/%lu\n", - (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP)); + (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));*/ + BA(sdis_estimator_buffer_get_definition(NULL, definition)); + BA(sdis_estimator_buffer_get_definition(buf, NULL)); + OK(sdis_estimator_buffer_get_definition(buf, definition)); + CHK(definition[0] == IMG_WIDTH); + CHK(definition[1] == IMG_HEIGHT); /* Write the image */ dump_image(buf); /* Release memory */ - OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_buffer_ref_put(buf)); OK(sdis_scene_ref_put(scn)); OK(sdis_camera_ref_put(cam)); OK(sdis_interface_ref_put(interf0)); OK(sdis_interface_ref_put(interf1)); OK(sdis_device_ref_put(dev)); - OK(sdis_accum_buffer_ref_put(buf)); geometry_release(&geom); check_memory_allocator(&allocator);