sgf_estimator.c (11228B)
1 /* Copyright (C) 2021, 2024 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2015-2018 EDF S.A., France (syrthes-support@edf.fr) 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* nextafterf support */ 18 19 #include "sgf.h" 20 #include "sgf_device_c.h" 21 #include "sgf_realisation.h" 22 23 #include <star/s3d.h> 24 #include <star/ssp.h> 25 26 #include <rsys/dynamic_array.h> 27 #include <rsys/logger.h> 28 #include <rsys/mem_allocator.h> 29 #include <rsys/ref_count.h> 30 31 /* A random walk may fail du to numerical inaccuracy. The following constant 32 * empirically defines the maximum number of "random walk" attempts before an 33 * error occurs. */ 34 #define MAX_FAILURES 10 35 36 /* Generate the accum dynamic array data type */ 37 #define DARRAY_NAME accum 38 #define DARRAY_DATA struct accum 39 #include <rsys/dynamic_array.h> 40 41 /* Generate the 2D realisation function */ 42 #define SGF_DIMENSIONALITY 2 43 #include "sgf_realisation.h" 44 45 /* Generate the 3D realisation function */ 46 #define SGF_DIMENSIONALITY 3 47 #include "sgf_realisation.h" 48 49 /* Estimator of the Gebhart Factor of a given face to all the other ones */ 50 struct sgf_estimator { 51 /* Per spectral band & per primitive radiative flux */ 52 struct darray_accum buf; 53 /* Per spectral band radiative flux going to the infinity */ 54 struct darray_accum buf_infinity; 55 /* Per spectral band radiative flux absorbed by the medium */ 56 struct darray_accum buf_medium; 57 58 size_t nsteps; /* # realisations */ 59 size_t nprims; /* # primitives */ 60 size_t nbands; /* # spectral bands */ 61 62 struct sgf_device* dev; 63 ref_T ref; 64 }; 65 66 /******************************************************************************* 67 * Helper functions 68 ******************************************************************************/ 69 static FINLINE void 70 setup_status 71 (const struct accum* acc, 72 const size_t nsteps, 73 struct sgf_status* status) 74 { 75 ASSERT(acc && status && nsteps); 76 77 status->E = acc->radiative_flux / (double)nsteps; /* Expected value */ 78 status->V = acc->sqr_radiative_flux / (double)nsteps - status->E*status->E; 79 status->V = MMAX(status->V, 0); 80 status->SE = sqrt(status->V / (double)nsteps); /* Standard error */ 81 status->nsteps = nsteps; /* # realisations */ 82 } 83 84 static res_T 85 estimator_create(struct sgf_device* dev, struct sgf_estimator** out_estimator) 86 { 87 struct sgf_estimator* estimator = NULL; 88 res_T res = RES_OK; 89 ASSERT(dev && out_estimator); 90 91 estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sgf_estimator)); 92 if(!estimator) { 93 log_error(dev, 94 "Not enough memory space: couldn't allocate the Gebhart Factor estimator.\n"); 95 res = RES_MEM_ERR; 96 goto error; 97 } 98 ref_init(&estimator->ref); 99 SGF(device_ref_get(dev)); 100 estimator->dev = dev; 101 darray_accum_init(dev->allocator, &estimator->buf); 102 darray_accum_init(dev->allocator, &estimator->buf_infinity); 103 darray_accum_init(dev->allocator, &estimator->buf_medium); 104 105 exit: 106 *out_estimator = estimator; 107 return res; 108 109 error: 110 if(estimator) { 111 SGF(estimator_ref_put(estimator)); 112 estimator = NULL; 113 } 114 goto exit; 115 } 116 117 static void 118 estimator_release(ref_T* ref) 119 { 120 struct sgf_device* dev; 121 struct sgf_estimator* estimator; 122 ASSERT(ref); 123 estimator = CONTAINER_OF(ref, struct sgf_estimator, ref); 124 darray_accum_release(&estimator->buf); 125 darray_accum_release(&estimator->buf_infinity); 126 darray_accum_release(&estimator->buf_medium); 127 dev = estimator->dev; 128 MEM_RM(dev->allocator, estimator); 129 SGF(device_ref_put(dev)); 130 } 131 132 /******************************************************************************* 133 * Exported functions 134 ******************************************************************************/ 135 res_T 136 sgf_integrate 137 (struct sgf_scene* scn, 138 const size_t iprim, 139 struct ssp_rng* rng, 140 const size_t steps_count, 141 struct sgf_estimator** out_estimator) 142 { 143 struct htable_bounce path; 144 struct sgf_estimator* estimator = NULL; 145 size_t istep; 146 size_t iband; 147 size_t nprims; 148 void* view; 149 res_T res = RES_OK; 150 gebhart_radiative_path_T gebhart_radiative_path; 151 152 if(!scn) return RES_BAD_ARG; 153 htable_bounce_init(scn->dev->allocator, &path); 154 155 if(!steps_count || !rng || !scn || !out_estimator) { 156 res = RES_BAD_ARG; 157 goto error; 158 } 159 160 res = estimator_create(scn->dev, &estimator); 161 if(res != RES_OK) goto error; 162 163 switch(scn->dimensionality) { 164 case 2: 165 view = scn->geometry.s2d.view; 166 gebhart_radiative_path = gebhart_radiative_path_2d; 167 if(view) S2D(scene_view_primitives_count(view, &nprims)); 168 break; 169 case 3: 170 view = scn->geometry.s3d.view; 171 gebhart_radiative_path = gebhart_radiative_path_3d; 172 if(view) S3D(scene_view_primitives_count(view, &nprims)); 173 break; 174 default: FATAL("Unreachable code\n"); break; 175 } 176 177 /* Check scene active sessions */ 178 if(!view) { 179 log_error(scn->dev, 180 "%s: no active integration on subimitted scene.\n", FUNC_NAME); 181 res = RES_BAD_OP; 182 goto error; 183 } 184 185 /* Check submitted primitive_id */ 186 if(iprim >= nprims) { 187 log_error(scn->dev, "%s: invalid primitive index `%lu'\n", 188 FUNC_NAME, (unsigned long)iprim); 189 res = RES_BAD_ARG; 190 goto error; 191 } 192 193 /* Allocate and init the accumulators of radiative flux */ 194 res = darray_accum_resize(&estimator->buf, nprims*scn->nbands); 195 if(res != RES_OK) { 196 log_error(scn->dev, "%s: couldn't allocate the Gebhart Factor result buffer.\n", 197 FUNC_NAME); 198 goto error; 199 } 200 memset(darray_accum_data_get(&estimator->buf), 0, 201 darray_accum_size_get(&estimator->buf)*sizeof(struct accum)); 202 203 if(scn->has_medium) { 204 /* Allocate and init the accumulators of absorbed radiative flux */ 205 res = darray_accum_resize(&estimator->buf_medium, scn->nbands); 206 if(res != RES_OK) { 207 log_error(scn->dev, 208 "%s: couldn't allocate the accumulators of the per spectral radiative flux " 209 "absorbed by the medium.\n", FUNC_NAME); 210 goto error; 211 } 212 memset(darray_accum_data_get(&estimator->buf_medium), 0, 213 darray_accum_size_get(&estimator->buf_medium)*sizeof(struct accum)); 214 } else { 215 /* Allocate and init the accumulators of infinite radiative flux */ 216 res = darray_accum_resize(&estimator->buf_infinity, scn->nbands); 217 if(res != RES_OK) { 218 log_error(scn->dev, 219 "%s: couldn't allocate the accumulators of the per spectral band radiative flux" 220 " that goes to the infinite.\n", FUNC_NAME); 221 goto error; 222 } 223 memset(darray_accum_data_get(&estimator->buf_infinity), 0, 224 darray_accum_size_get(&estimator->buf_infinity)*sizeof(struct accum)); 225 } 226 227 /* Invoke the Gebhart factor integration. */ 228 FOR_EACH(iband, 0, scn->nbands) { 229 struct accum* accums = NULL; 230 struct accum* accum_infinity = NULL; 231 struct accum* accum_medium = NULL; 232 double ka = -1; /* Absorption coefficient of the medium */ 233 234 accums = darray_accum_data_get(&estimator->buf) + iband * nprims; 235 if(!scn->has_medium) { 236 accum_infinity = darray_accum_data_get(&estimator->buf_infinity) + iband; 237 } else { 238 accum_medium = darray_accum_data_get(&estimator->buf_medium) + iband; 239 ka = scene_get_absorption(scn, iprim, iband); 240 if(ka < 0) { 241 log_error(scn->dev, "%s: invalid absorption coefficient `%g'.\n", 242 FUNC_NAME, ka); 243 res = RES_BAD_ARG; 244 goto error; 245 } 246 } 247 248 FOR_EACH(istep, 0, steps_count) { 249 size_t nfailures = 0; 250 do { 251 res = gebhart_radiative_path(scn->dev, accums, accum_infinity, 252 accum_medium, rng, &path, ka, iband, iprim, scn); 253 if(res == RES_BAD_OP) { 254 log_error(scn->dev, "%s: reject radiative random walk.\n", FUNC_NAME); 255 ++nfailures; 256 } else if(res != RES_OK) { 257 goto error; 258 } 259 } while(res != RES_OK && nfailures < MAX_FAILURES); 260 261 if(++nfailures > MAX_FAILURES) { 262 log_error(scn->dev, "%s: too many numerical issues.\n", FUNC_NAME); 263 goto error; 264 } 265 } 266 } 267 estimator->nsteps = steps_count; 268 estimator->nprims = nprims; 269 estimator->nbands = scn->nbands; 270 271 exit: 272 if(out_estimator) *out_estimator = estimator; 273 htable_bounce_release(&path); 274 return res; 275 error: 276 if(estimator) SGF(estimator_ref_put(estimator)); 277 goto exit; 278 } 279 280 res_T 281 sgf_estimator_ref_get(struct sgf_estimator* estimator) 282 { 283 if(!estimator) return RES_BAD_ARG; 284 ref_get(&estimator->ref); 285 return RES_OK; 286 } 287 288 res_T 289 sgf_estimator_ref_put(struct sgf_estimator* estimator) 290 { 291 if(!estimator) return RES_BAD_ARG; 292 ref_put(&estimator->ref, estimator_release); 293 return RES_OK; 294 } 295 296 res_T 297 sgf_estimator_get_status 298 (struct sgf_estimator* estimator, 299 const size_t iprim, 300 const size_t iband, 301 struct sgf_status* status) 302 { 303 const struct accum* acc; 304 size_t iacc; 305 306 if(!estimator || !status) 307 return RES_BAD_ARG; 308 309 if(iprim >= estimator->nprims) { 310 log_error(estimator->dev, "%s: out of bound primitive index `%lu'.\n", 311 FUNC_NAME, (unsigned long)iprim); 312 return RES_BAD_ARG; 313 } 314 315 if(iband >= estimator->nbands) { 316 log_error(estimator->dev, "%s: out of bound spectral band index `%lu'.\n", 317 FUNC_NAME, (unsigned long)iband); 318 return RES_BAD_ARG; 319 } 320 321 iacc = iband * estimator->nprims + iprim; 322 acc = darray_accum_cdata_get(&estimator->buf) + iacc; 323 setup_status(acc, estimator->nsteps, status); 324 return RES_OK; 325 } 326 327 res_T 328 sgf_estimator_get_status_infinity 329 (struct sgf_estimator* estimator, 330 const size_t iband, 331 struct sgf_status* status) 332 { 333 const struct accum* acc; 334 335 if(!estimator || !status) 336 return RES_BAD_ARG; 337 338 if(iband >= estimator->nbands) { 339 log_error(estimator->dev, "%s: out of bound spectral band index `%lu'.\n", 340 FUNC_NAME, (unsigned long)iband); 341 return RES_BAD_ARG; 342 } 343 if(darray_accum_size_get(&estimator->buf_infinity) != 0) { 344 acc = darray_accum_cdata_get(&estimator->buf_infinity) + iband; 345 setup_status(acc, estimator->nsteps, status); 346 } else { 347 status->E = status->V = status->SE = 0; 348 status->nsteps = estimator->nsteps; 349 } 350 return RES_OK; 351 } 352 353 res_T 354 sgf_estimator_get_status_medium 355 (struct sgf_estimator* estimator, 356 const size_t iband, 357 struct sgf_status* status) 358 { 359 const struct accum* acc; 360 361 if(!estimator || !status) 362 return RES_BAD_ARG; 363 364 if(iband >= estimator->nbands) { 365 log_error(estimator->dev, "%s: out of bound spectral band index `%lu'\n.", 366 FUNC_NAME, (unsigned long)iband); 367 return RES_BAD_ARG; 368 } 369 if(darray_accum_size_get(&estimator->buf_medium) != 0) { 370 acc = darray_accum_cdata_get(&estimator->buf_medium) + iband; 371 setup_status(acc, estimator->nsteps, status); 372 } else { 373 status->E = status->V = status->SE = 0; 374 status->nsteps = estimator->nsteps; 375 } 376 return RES_OK; 377 } 378