atrstm.c (13313B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique 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 /* snprintf & lround support */ 18 19 #include "atrstm.h" 20 #include "atrstm_c.h" 21 #include "atrstm_cache.h" 22 #include "atrstm_log.h" 23 #include "atrstm_setup_octrees.h" 24 25 #include <astoria/atrtp.h> 26 27 #include <star/suvm.h> 28 #include <star/svx.h> 29 30 #include <rsys/cstr.h> 31 #include <rsys/double3.h> 32 33 #include <omp.h> 34 35 /******************************************************************************* 36 * Helper functions 37 ******************************************************************************/ 38 static int 39 check_args(const struct atrstm_args* args) 40 { 41 if(!args 42 || !args->sth_filename 43 || (args->spectral_type == ATRSTM_SPECTRAL_LW && !args->atrck_filename) 44 || !args->atrtp_filename 45 || !args->atrri_filename 46 || !args->name 47 || (unsigned)args->spectral_type >= ATRSTM_SPECTRAL_TYPES_COUNT__ 48 || args->wlen_range[0] > args->wlen_range[1] 49 || args->fractal_prefactor <= 0 50 || args->fractal_dimension <= 0 51 || args->optical_thickness < 0 52 || !args->nthreads) { 53 return 0; 54 } 55 56 if(args->auto_grid_definition) { 57 if(!args->auto_grid_definition_hint) 58 return 0; 59 } else { 60 if(!args->grid_max_definition[0] 61 || !args->grid_max_definition[1] 62 || !args->grid_max_definition[2]) 63 return 0; 64 } 65 66 return 1; 67 } 68 69 static FINLINE unsigned 70 round_pow2(const unsigned val) 71 { 72 const unsigned next_pow2 = (unsigned)round_up_pow2(val); 73 if(next_pow2 - val <= next_pow2/4) { 74 return next_pow2; 75 } else { 76 return next_pow2/2; 77 } 78 } 79 80 static void 81 compute_default_grid_definition 82 (const struct atrstm* atrstm, 83 const unsigned definition_hint, 84 unsigned def[3]) 85 { 86 double low[3]; 87 double upp[3]; 88 double sz[3]; 89 double sz_max; 90 double vxsz; 91 int iaxis_max; 92 int iaxis_remain[2]; 93 int i; 94 ASSERT(atrstm && def); 95 96 SUVM(volume_get_aabb(atrstm->volume, low, upp)); 97 sz[0] = upp[0] - low[0]; 98 sz[1] = upp[1] - low[1]; 99 sz[2] = upp[2] - low[2]; 100 101 /* Define the dimension along which the volume AABB is the greatest */ 102 sz_max = -DBL_MAX; 103 FOR_EACH(i, 0, 3) { 104 if(sz[i] > sz_max) { iaxis_max = i; sz_max = sz[i]; } 105 } 106 107 /* Define the other axes */ 108 iaxis_remain[0] = (iaxis_max + 1) % 3; 109 iaxis_remain[1] = (iaxis_max + 2) % 3; 110 111 /* Fix the definition to along the maximum axis and compute the voxel size 112 * along this axis */ 113 def[iaxis_max] = round_pow2(definition_hint); 114 vxsz = sz[iaxis_max] / def[iaxis_max]; 115 116 /* Compute the definition along the 2 remaining axes. First, compute them by 117 * assuming that the voxels are cubics and then round these definitions to 118 * their nearest power of two. */ 119 def[iaxis_remain[0]] = (unsigned)lround(sz[iaxis_remain[0]]/vxsz); 120 def[iaxis_remain[1]] = (unsigned)lround(sz[iaxis_remain[1]]/vxsz); 121 def[iaxis_remain[0]] = round_pow2(def[iaxis_remain[0]]); 122 def[iaxis_remain[1]] = round_pow2(def[iaxis_remain[1]]); 123 } 124 125 static void 126 release_atrstm(ref_T* ref) 127 { 128 struct atrstm* atrstm; 129 ASSERT(ref); 130 atrstm = CONTAINER_OF(ref, struct atrstm, ref); 131 octrees_clean(atrstm); 132 if(atrstm->atrtp) ATRTP(ref_put(atrstm->atrtp)); 133 if(atrstm->atrri) ATRRI(ref_put(atrstm->atrri)); 134 if(atrstm->suvm) SUVM(device_ref_put(atrstm->suvm)); 135 if(atrstm->volume) SUVM(volume_ref_put(atrstm->volume)); 136 if(atrstm->svx) SVX(device_ref_put(atrstm->svx)); 137 if(atrstm->cache) cache_ref_put(atrstm->cache); 138 if(atrstm->logger == &atrstm->logger__) logger_release(&atrstm->logger__); 139 str_release(&atrstm->name); 140 ASSERT(MEM_ALLOCATED_SIZE(&atrstm->svx_allocator) == 0); 141 mem_shutdown_proxy_allocator(&atrstm->svx_allocator); 142 MEM_RM(atrstm->allocator, atrstm); 143 } 144 145 /******************************************************************************* 146 * Exported functions 147 ******************************************************************************/ 148 res_T 149 atrstm_create 150 (struct logger* logger, /* NULL <=> use default logger */ 151 struct mem_allocator* mem_allocator, /* NULL <=> use default allocator */ 152 const struct atrstm_args* args, 153 struct atrstm** out_atrstm) 154 { 155 struct atrstm* atrstm = NULL; 156 struct mem_allocator* allocator = NULL; 157 int nthreads_max; 158 res_T res = RES_OK; 159 160 if(!out_atrstm || !check_args(args)) { 161 res = RES_BAD_ARG; 162 goto error; 163 } 164 165 allocator = mem_allocator ? mem_allocator : &mem_default_allocator; 166 atrstm = MEM_CALLOC(allocator, 1, sizeof(*atrstm)); 167 if(!atrstm) { 168 if(args->verbose) { 169 #define ERR_STR "Could not allocate the AtrSTM data structure.\n" 170 if(logger) { 171 logger_print(logger, LOG_ERROR, ERR_STR); 172 } else { 173 fprintf(stderr, MSG_ERROR_PREFIX ERR_STR); 174 } 175 #undef ERR_STR 176 } 177 res = RES_MEM_ERR; 178 goto error; 179 } 180 nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs()); 181 ref_init(&atrstm->ref); 182 atrstm->allocator = allocator; 183 atrstm->verbose = args->verbose; 184 atrstm->nthreads = MMIN(args->nthreads, (unsigned)nthreads_max); 185 atrstm->fractal_prefactor = args->fractal_prefactor; 186 atrstm->fractal_dimension = args->fractal_dimension; 187 atrstm->spectral_type = args->spectral_type; 188 atrstm->wlen_range[0] = args->wlen_range[0]; 189 atrstm->wlen_range[1] = args->wlen_range[1]; 190 atrstm->optical_thickness = args->optical_thickness; 191 192 str_init(allocator, &atrstm->name); 193 194 if(logger) { 195 atrstm->logger = logger; 196 } else { 197 setup_log_default(atrstm); 198 } 199 200 if(args->use_simd) { 201 #ifdef ATRSTM_USE_SIMD 202 atrstm->use_simd = 1; 203 #else 204 log_warn(atrstm, 205 "AtrSTM cannot use the SIMD instruction set: " 206 "it was compiled without SIMD support.\n"); 207 atrstm->use_simd = 0; 208 #endif 209 } 210 211 if(args->spectral_type == ATRSTM_SPECTRAL_SW 212 && args->wlen_range[0] != args->wlen_range[1]) { 213 log_err(atrstm, 214 "Invalid spectral range [%g, %g], only monochromatic computations are " 215 "supported in shortwave.\n", 216 args->wlen_range[0], 217 args->wlen_range[1]); 218 res = RES_BAD_ARG; 219 goto error; 220 } 221 222 res = str_set(&atrstm->name, args->name); 223 if(res != RES_OK) { 224 log_err(atrstm, "Cannot setup the gas mixture name to `%s' -- %s.\n", 225 args->name, res_to_cstr(res)); 226 goto error; 227 } 228 229 /* Setup the allocator used by the SVX library */ 230 res = mem_init_proxy_allocator(&atrstm->svx_allocator, atrstm->allocator); 231 if(res != RES_OK) { 232 log_err(atrstm, 233 "Cannot initialise the allocator used to manager the Star-VoXel memory " 234 "-- %s.\n", res_to_cstr(res)); 235 goto error; 236 } 237 238 /* Load the refractive index */ 239 res = atrri_create 240 (atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->atrri); 241 if(res != RES_OK) goto error; 242 res = atrri_load(atrstm->atrri, args->atrri_filename); 243 if(res != RES_OK) goto error; 244 245 /* In shortwave, pre-fetched the refractive index */ 246 if(atrstm->spectral_type != ATRSTM_SPECTRAL_SW) { 247 atrstm->refract_id = ATRRI_REFRACTIVE_INDEX_NULL; 248 } else { 249 const double wlen = atrstm->wlen_range[0]; 250 ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]); 251 res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &atrstm->refract_id); 252 if(res != RES_OK) { 253 log_err(atrstm, 254 "Could not fetch the refractive index of the shortwave wavelength %g " 255 "-- %s.\n", wlen, res_to_cstr(res)); 256 goto error; 257 } 258 } 259 260 /* Create the Star-UnstructuredVolumetricMesh device */ 261 res = suvm_device_create 262 (atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->suvm); 263 if(res != RES_OK) goto error; 264 265 /* Structure the volumetric mesh */ 266 res = setup_unstructured_volumetric_mesh 267 (atrstm, args->precompute_normals, args->sth_filename, &atrstm->volume); 268 if(res != RES_OK) goto error; 269 270 /* Define the grid definition */ 271 if(args->auto_grid_definition) { 272 compute_default_grid_definition 273 (atrstm, args->auto_grid_definition_hint, atrstm->grid_max_definition); 274 } else { 275 atrstm->grid_max_definition[0] = args->grid_max_definition[0]; 276 atrstm->grid_max_definition[1] = args->grid_max_definition[1]; 277 atrstm->grid_max_definition[2] = args->grid_max_definition[2]; 278 } 279 280 /* Load the thermodynamic properties of the volumetric mesh */ 281 res = atrtp_create 282 (atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->atrtp); 283 if(res != RES_OK) goto error; 284 res = atrtp_load(atrstm->atrtp, args->atrtp_filename); 285 if(res != RES_OK) goto error; 286 287 /* Create the Star-VoXel device */ 288 res = svx_device_create 289 (atrstm->logger, &atrstm->svx_allocator, atrstm->verbose, &atrstm->svx); 290 if(res != RES_OK) goto error; 291 292 /* Create the cache if necessary */ 293 if(args->cache_filename) { 294 res = cache_create(atrstm, args->cache_filename, &atrstm->cache); 295 if(res != RES_OK) goto error; 296 } 297 298 /* Build the octree */ 299 res = setup_octrees(atrstm); 300 if(res != RES_OK) goto error; 301 302 exit: 303 if(out_atrstm) *out_atrstm = atrstm; 304 return res; 305 error: 306 if(atrstm) { ATRSTM(ref_put(atrstm)); atrstm = NULL; } 307 goto exit; 308 } 309 310 res_T 311 atrstm_ref_get(struct atrstm* atrstm) 312 { 313 if(!atrstm) return RES_BAD_ARG; 314 ref_get(&atrstm->ref); 315 return RES_OK; 316 } 317 318 res_T 319 atrstm_ref_put(struct atrstm* atrstm) 320 { 321 if(!atrstm) return RES_BAD_ARG; 322 ref_put(&atrstm->ref, release_atrstm); 323 return RES_OK; 324 } 325 326 const char* 327 atrstm_get_name(const struct atrstm* atrstm) 328 { 329 ASSERT(atrstm); 330 return str_cget(&atrstm->name); 331 } 332 333 void 334 atrstm_get_aabb(const struct atrstm* atrstm, double lower[3], double upper[3]) 335 { 336 struct svx_tree_desc tree_desc = SVX_TREE_DESC_NULL; 337 ASSERT(atrstm && lower && upper); 338 SVX(tree_get_desc(atrstm->octree, &tree_desc)); 339 d3_set(lower, tree_desc.lower); 340 d3_set(upper, tree_desc.upper); 341 } 342 343 res_T 344 atrstm_at 345 (const struct atrstm* atrstm, 346 const double pos[3], 347 struct suvm_primitive* prim, /* Volumetric primitive */ 348 double bcoords[4]) /* `pos' in `prim' */ 349 { 350 res_T res = RES_OK; 351 352 if(!atrstm || !pos || !prim || !bcoords) { 353 res = RES_BAD_ARG; 354 goto error; 355 } 356 357 /* Find the primitive into which pos lies */ 358 res = suvm_volume_at(atrstm->volume, pos, prim, bcoords); 359 if(res != RES_OK) { 360 log_err(atrstm, 361 "Error accessing the volumetric mesh at {%g, %g, %g} -- %s.\n", 362 SPLIT3(pos), res_to_cstr(res)); 363 goto error; 364 } 365 366 exit: 367 return res; 368 error: 369 goto exit; 370 } 371 372 /******************************************************************************* 373 * Local functions 374 ******************************************************************************/ 375 void 376 dump_memory_size 377 (const size_t size, /* In Bytes */ 378 size_t* real_dump_len, /* May be NULL */ 379 char* dump, /* May be NULL */ 380 size_t max_dump_len) 381 { 382 size_t available_dump_space = max_dump_len; 383 char* dst = dump; 384 const size_t KILO_BYTE = 1024; 385 const size_t MEGA_BYTE = 1024*KILO_BYTE; 386 const size_t GIGA_BYTE = 1024*MEGA_BYTE; 387 size_t ngigas, nmegas, nkilos; 388 size_t nbytes = size; 389 390 #define DUMP(Size, Suffix) \ 391 { \ 392 const int len = snprintf \ 393 (dst, available_dump_space, \ 394 "%li %s", (long)Size, Size > 1 ? Suffix "s ": Suffix " "); \ 395 ASSERT(len >= 0); \ 396 if(real_dump_len) { \ 397 *real_dump_len += (size_t)len; \ 398 } \ 399 if((size_t)len < available_dump_space) { \ 400 dst += len; \ 401 available_dump_space -= (size_t)len; \ 402 } else if(dst) { \ 403 available_dump_space = 0; \ 404 dst = NULL; \ 405 } \ 406 } (void) 0 407 408 if((ngigas = nbytes / GIGA_BYTE) != 0) { 409 DUMP(ngigas, "GB"); 410 nbytes -= ngigas * GIGA_BYTE; 411 } 412 if((nmegas = nbytes / MEGA_BYTE) != 0) { 413 DUMP(nmegas, "MB"); 414 nbytes -= nmegas * MEGA_BYTE; 415 } 416 if((nkilos = nbytes / KILO_BYTE) != 0) { 417 DUMP(nkilos, "kB"); 418 nbytes -= nkilos * KILO_BYTE; 419 } 420 if(nbytes) { 421 DUMP(nbytes, "Byte"); 422 } 423 424 /* Remove last space */ 425 if(real_dump_len) *real_dump_len -= 1; 426 if(dst != dump && dst[-1] == ' ') dst[-1] = '\0'; 427 428 #undef DUMP 429 }