htsky_cloud.c (37493B)
1 /* Copyright (C) 2018, 2019, 2020, 2021 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018, 2019 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018, 2019 Université Paul Sabatier 4 * 5 * This program is free software: you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation, either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 17 18 #define _POSIX_C_SOURCE 200809L /* nextafterf */ 19 20 #include "htsky_c.h" 21 #include "htsky_cloud.h" 22 #include "htsky_log.h" 23 #include "htsky_svx.h" 24 25 #include <high_tune/htgop.h> 26 #include <rsys/cstr.h> 27 #include <rsys/dynamic_array.h> 28 #include <star/svx.h> 29 30 #include <omp.h> 31 32 struct spectral_data { 33 size_t iband; /* Index of the spectral band */ 34 size_t iquad; /* Quadrature point into the band */ 35 }; 36 37 /* Define the dynamic array of spectral data */ 38 #define DARRAY_NAME specdata 39 #define DARRAY_DATA struct spectral_data 40 #include <rsys/dynamic_array.h> 41 42 /******************************************************************************* 43 * Helper functions 44 ******************************************************************************/ 45 static INLINE int 46 aabb_intersect 47 (const double aabb0_low[3], 48 const double aabb0_upp[3], 49 const double aabb1_low[3], 50 const double aabb1_upp[3]) 51 { 52 ASSERT(aabb0_low[0] < aabb0_upp[0] && aabb1_low[0] < aabb1_upp[0]); 53 ASSERT(aabb0_low[1] < aabb0_upp[1] && aabb1_low[1] < aabb1_upp[1]); 54 ASSERT(aabb0_low[2] < aabb0_upp[2] && aabb1_low[2] < aabb1_upp[2]); 55 return !(aabb0_upp[0] < aabb1_low[0]) && !(aabb0_low[0] > aabb1_upp[0]) 56 && !(aabb0_upp[1] < aabb1_low[1]) && !(aabb0_low[1] > aabb1_upp[1]) 57 && !(aabb0_upp[2] < aabb1_low[2]) && !(aabb0_low[2] > aabb1_upp[2]); 58 } 59 60 static res_T 61 setup_svx2htcp_z_lut(struct htsky* sky, double* cloud_upp_z) 62 { 63 double org, height; 64 double len_z; 65 size_t nsplits; 66 size_t iz, iz2; 67 res_T res = RES_OK; 68 ASSERT(sky && sky->htcp_desc.irregular_z && cloud_upp_z); 69 70 /* Find the min voxel size along Z and compute the length of a Z column */ 71 len_z = 0; 72 sky->lut_cell_sz = DBL_MAX; 73 FOR_EACH(iz, 0, sky->htcp_desc.spatial_definition[2]) { 74 len_z += sky->htcp_desc.vxsz_z[iz]; 75 sky->lut_cell_sz = MMIN(sky->lut_cell_sz, sky->htcp_desc.vxsz_z[iz]); 76 } 77 78 /* Allocate the svx2htcp LUT. This LUT is a regular table whose absolute 79 * size is greater or equal to a Z column in the htcp file. The size of its 80 * cells is the minimal voxel size in Z of the htcp file */ 81 nsplits = (size_t)ceil(len_z / sky->lut_cell_sz); 82 res = darray_split_resize(&sky->svx2htcp_z, nsplits); 83 if(res != RES_OK) { 84 log_err(sky, 85 "Could not allocate the table mapping regular to irregular Z.\n"); 86 goto error; 87 } 88 89 /* Setup the svx2htcp LUT. Each LUT entry stores the index of the current Z 90 * voxel in the htcp file that overlaps the entry lower bound as well as the 91 * lower bound in Z of the next htcp voxel. */ 92 iz2 = 0; 93 org = sky->htcp_desc.lower[2]; 94 height = org + sky->htcp_desc.vxsz_z[iz2]; 95 FOR_EACH(iz, 0, nsplits) { 96 const double upp_z = (double)(iz + 1) * sky->lut_cell_sz + org; 97 darray_split_data_get(&sky->svx2htcp_z)[iz].index = iz2; 98 darray_split_data_get(&sky->svx2htcp_z)[iz].height = height; 99 if(upp_z >= height && iz + 1 < nsplits) { 100 ASSERT(iz2 + 1 < sky->htcp_desc.spatial_definition[2]); 101 height += sky->htcp_desc.vxsz_z[++iz2]; 102 } 103 } 104 ASSERT(eq_eps(height - org, len_z, 1.e-6)); 105 106 exit: 107 *cloud_upp_z = height; 108 return res; 109 error: 110 darray_split_clear(&sky->svx2htcp_z); 111 height = -1; 112 goto exit; 113 } 114 115 static INLINE void 116 compute_k_bounds_regular_z 117 (const struct htsky* sky, 118 const size_t iband, 119 const double vox_svx_low[3], /* Lower bound of the SVX voxel */ 120 const double vox_svx_upp[3], /* Upper bound of the SVX voxel */ 121 double ka_bounds[2], 122 double ks_bounds[2], 123 double kext_bounds[2]) 124 { 125 size_t ivox[3]; 126 size_t igrid_low[3], igrid_upp[3]; 127 size_t i; 128 double ka, ks, kext; 129 double low[3], upp[3]; 130 double ipart; 131 double rho_da; /* Dry air density */ 132 double rct; 133 double Ca_avg; 134 double Cs_avg; 135 ASSERT(!sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp); 136 ASSERT(ka_bounds && ks_bounds && kext_bounds); 137 138 i = iband - sky->bands_range[0]; 139 140 /* Fetch the optical properties of the spectral band */ 141 Ca_avg = sky->bands[i].Ca_avg; 142 Cs_avg = sky->bands[i].Cs_avg; 143 144 /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by 145 * the SVX voxel */ 146 low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x; 147 low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x; 148 low[2] = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0]; 149 upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y; 150 upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y; 151 upp[2] = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0]; 152 igrid_low[0] = (size_t)low[0]; 153 igrid_low[1] = (size_t)low[1]; 154 igrid_low[2] = (size_t)low[2]; 155 igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0); 156 igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0); 157 igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0); 158 ASSERT(igrid_low[0] <= igrid_upp[0]); 159 ASSERT(igrid_low[1] <= igrid_upp[1]); 160 ASSERT(igrid_low[2] <= igrid_upp[2]); 161 162 /* Prepare the iteration over the grid voxels overlapped by the SVX voxel */ 163 ka_bounds[0] = ks_bounds[0] = kext_bounds[0] = DBL_MAX; 164 ka_bounds[1] = ks_bounds[1] = kext_bounds[1] =-DBL_MAX; 165 166 /* Iterate over the grid voxels overlapped by the SVX voxel */ 167 FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) { 168 FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) { 169 FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) { 170 /* Compute the radiative properties */ 171 rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox); 172 rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0); 173 ka = Ca_avg * rho_da * rct; 174 ks = Cs_avg * rho_da * rct; 175 kext = ka + ks; 176 /* Update the boundaries of the radiative properties */ 177 ka_bounds[0] = MMIN(ka_bounds[0], ka); 178 ka_bounds[1] = MMAX(ka_bounds[1], ka); 179 ks_bounds[0] = MMIN(ks_bounds[0], ks); 180 ks_bounds[1] = MMAX(ks_bounds[1], ks); 181 kext_bounds[0] = MMIN(kext_bounds[0], kext); 182 kext_bounds[1] = MMAX(kext_bounds[1], kext); 183 #ifndef NDEBUG 184 { 185 double tmp_low[3], tmp_upp[3]; 186 htcp_desc_get_voxel_aabb 187 (&sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp); 188 ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_svx_low, vox_svx_upp)); 189 } 190 #endif 191 } 192 } 193 } 194 } 195 196 static void 197 compute_k_bounds_irregular_z 198 (const struct htsky* sky, 199 const size_t iband, 200 const double vox_svx_low[3], /* Lower bound of the SVX voxel */ 201 const double vox_svx_upp[3], /* Upper bound of the SVX voxel */ 202 double ka_bounds[2], 203 double ks_bounds[2], 204 double kext_bounds[2]) 205 { 206 size_t ivox[3]; 207 size_t igrid_low[3], igrid_upp[3]; 208 size_t i; 209 double ka, ks, kext; 210 double low[2], upp[2]; 211 double ipart; 212 double rho_da; /* Dry air density */ 213 double rct; 214 double Ca_avg; 215 double Cs_avg; 216 double pos_z; 217 double lut_low, lut_upp; 218 size_t ilut_low, ilut_upp; 219 size_t ilut; 220 size_t ivox_z_prev = SIZE_MAX; 221 ASSERT(sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp); 222 ASSERT(ka_bounds && ks_bounds && kext_bounds); 223 224 i = iband - sky->bands_range[0]; 225 226 /* Fetch the optical properties of the spectral band */ 227 Ca_avg = sky->bands[i].Ca_avg; 228 Cs_avg = sky->bands[i].Cs_avg; 229 230 /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by 231 * the SVX voxel along the X and Y axes */ 232 low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x; 233 low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x; 234 upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y; 235 upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y; 236 igrid_low[0] = (size_t)low[0]; 237 igrid_low[1] = (size_t)low[1]; 238 igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0); 239 igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0); 240 ASSERT(igrid_low[0] <= igrid_upp[0]); 241 ASSERT(igrid_low[1] <= igrid_upp[1]); 242 243 /* Compute the *inclusive* bounds of the indices of the LUT cells 244 * overlapped by the SVX voxel */ 245 lut_low = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz; 246 lut_upp = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz; 247 ilut_low = (size_t)lut_low; 248 ilut_upp = (size_t)lut_upp; 249 250 /* A LUT coordinate equal to its search index means that this coordinate lies 251 * strictly on the boundary of a LUT cell. If this coordinate is the upper 252 * limit of a coordinate range, we subtract one from the LUT index to ensure 253 * that the corresponding cell is indeed included in the submitted range 254 * coordinates */ 255 if(lut_upp == ilut_upp) { 256 --ilut_upp; 257 } 258 259 ASSERT(ilut_low <= ilut_upp); 260 261 /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */ 262 ka_bounds[0] = ks_bounds[0] = kext_bounds[0] = DBL_MAX; 263 ka_bounds[1] = ks_bounds[1] = kext_bounds[1] =-DBL_MAX; 264 ivox_z_prev = SIZE_MAX; 265 pos_z = vox_svx_low[2]; 266 ASSERT(vox_svx_low[2] >= 267 sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_low)); 268 ASSERT(vox_svx_upp[2] <= 269 sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_upp+1)); 270 271 /* Iterate over the LUT cells overlapped by the voxel */ 272 FOR_EACH(ilut, ilut_low, ilut_upp+1) { 273 const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z)+ilut; 274 ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z)); 275 276 /* Compute the upper bound of the current LUT cell clamped to the voxel 277 * upper bound. */ 278 pos_z = MMIN((double)(ilut+1)*sky->lut_cell_sz, vox_svx_upp[2]); 279 280 igrid_low[2] = split->index; 281 igrid_upp[2] = split->index + (pos_z > split->height); 282 283 if(igrid_upp[2] >= sky->htcp_desc.spatial_definition[2] 284 && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */ 285 igrid_upp[2] = split->index; 286 } 287 288 FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) { 289 FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) { 290 FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) { 291 292 /* Does the LUT cell overlap an already handled LES voxel? */ 293 if(ivox[2] == ivox_z_prev) continue; 294 ivox_z_prev = ivox[2]; 295 296 /* Compute the radiative properties */ 297 rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox); 298 rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0); 299 ka = Ca_avg * rho_da * rct; 300 ks = Cs_avg * rho_da * rct; 301 kext = ka + ks; 302 303 /* Update the boundaries of the radiative properties */ 304 ka_bounds[0] = MMIN(ka_bounds[0], ka); 305 ka_bounds[1] = MMAX(ka_bounds[1], ka); 306 ks_bounds[0] = MMIN(ks_bounds[0], ks); 307 ks_bounds[1] = MMAX(ks_bounds[1], ks); 308 kext_bounds[0] = MMIN(kext_bounds[0], kext); 309 kext_bounds[1] = MMAX(kext_bounds[1], kext); 310 } 311 } 312 } 313 } 314 } 315 316 static void 317 cloud_vox_get_particle 318 (const size_t xyz[3], 319 float vox[], 320 const struct build_tree_context* ctx) 321 { 322 double vox_low[3], vox_upp[3]; 323 double ka[2], ks[2], kext[2]; 324 ASSERT(xyz && vox && ctx); 325 326 /* Compute the AABB of the SVX voxel */ 327 vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + ctx->sky->htcp_desc.lower[0]; 328 vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + ctx->sky->htcp_desc.lower[1]; 329 vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + ctx->sky->htcp_desc.lower[2]; 330 vox_upp[0] = vox_low[0] + ctx->vxsz[0]; 331 vox_upp[1] = vox_low[1] + ctx->vxsz[1]; 332 vox_upp[2] = vox_low[2] + ctx->vxsz[2]; 333 334 if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */ 335 compute_k_bounds_regular_z 336 (ctx->sky, ctx->iband, vox_low, vox_upp, ka, ks, kext); 337 } else { 338 ASSERT(xyz[2] < darray_split_size_get(&ctx->sky->svx2htcp_z)); 339 compute_k_bounds_irregular_z 340 (ctx->sky, ctx->iband, vox_low, vox_upp, ka, ks, kext); 341 } 342 343 /* Ensure that the single precision bounds include their double precision 344 * version. */ 345 if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX); 346 if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX); 347 if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX); 348 if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX); 349 if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX); 350 if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX); 351 352 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka[0]); 353 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka[1]); 354 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks[0]); 355 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks[1]); 356 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext[0]); 357 vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext[1]); 358 } 359 360 static void 361 compute_xh2o_range_regular_z 362 (const struct htsky* sky, 363 const double vox_svx_low[3], /* Lower bound of the SVX voxel */ 364 const double vox_svx_upp[3], /* Upper bound of the SVX voxel */ 365 double x_h2o_range[2]) 366 { 367 size_t ivox[3]; 368 size_t igrid_low[3], igrid_upp[3]; 369 double low[3], upp[3]; 370 double ipart; 371 ASSERT(!sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp); 372 ASSERT(x_h2o_range); 373 374 /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by 375 * the SVX voxel in X and Y */ 376 low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x; 377 low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x; 378 low[2] = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0]; 379 upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y; 380 upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y; 381 upp[2] = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0]; 382 igrid_low[0] = (size_t)low[0]; 383 igrid_low[1] = (size_t)low[1]; 384 igrid_low[2] = (size_t)low[2]; 385 igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0); 386 igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0); 387 igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0); 388 ASSERT(igrid_low[0] <= igrid_upp[0]); 389 ASSERT(igrid_low[1] <= igrid_upp[1]); 390 ASSERT(igrid_low[2] <= igrid_upp[2]); 391 392 /* Prepare the iteration overt the grid voxels overlapped by the SVX voxel */ 393 x_h2o_range[0] = DBL_MAX; 394 x_h2o_range[1] =-DBL_MAX; 395 396 FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) { 397 FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) { 398 FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) { 399 double x_h2o; 400 401 /* Compute the xH2O for the current LES voxel */ 402 x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox); 403 404 /* Update the xH2O range of the SVX voxel */ 405 x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]); 406 x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]); 407 #ifndef NDEBUG 408 { 409 double tmp_low[3], tmp_upp[3]; 410 htcp_desc_get_voxel_aabb 411 (&sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp); 412 ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_svx_low, vox_svx_upp)); 413 } 414 #endif 415 } 416 } 417 } 418 } 419 420 static void 421 compute_xh2o_range_irregular_z 422 (const struct htsky* sky, 423 const double vox_svx_low[3], /* Lower bound of the SVX voxel */ 424 const double vox_svx_upp[3], /* Upper bound of the SVX voxel */ 425 double x_h2o_range[2]) 426 { 427 size_t ivox[3]; 428 size_t igrid_low[3], igrid_upp[3]; 429 size_t ilut_low, ilut_upp; 430 size_t ilut; 431 size_t ivox_z_prev; 432 double lut_low, lut_upp; 433 double low[2], upp[2]; 434 double pos_z; 435 double ipart; 436 ASSERT(sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp); 437 ASSERT(x_h2o_range); 438 439 /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by 440 * the SVX voxel in X and Y */ 441 low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x; 442 low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x; 443 upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y; 444 upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y; 445 igrid_low[0] = (size_t)low[0]; 446 igrid_low[1] = (size_t)low[1]; 447 igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0); 448 igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0); 449 ASSERT(igrid_low[0] <= igrid_upp[0]); 450 ASSERT(igrid_low[1] <= igrid_upp[1]); 451 452 /* Compute the *inclusive* bounds of the indices of the LUT cells 453 * overlapped by the SVX voxel */ 454 lut_low = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz; 455 lut_upp = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz; 456 ilut_low = (size_t)lut_low; 457 ilut_upp = (size_t)lut_upp; 458 459 /* A LUT coordinate equal to its search index means that this coordinate lies 460 * strictly on the boundary of a LUT cell. If this coordinate is the upper 461 * limit of a coordinate range, we subtract one from the LUT index to ensure 462 * that the corresponding cell is indeed included in the submitted range 463 * coordinates */ 464 if(lut_upp == ilut_upp) { 465 --ilut_upp; 466 } 467 468 ASSERT(ilut_low <= ilut_upp); 469 470 /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */ 471 x_h2o_range[0] = DBL_MAX; 472 x_h2o_range[1] =-DBL_MAX; 473 ivox_z_prev = SIZE_MAX; 474 ASSERT(vox_svx_low[2] >= 475 sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_low)); 476 ASSERT(vox_svx_upp[2] <= 477 sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_upp+1)); 478 479 /* Iterate over the LUT cells overlapped by the voxel */ 480 FOR_EACH(ilut, ilut_low, ilut_upp+1) { 481 const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z) + ilut; 482 ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z)); 483 484 /* Compute the upper bound of the current LUT cell clamped to the voxel 485 * upper bound.*/ 486 pos_z = MMIN((double)(ilut+1)*sky->lut_cell_sz, vox_svx_upp[2]); 487 488 igrid_low[2] = split->index; 489 igrid_upp[2] = split->index + (pos_z > split->height); 490 491 if(igrid_upp[2] >= sky->htcp_desc.spatial_definition[2] 492 && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */ 493 igrid_upp[2] = split->index; 494 } 495 496 FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) { 497 FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) { 498 FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) { 499 double x_h2o; 500 501 /* Does the LUT cell overlap an already handled LES voxel? */ 502 if(ivox[2] == ivox_z_prev) continue; 503 ivox_z_prev = ivox[2]; 504 505 /* Compute the xH2O for the current LES voxel */ 506 x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox); 507 508 /* Update the xH2O range of the SVX voxel */ 509 x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]); 510 x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]); 511 } 512 } 513 } 514 } 515 } 516 517 static void 518 cloud_vox_get_gas 519 (const size_t xyz[3], 520 float vox[], 521 const struct build_tree_context* ctx) 522 { 523 const struct htcp_desc* htcp_desc; 524 struct htgop_layer layer; 525 size_t ilayer; 526 size_t layer_range[2]; 527 double x_h2o_range[2]; 528 double vox_low[3], vox_upp[3]; /* Voxel AABB */ 529 double ka[2], ks[2], kext[2]; 530 ASSERT(xyz && vox && ctx); 531 532 htcp_desc = &ctx->sky->htcp_desc; 533 534 /* Compute the AABB of the SVX voxel */ 535 vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + htcp_desc->lower[0]; 536 vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + htcp_desc->lower[1]; 537 vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + htcp_desc->lower[2]; 538 vox_upp[0] = vox_low[0] + ctx->vxsz[0]; 539 vox_upp[1] = vox_low[1] + ctx->vxsz[1]; 540 vox_upp[2] = vox_low[2] + ctx->vxsz[2]; 541 542 /* Define the xH2O range from the LES data */ 543 if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */ 544 compute_xh2o_range_regular_z(ctx->sky, vox_low, vox_upp, x_h2o_range); 545 } else { /* A SVX voxel can be overlapped by 2 LES voxels */ 546 ASSERT(xyz[2] < darray_split_size_get(&ctx->sky->svx2htcp_z)); 547 compute_xh2o_range_irregular_z(ctx->sky, vox_low, vox_upp, x_h2o_range); 548 } 549 550 /* Define the atmospheric layers overlapped by the SVX voxel */ 551 HTGOP(position_to_layer_id(ctx->sky->htgop, vox_low[2], &layer_range[0])); 552 HTGOP(position_to_layer_id(ctx->sky->htgop, vox_upp[2], &layer_range[1])); 553 554 ka[0] = ks[0] = kext[0] = DBL_MAX; 555 ka[1] = ks[1] = kext[1] =-DBL_MAX; 556 557 /* For each atmospheric layer that overlaps the SVX voxel ... */ 558 if(ctx->sky->spectral_type == HTSKY_SPECTRAL_LW) { 559 FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { 560 struct htgop_layer_lw_spectral_interval band; 561 double k[2]; 562 563 HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer)); 564 565 /* ... retrieve the considered spectral interval */ 566 HTGOP(layer_get_lw_spectral_interval(&layer, ctx->iband, &band)); 567 ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]); 568 ASSERT(ctx->quadrature_range[1] < band.quadrature_length); 569 570 /* ... and compute the radiative properties and upd their bounds */ 571 HTGOP(layer_lw_spectral_interval_quadpoints_get_ka_bounds 572 (&band, ctx->quadrature_range, x_h2o_range, k)); 573 ka[0] = MMIN(ka[0], k[0]); 574 ka[1] = MMAX(ka[1], k[1]); 575 } 576 ks[0] = 0; 577 ks[1] = 0; 578 kext[0] = ka[0]; 579 kext[1] = ka[1]; 580 } else { 581 ASSERT(ctx->sky->spectral_type == HTSKY_SPECTRAL_SW); 582 583 FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { 584 struct htgop_layer_sw_spectral_interval band; 585 double k[2]; 586 587 HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer)); 588 589 /* ... retrieve the considered spectral interval */ 590 HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band)); 591 ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]); 592 ASSERT(ctx->quadrature_range[1] < band.quadrature_length); 593 594 /* ... and compute the radiative properties and upd their bounds */ 595 HTGOP(layer_sw_spectral_interval_quadpoints_get_ka_bounds 596 (&band, ctx->quadrature_range, x_h2o_range, k)); 597 ka[0] = MMIN(ka[0], k[0]); 598 ka[1] = MMAX(ka[1], k[1]); 599 HTGOP(layer_sw_spectral_interval_quadpoints_get_ks_bounds 600 (&band, ctx->quadrature_range, x_h2o_range, k)); 601 ks[0] = MMIN(ks[0], k[0]); 602 ks[1] = MMAX(ks[1], k[1]); 603 HTGOP(layer_sw_spectral_interval_quadpoints_get_kext_bounds 604 (&band, ctx->quadrature_range, x_h2o_range, k)); 605 kext[0] = MMIN(kext[0], k[0]); 606 kext[1] = MMAX(kext[1], k[1]); 607 } 608 } 609 610 /* Ensure that the single precision bounds include their double precision 611 * version. */ 612 if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX); 613 if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX); 614 if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX); 615 if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX); 616 if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX); 617 if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX); 618 619 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka[0]); 620 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka[1]); 621 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks[0]); 622 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks[1]); 623 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext[0]); 624 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext[1]); 625 } 626 627 static void 628 cloud_vox_get 629 (const size_t xyz[3], 630 const uint64_t mcode, 631 void* dst, 632 void* context) 633 { 634 struct build_tree_context* ctx = context; 635 ASSERT(context); 636 (void)mcode; 637 cloud_vox_get_particle(xyz, dst, ctx); 638 cloud_vox_get_gas(xyz, dst, ctx); 639 } 640 641 static void 642 cloud_vox_merge(void* dst, const void* voxels[], const size_t nvoxs, void* context) 643 { 644 ASSERT(dst && voxels && nvoxs); 645 (void)context; 646 vox_merge_component(dst, HTSKY_CPNT_PARTICLES, (const float**)voxels, nvoxs); 647 vox_merge_component(dst, HTSKY_CPNT_GAS, (const float**)voxels, nvoxs); 648 } 649 650 static int 651 cloud_vox_challenge_merge 652 (const struct svx_voxel voxels[], const size_t nvoxs, void* ctx) 653 { 654 ASSERT(voxels); 655 return vox_challenge_merge_component(HTSKY_CPNT_PARTICLES, voxels, nvoxs, ctx) 656 && vox_challenge_merge_component(HTSKY_CPNT_GAS, voxels, nvoxs, ctx); 657 } 658 659 static INLINE void 660 print_progress(struct htsky* sky, const int progress) 661 { 662 ASSERT(sky && progress >= 0 && progress <= 100); 663 log_info(sky, "\033[2K\r"); /* Erase the line */ 664 log_info(sky,"Computing clouds data & building octrees: %3d%%\r" , progress); 665 } 666 667 668 static res_T 669 cloud_build_octrees 670 (struct htsky* sky, 671 const double low[3], /* Lower bound of the clouds */ 672 const double upp[3], /* Upper bound of the clouds */ 673 const unsigned grid_max_definition[3], 674 const double optical_thickness_threshold, 675 struct darray_specdata* specdata) 676 { 677 double vxsz[3] = {0, 0, 0}; 678 size_t nvoxs[3] = {0, 0, 0}; 679 const size_t* raw_def = NULL; 680 int64_t ispecdata; 681 int64_t nspecdata; 682 int32_t progress; 683 ATOMIC nbuilt_octrees = 0; 684 ATOMIC res = RES_OK; 685 ASSERT(sky && specdata); 686 ASSERT(low && upp && grid_max_definition); 687 ASSERT(low[0] < upp[0]); 688 ASSERT(low[1] < upp[1]); 689 ASSERT(low[2] < upp[2]); 690 ASSERT(grid_max_definition[0]); 691 ASSERT(grid_max_definition[1]); 692 ASSERT(grid_max_definition[2]); 693 ASSERT(optical_thickness_threshold >= 0); 694 695 progress = 0; 696 nspecdata = (int64_t)darray_specdata_size_get(specdata); 697 698 /* Define the number of voxels */ 699 raw_def = sky->htcp_desc.spatial_definition; 700 nvoxs[0] = MMIN(raw_def[0], grid_max_definition[0]); 701 nvoxs[1] = MMIN(raw_def[1], grid_max_definition[1]); 702 nvoxs[2] = MMIN(raw_def[2], grid_max_definition[2]); 703 704 /* Setup the build context */ 705 vxsz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0]; 706 vxsz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1]; 707 vxsz[2] = sky->htcp_desc.upper[2] - sky->htcp_desc.lower[2]; 708 vxsz[0] = vxsz[0] / (double)nvoxs[0]; 709 vxsz[1] = vxsz[1] / (double)nvoxs[1]; 710 vxsz[2] = vxsz[2] / (double)nvoxs[2]; 711 712 omp_set_num_threads((int)sky->nthreads); 713 print_progress(sky, 0); 714 #pragma omp parallel for schedule(dynamic, 1/*chunksize*/) 715 for(ispecdata=0; ispecdata < nspecdata; ++ispecdata) { 716 struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL; 717 struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL; 718 const size_t iband = darray_specdata_data_get(specdata)[ispecdata].iband; 719 const size_t iquad = darray_specdata_data_get(specdata)[ispecdata].iquad; 720 const size_t id = iband - sky->bands_range[0]; 721 int32_t pcent; 722 size_t n; 723 res_T res_local = RES_OK; 724 725 if(ATOMIC_GET(&res) != RES_OK) continue; 726 727 /* The current octree could be already setup from a cache and the 728 * invocation of the current function is still valid. Indeed, if some 729 * octrees could not be restore from the cache, one can choose to build 730 * them from scratch rather than rejecting the whole process. */ 731 if(sky->clouds[id][iquad].octree) continue; 732 733 /* Setup the build context */ 734 ctx.sky = sky; 735 ctx.vxsz[0] = vxsz[0]; 736 ctx.vxsz[1] = vxsz[1]; 737 ctx.vxsz[2] = vxsz[2]; 738 ctx.tau_threshold = optical_thickness_threshold; 739 ctx.iband = iband; 740 ctx.quadrature_range[0] = iquad; 741 ctx.quadrature_range[1] = iquad; 742 743 /* Setup the voxel descriptor */ 744 vox_desc.get = cloud_vox_get; 745 vox_desc.merge = cloud_vox_merge; 746 vox_desc.challenge_merge = cloud_vox_challenge_merge; 747 vox_desc.context = &ctx; 748 vox_desc.size = sizeof(float) * NFLOATS_PER_VOXEL; 749 750 /* Create the octree */ 751 res_local = svx_octree_create 752 (sky->svx, low, upp, nvoxs, &vox_desc, &sky->clouds[id][iquad].octree); 753 754 if(res_local != RES_OK) { 755 log_err(sky, 756 "Could not create the octree of the cloud properties for the band %lu.\n", 757 (unsigned long)ctx.iband); 758 ATOMIC_SET(&res, res_local); 759 continue; 760 } 761 762 /* Fetch the octree descriptor for future use */ 763 SVX(tree_get_desc 764 (sky->clouds[id][iquad].octree, &sky->clouds[id][iquad].octree_desc)); 765 766 /* Update the progress message */ 767 n = (size_t)ATOMIC_INCR(&nbuilt_octrees); 768 pcent = (int32_t)(n * 100 / darray_specdata_size_get(specdata)); 769 770 #pragma omp critical 771 if(pcent > progress) { 772 progress = pcent; 773 print_progress(sky, pcent); 774 } 775 } 776 777 if(res != RES_OK) 778 goto error; 779 780 print_progress(sky, 100); 781 log_info(sky, "\n"); /* New line */ 782 783 exit: 784 return (res_T)res; 785 error: 786 goto exit; 787 } 788 789 static res_T 790 cloud_write_octrees(const struct htsky* sky, FILE* stream) 791 { 792 size_t nbands; 793 size_t i; 794 res_T res = RES_OK; 795 ASSERT(sky && sky->clouds); 796 797 nbands = htsky_get_spectral_bands_count(sky); 798 FOR_EACH(i, 0, nbands) { 799 struct htgop_spectral_interval band; 800 size_t iband; 801 size_t iquad; 802 803 iband = sky->bands_range[0] + i; 804 switch(sky->spectral_type) { 805 case HTSKY_SPECTRAL_LW: 806 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 807 break; 808 case HTSKY_SPECTRAL_SW: 809 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 810 break; 811 default: FATAL("Unreachable code.\n"); break; 812 } 813 ASSERT(sky->clouds[i]); 814 815 FOR_EACH(iquad, 0, band.quadrature_length) { 816 ASSERT(sky->clouds[i][iquad].octree); 817 res = svx_tree_write(sky->clouds[i][iquad].octree, stream); 818 if(res != RES_OK) { 819 log_err(sky, "Could not write the octrees of the clouds.\n"); 820 goto error; 821 } 822 } 823 } 824 825 exit: 826 return res; 827 error: 828 goto exit; 829 } 830 831 static res_T 832 cloud_read_octrees(struct htsky* sky, const char* stream_name, FILE* stream) 833 { 834 size_t nbands; 835 size_t i; 836 res_T res = RES_OK; 837 ASSERT(sky && sky->clouds && stream_name); 838 839 log_info(sky, "Loading octrees from `%s'.\n", stream_name); 840 nbands = htsky_get_spectral_bands_count(sky); 841 FOR_EACH(i, 0, nbands) { 842 struct htgop_spectral_interval band; 843 size_t iband; 844 size_t iquad; 845 846 iband = sky->bands_range[0] + i; 847 switch(sky->spectral_type) { 848 case HTSKY_SPECTRAL_LW: 849 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 850 break; 851 case HTSKY_SPECTRAL_SW: 852 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 853 break; 854 default: FATAL("Unreachable code.\n"); break; 855 } 856 ASSERT(sky->clouds[i]); 857 858 FOR_EACH(iquad, 0, band.quadrature_length) { 859 /* The octree was already setup */ 860 if(sky->clouds[i][iquad].octree != NULL) continue; 861 862 res = svx_tree_create_from_stream 863 (sky->svx, stream, &sky->clouds[i][iquad].octree); 864 if(res != RES_OK) { 865 log_err(sky, "Could not read the octree %lu:%lu from `%s' -- %s.\n", 866 (unsigned long)iband, (unsigned long)iquad, stream_name, 867 res_to_cstr(res)); 868 goto error; 869 } 870 871 /* Fetch the octree descriptor for future use */ 872 SVX(tree_get_desc 873 (sky->clouds[i][iquad].octree, &sky->clouds[i][iquad].octree_desc)); 874 } 875 } 876 877 exit: 878 return res; 879 error: 880 goto exit; 881 } 882 883 /******************************************************************************* 884 * Local functions 885 ******************************************************************************/ 886 res_T 887 cloud_setup 888 (struct htsky* sky, 889 const unsigned grid_max_definition[3], 890 const double optical_thickness_threshold, 891 const char* cache_name, 892 const int force_cache_update, 893 FILE* cache) 894 { 895 struct darray_specdata specdata; 896 const size_t* raw_def; 897 double low[3]; 898 double upp[3]; 899 size_t nbands; 900 size_t i; 901 ATOMIC res = RES_OK; 902 ASSERT(sky && sky->bands && optical_thickness_threshold >= 0); 903 ASSERT(grid_max_definition); 904 ASSERT(grid_max_definition[0]); 905 ASSERT(grid_max_definition[1]); 906 ASSERT(grid_max_definition[2]); 907 908 darray_specdata_init(sky->allocator, &specdata); 909 910 res = htcp_get_desc(sky->htcp, &sky->htcp_desc); 911 if(res != RES_OK) { 912 log_err(sky, "Could not retrieve the HTCP descriptor.\n"); 913 goto error; 914 } 915 916 log_info(sky, "Clouds bounding box: {%g, %g, %g} / {%g, %g, %g}.\n", 917 SPLIT3(sky->htcp_desc.lower), SPLIT3(sky->htcp_desc.upper)); 918 919 /* Define the number of voxels */ 920 raw_def = sky->htcp_desc.spatial_definition; 921 922 /* Define the octree AABB excepted for the Z dimension */ 923 low[0] = sky->htcp_desc.lower[0]; 924 low[1] = sky->htcp_desc.lower[1]; 925 low[2] = sky->htcp_desc.lower[2]; 926 upp[0] = low[0] + (double)raw_def[0] * sky->htcp_desc.vxsz_x; 927 upp[1] = low[1] + (double)raw_def[1] * sky->htcp_desc.vxsz_y; 928 929 if(!sky->htcp_desc.irregular_z) { 930 /* Regular voxel size along the Z dimension: compute its upper boundary as 931 * the others dimensions */ 932 upp[2] = low[2] + (double)raw_def[2] * sky->htcp_desc.vxsz_z[0]; 933 } else { /* Irregular voxel size along Z */ 934 res = setup_svx2htcp_z_lut(sky, &upp[2]); 935 if(res != RES_OK) goto error; 936 } 937 938 /* Create as many cloud data structure than considered SW spectral bands */ 939 nbands = htsky_get_spectral_bands_count(sky); 940 sky->clouds = MEM_CALLOC(sky->allocator, nbands, sizeof(*sky->clouds)); 941 if(!sky->clouds) { 942 log_err(sky, 943 "Could not create the list of per SW band cloud data structure.\n"); 944 res = RES_MEM_ERR; 945 goto error; 946 } 947 948 /* Compute how many octrees are going to be built */ 949 FOR_EACH(i, 0, nbands) { 950 struct htgop_spectral_interval band; 951 const size_t iband = i + sky->bands_range[0]; 952 size_t iquad; 953 954 switch(sky->spectral_type) { 955 case HTSKY_SPECTRAL_LW: 956 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 957 break; 958 case HTSKY_SPECTRAL_SW: 959 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 960 break; 961 default: FATAL("Unreachable code.\n"); break; 962 } 963 964 sky->clouds[i] = MEM_CALLOC(sky->allocator, 965 band.quadrature_length, sizeof(*sky->clouds[i])); 966 if(!sky->clouds[i]) { 967 log_err(sky, 968 "Could not create the list of per quadrature point cloud data " 969 "for the band %lu.\n", (unsigned long)iband); 970 res = RES_MEM_ERR; 971 goto error; 972 } 973 974 FOR_EACH(iquad, 0, band.quadrature_length) { 975 struct spectral_data spectral_data; 976 spectral_data.iband = iband; 977 spectral_data.iquad = iquad; 978 res = darray_specdata_push_back(&specdata, &spectral_data); 979 if(res != RES_OK) { 980 log_err(sky, 981 "Could not register the quadrature point %lu of the spectral band " 982 "%lu .\n", (unsigned long)iband, (unsigned long)iquad); 983 goto error; 984 } 985 } 986 } 987 988 if(!cache) { /* No cache => build the octrees from scratch */ 989 res = cloud_build_octrees(sky, low, upp, grid_max_definition, 990 optical_thickness_threshold, &specdata); 991 if(res != RES_OK) goto error; 992 993 } else if(force_cache_update) { 994 /* Build the octrees from scratch */ 995 res = cloud_build_octrees(sky, low, upp, grid_max_definition, 996 optical_thickness_threshold, &specdata); 997 if(res != RES_OK) goto error; 998 999 /* Write the octrees into the cache stream */ 1000 res = cloud_write_octrees(sky, cache); 1001 if(res != RES_OK) goto error; 1002 1003 } else { 1004 const long offset = ftell(cache); 1005 1006 /* Try to load the octrees from the cache */ 1007 res = cloud_read_octrees(sky, cache_name, cache); 1008 if(res != RES_OK) { 1009 const int err = fseek(cache, offset, SEEK_SET); /* Restore the stream */ 1010 if(err) { 1011 log_err(sky, "Could not restore the cache stream.\n"); 1012 res = RES_IO_ERR; 1013 goto error; 1014 } 1015 1016 /* Build the octrees from scratch */ 1017 res = cloud_build_octrees(sky, low, upp, grid_max_definition, 1018 optical_thickness_threshold, &specdata); 1019 if(res != RES_OK) goto error; 1020 1021 /* Write the octrees toward the cache stream */ 1022 res = cloud_write_octrees(sky, cache); 1023 if(res != RES_OK) goto error; 1024 } 1025 } 1026 1027 exit: 1028 darray_specdata_release(&specdata); 1029 return (res_T)res; 1030 error: 1031 cloud_clean(sky); 1032 darray_split_clear(&sky->svx2htcp_z); 1033 goto exit; 1034 } 1035 1036 void 1037 cloud_clean(struct htsky* sky) 1038 { 1039 size_t nbands; 1040 size_t i; 1041 ASSERT(sky); 1042 1043 if(!sky->clouds) return; 1044 1045 nbands = htsky_get_spectral_bands_count(sky); 1046 FOR_EACH(i, 0, nbands) { 1047 struct htgop_spectral_interval band; 1048 size_t iband; 1049 size_t iquad; 1050 1051 iband = sky->bands_range[0] + i; 1052 switch(sky->spectral_type) { 1053 case HTSKY_SPECTRAL_LW: 1054 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 1055 break; 1056 case HTSKY_SPECTRAL_SW: 1057 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 1058 break; 1059 default: FATAL("Unreachable code.\n"); break; 1060 } 1061 1062 if(!sky->clouds[i]) continue; 1063 1064 FOR_EACH(iquad, 0, band.quadrature_length) { 1065 if(sky->clouds[i][iquad].octree) { 1066 SVX(tree_ref_put(sky->clouds[i][iquad].octree)); 1067 sky->clouds[i][iquad].octree = NULL; 1068 } 1069 } 1070 MEM_RM(sky->allocator, sky->clouds[i]); 1071 } 1072 MEM_RM(sky->allocator, sky->clouds); 1073 sky->clouds = NULL; 1074 } 1075