commit 29a0cf998c691e1221fc13db205e211f06e91a0a
parent c5300900828bcb39f5e83bfbc5609c0421a726a5
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 8 Jun 2020 18:18:24 +0200
Fix data access during tree construction
Some cloud data were not into taken account during the construction of
the hierarchical data structure when they were irregularly structured
along the Z axis. This might lead to k_min/k_max that did not really
encompass the underlying k coefficients.
Diffstat:
2 files changed, 54 insertions(+), 52 deletions(-)
diff --git a/src/htsky.c b/src/htsky.c
@@ -825,7 +825,7 @@ htsky_fetch_raw_property
} else {
/* Irregular voxel size along the Z dimension. Compute the index of the Z
* position in the svx2htcp_z Look Up Table and use the LUT to define the
- * voxel index into the HTCP descripptor */
+ * voxel index into the HTCP descriptor */
const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
const size_t ilut = (size_t)
((pos_cs[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
diff --git a/src/htsky_cloud.c b/src/htsky_cloud.c
@@ -203,7 +203,7 @@ compute_k_bounds_irregular_z
double kext_bounds[2])
{
size_t ivox[3];
- size_t igrid_low[2], igrid_upp[2];
+ size_t igrid_low[3], igrid_upp[3];
size_t i;
double ka, ks, kext;
double low[2], upp[2];
@@ -257,40 +257,41 @@ compute_k_bounds_irregular_z
const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z)+ilut;
ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z));
- ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
- if(ivox[2] >= sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
- ivox[2] = split->index;
- }
+ /* Compute the upper bound of the current LUT cell clamped to the voxel
+ * upper bound. */
+ pos_z = MMIN((double)(ilut+1)*sky->lut_cell_sz, vox_svx_upp[2]);
- /* Compute the upper bound of the *next* LUT cell clamped to the voxel
- * upper bound. Note that the upper bound of the current LUT cell is
- * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
- * upper bound of the next cell is thus the lower bound of the cell
- * following the next cell, i.e. (ilut+2)*lut_cell_sz */
- pos_z = MMIN((double)(ilut+2)*sky->lut_cell_sz, vox_svx_upp[2]);
+ igrid_low[2] = split->index;
+ igrid_upp[2] = split->index + (pos_z > split->height);
- /* Does the LUT cell overlap an already handled LES voxel? */
- if(ivox[2] == ivox_z_prev) continue;
- ivox_z_prev = ivox[2];
+ if(igrid_upp[2] >= sky->htcp_desc.spatial_definition[2]
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ igrid_upp[2] = split->index;
+ }
FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
-
- /* Compute the radiative properties */
- rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
- rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ka = Ca_avg * rho_da * rct;
- ks = Cs_avg * rho_da * rct;
- kext = ka + ks;
-
- /* Update the boundaries of the radiative properties */
- ka_bounds[0] = MMIN(ka_bounds[0], ka);
- ka_bounds[1] = MMAX(ka_bounds[1], ka);
- ks_bounds[0] = MMIN(ks_bounds[0], ks);
- ks_bounds[1] = MMAX(ks_bounds[1], ks);
- kext_bounds[0] = MMIN(kext_bounds[0], kext);
- kext_bounds[1] = MMAX(kext_bounds[1], kext);
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+
+ /* Does the LUT cell overlap an already handled LES voxel? */
+ if(ivox[2] == ivox_z_prev) continue;
+ ivox_z_prev = ivox[2];
+
+ /* Compute the radiative properties */
+ rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
+ rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
+ kext = ka + ks;
+
+ /* Update the boundaries of the radiative properties */
+ ka_bounds[0] = MMIN(ka_bounds[0], ka);
+ ka_bounds[1] = MMAX(ka_bounds[1], ka);
+ ks_bounds[0] = MMIN(ks_bounds[0], ks);
+ ks_bounds[1] = MMAX(ks_bounds[1], ks);
+ kext_bounds[0] = MMIN(kext_bounds[0], kext);
+ kext_bounds[1] = MMAX(kext_bounds[1], kext);
+ }
}
}
}
@@ -408,7 +409,7 @@ compute_xh2o_range_irregular_z
double x_h2o_range[2])
{
size_t ivox[3];
- size_t igrid_low[2], igrid_upp[2];
+ size_t igrid_low[2], igrid_upp[3];
size_t ilut_low, ilut_upp;
size_t ilut;
size_t ivox_z_prev;
@@ -450,33 +451,34 @@ compute_xh2o_range_irregular_z
const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z) + ilut;
ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z));
- ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
- if(ivox[2] >= sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
- ivox[2] = split->index;
- }
+ /* Compute the upper bound of the current LUT cell clamped to the voxel
+ * upper bound.*/
+ pos_z = MMIN((double)(ilut+1)*sky->lut_cell_sz, vox_svx_upp[2]);
- /* Compute the upper bound of the *next* LUT cell clamped to the voxel
- * upper bound. Note that the upper bound of the current LUT cell is
- * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
- * upper bound of the next cell is thus the lower bound of the cell
- * following the next cell, i.e. (ilut+2)*lut_cell_sz */
- pos_z = MMIN((double)(ilut+2)*sky->lut_cell_sz, vox_svx_upp[2]);
+ igrid_low[2] = split->index;
+ igrid_upp[2] = split->index + (pos_z > split->height);
- /* Does the LUT voxel overlap an already handled LES voxel? */
- if(ivox[2] == ivox_z_prev) continue;
- ivox_z_prev = ivox[2];
+ if(igrid_upp[2] >= sky->htcp_desc.spatial_definition[2]
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ igrid_upp[2] = split->index;
+ }
FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
- double x_h2o;
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+ double x_h2o;
- /* Compute the xH2O for the current LES voxel */
- x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+ /* Does the LUT cell overlap an already handled LES voxel? */
+ if(ivox[2] == ivox_z_prev) continue;
+ ivox_z_prev = ivox[2];
- /* Update the xH2O range of the SVX voxel */
- x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
- x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ /* Compute the xH2O for the current LES voxel */
+ x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+
+ /* Update the xH2O range of the SVX voxel */
+ x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
+ x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ }
}
}
}
