commit 9157d07191280f7915ec2a6947b44d00f864ba77
parent c91f8ebe300d244447dc828766a411c06d0e3002
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 30 Aug 2022 16:18:05 +0200
Kext<min|max> are no longer stored in voxels
They are calculated at runtime from ka<min|max> and ks<min|max> to
reduce voxel memory space.
Diffstat:
4 files changed, 59 insertions(+), 86 deletions(-)
diff --git a/src/rnatm.c b/src/rnatm.c
@@ -259,24 +259,39 @@ rnatm_get_k_svx_voxel
const enum rnatm_svx_op op)
{
const float* vx;
- float k_min;
- float k_max;
ASSERT(atm && voxel && radcoef);
ASSERT((unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
ASSERT((unsigned)op < RNATM_SVX_OPS_COUNT__);
(void)atm;
vx = voxel->data;
- k_min = vx[voxel_idata(radcoef, RNATM_SVX_OP_MIN)];
- k_max = vx[voxel_idata(radcoef, RNATM_SVX_OP_MAX)];
- /* Return a zero radiative coefficient for empty voxels */
- if(k_min > k_max) return 0;
+ /* Absorption/Scattering coefficient */
+ if(radcoef != RNATM_RADCOEF_Kext) {
+ const float k_min = vx[voxel_idata(radcoef, RNATM_SVX_OP_MIN)];
+ const float k_max = vx[voxel_idata(radcoef, RNATM_SVX_OP_MAX)];
+ if(k_min > k_max) return 0; /* Empty voxel => null radiative coefficient */
+ switch(op) {
+ case RNATM_SVX_OP_MIN: return k_min;
+ case RNATM_SVX_OP_MAX: return k_max;
+ default: FATAL("Unreachable code\n"); break;
+ }
- switch(op) {
- case RNATM_SVX_OP_MIN: return k_min;
- case RNATM_SVX_OP_MAX: return k_max;
- default: FATAL("Unreachable code\n"); break;
+ /* Extinction coefficient */
+ } else {
+ const float ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
+ const float ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
+ const float ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
+ const float ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
+
+ /* Empty voxel => null radiative coefficient */
+ if(ka_min > ka_max) { ASSERT(ks_min > ks_max); return 0; }
+
+ switch(op) {
+ case RNATM_SVX_OP_MIN: return ka_min + ks_min;
+ case RNATM_SVX_OP_MAX: return ka_max + ks_max;
+ default: FATAL("Unreachable code\n"); break;
+ }
}
}
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -70,7 +70,6 @@ static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
struct radcoefs {
float ka_min, ka_max;
float ks_min, ks_max;
- float kext_min, kext_max;
};
#define RADCOEFS_NULL__ { \
FLT_MAX, -FLT_MAX, \
@@ -266,20 +265,20 @@ find_band_range
nbands = sck_get_bands_count(atm->gas.ck);
- /* Find the lower bound */
+ /* Find the lower bound (inclusive) */
FOR_EACH(ilow, 0, nbands) {
SCK(get_band(atm->gas.ck, ilow, &band_low));
if(band_low.upper >= range[0]) break;
}
- /* Find the upper bound */
+ /* Find the upper bound (inclusive) */
FOR_EACH(iupp, ilow, nbands) {
SCK(get_band(atm->gas.ck, iupp, &band_upp));
- if(band_upp.lower > range[1]) break;
+ if(band_upp.upper > range[1]) break;
}
bands[0] = ilow;
- bands[1] = iupp - 1; /* Make the boundary inclusive */
+ bands[1] = iupp; /* Make the boundary inclusive */
if(bands[0] > bands[1]) {
log_err(atm,
@@ -462,7 +461,6 @@ setup_tetra_radcoefs_gas
{
float ka[4];
float ks[4];
- float kext[4];
struct sck_band band;
struct sck_quad_pt quad_pt;
ASSERT(atm && tetra && radcoefs);
@@ -485,14 +483,6 @@ setup_tetra_radcoefs_gas
ka[3] = quad_pt.ka_list[tetra->indices[3]];
radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
-
- /* Compute the extinction coefficient range of the tetrahedron */
- kext[0] = ka[0] + ks[0];
- kext[1] = ka[1] + ks[1];
- kext[2] = ka[2] + ks[2];
- kext[3] = ka[3] + ks[3];
- radcoefs->kext_min = MMIN(MMIN(kext[0], kext[1]), MMIN(kext[2], kext[3]));
- radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
}
static void
@@ -507,7 +497,7 @@ setup_tetra_radcoefs_aerosol
struct sars_band ars_band;
double gas_spectral_range[2]; /* In nm */
size_t ars_ibands[2];
- float ka[4], ks[4], kext[4];
+ float ka[4], ks[4];
ASSERT(atm && aerosol && tetra && radcoefs);
/* Look for the aerosol bands covered by the gas band */
@@ -522,8 +512,6 @@ setup_tetra_radcoefs_aerosol
radcoefs->ks_max =-FLT_MAX;
radcoefs->ka_min = FLT_MAX;
radcoefs->ka_max =-FLT_MAX;
- radcoefs->kext_min = FLT_MAX;
- radcoefs->kext_max =-FLT_MAX;
return;
}
@@ -551,29 +539,22 @@ setup_tetra_radcoefs_aerosol
radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
- /* Compute the extinction coefficient range of the tetrahedron */
- kext[0] = ka[0] + ks[0];
- kext[1] = ka[1] + ks[1];
- kext[2] = ka[2] + ks[2];
- kext[3] = ka[3] + ks[3];
- radcoefs->kext_min = MMIN(MMIN(kext[0], kext[1]), MMIN(kext[2], kext[3]));
- radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
-
/* The gas band overlaid N aerosol bands (N >= 1) */
} else {
float tau_ka[4] = {0, 0, 0, 0};
float tau_ks[4] = {0, 0, 0, 0};
+ double lambda_min;
+ double lambda_max;
float rcp_gas_band_len;
size_t iars_band;
FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) {
- double lambda_min;
- double lambda_max;
- float lambda_len;
+ float lambda_len;
SARS(get_band(aerosol->sars, iars_band, &ars_band));
lambda_min = MMAX(gas_band.lower, ars_band.lower);
- lambda_max = MMIN(gas_band.upper, ars_band.upper);
+ lambda_max = MMIN(gas_band.upper, ars_band.upper); /* exclusive */
+ lambda_max = nextafter(lambda_max, 0); /* inclusive */
lambda_len = (float)(lambda_max - lambda_min);
ASSERT(lambda_len > 0);
@@ -584,24 +565,22 @@ setup_tetra_radcoefs_aerosol
tau_ka[0] += ars_band.ka_list[tetra->indices[0]] * lambda_len;
tau_ka[1] += ars_band.ka_list[tetra->indices[1]] * lambda_len;
- tau_ka[2] += ars_band.ka_list[tetra->indices[2]] * lambda_len;
- tau_ka[3] += ars_band.ka_list[tetra->indices[3]] * lambda_len;
+ tau_ka[2] += ars_band.ka_list[tetra->indices[2]] * lambda_len;
+ tau_ka[3] += ars_band.ka_list[tetra->indices[3]] * lambda_len;
}
/* Compute the radiative coefficients of the tetrahedron */
- ASSERT(gas_band.upper > gas_band.lower);
- rcp_gas_band_len = 1.f / (float)(gas_band.upper - gas_band.lower);
+ lambda_min = gas_band.lower;
+ lambda_max = nextafter(gas_band.upper, 0); /* inclusive */
+ rcp_gas_band_len = 1.f/(float)(lambda_max - lambda_min);
f4_mulf(ks, tau_ks, rcp_gas_band_len);
f4_mulf(ka, tau_ka, rcp_gas_band_len);
- f4_add(kext, ka, ks);
/* Compute the radiative coefficients range of the tetrahedron */
radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
- radcoefs->kext_min = MMIN(MMIN(kext[0], kext[1]), MMIN(kext[2], kext[3]));
- radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
}
}
@@ -635,7 +614,6 @@ update_voxel
{
float vx_ka_min, vx_ka_max;
float vx_ks_min, vx_ks_max;
- float vx_kext_min, vx_kext_max;
float* vx = NULL;
ASSERT(atm && radcoefs && part);
@@ -648,20 +626,14 @@ update_voxel
vx_ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
vx_ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
vx_ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
- vx_kext_min = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)];
- vx_kext_max = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)];
vx_ka_min = MMIN(vx_ka_min, radcoefs->ka_min);
vx_ka_max = MMAX(vx_ka_max, radcoefs->ka_max);
vx_ks_min = MMIN(vx_ks_min, radcoefs->ks_min);
vx_ks_max = MMAX(vx_ks_max, radcoefs->ks_max);
- vx_kext_min = MMIN(vx_kext_min, radcoefs->kext_min);
- vx_kext_max = MMAX(vx_kext_max, radcoefs->kext_max);
vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)] = vx_ka_min;
vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)] = vx_ka_max;
vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)] = vx_ks_min;
vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] = vx_ks_max;
- vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)] = vx_kext_min;
- vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)] = vx_kext_max;
}
static res_T
@@ -1122,7 +1094,6 @@ vx_merge(void* dst, const void* vxs[], const size_t nvxs, void* ctx)
{
float ka_min = FLT_MAX, ka_max = -FLT_MAX;
float ks_min = FLT_MAX, ks_max = -FLT_MAX;
- float kext_min = FLT_MAX, kext_max = -FLT_MAX;
float* merged_vx = dst;
size_t ivx;
ASSERT(dst && vxs && nvxs);
@@ -1134,15 +1105,12 @@ vx_merge(void* dst, const void* vxs[], const size_t nvxs, void* ctx)
ka_max = MMAX(ka_max, vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)]);
ks_min = MMIN(ks_min, vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)]);
ks_max = MMAX(ks_max, vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)]);
- kext_min = MMIN(kext_min, vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)]);
- kext_max = MMAX(kext_max, vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)]);
}
+
merged_vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)] = ka_min;
merged_vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)] = ka_max;
merged_vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)] = ks_min;
merged_vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] = ks_max;
- merged_vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)] = kext_min;
- merged_vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)] = kext_max;
}
static int
@@ -1163,8 +1131,13 @@ vx_challenge_merge
/* Compute the range of the extinction coefficients of the submitted voxels */
FOR_EACH(ivx, 0, nvxs) {
const float* vx = vxs[ivx].data;
- kext_min = MMIN(kext_min, vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)]);
- kext_max = MMAX(kext_max, vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)]);
+ const float ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
+ const float ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
+ const float ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
+ const float ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
+ if(ka_min > ka_max) { ASSERT(ks_min > ks_max); continue; } /* Empty voxel */
+ kext_min = MMIN(kext_min, ka_min + ks_min);
+ kext_max = MMAX(kext_max, ka_max + ks_max);
}
/* Always merge empty voxels */
diff --git a/src/rnatm_voxel.h b/src/rnatm_voxel.h
@@ -27,22 +27,16 @@
* Memory layout of a voxel
* ------------------------
*
- * The data of a voxel are stored in a list of N single-precision floating-point
- * numbers, with N defined as below:
- *
- * N = RNATM_RADCOEFS_COUNT__ * RNATM_SVX_OPS_COUNT__
- *
- * RNATM_RADCOEFS_COUNT__ is the number of radiative coefficients to be stored
- * per voxel and RNATM_SVX_OPS_COUNT__ the number of operations for these
- * coefficients. For a given voxel, the data corresponding to the operation 'O'
- * on the coefficient 'C' are stored at the index 'id' between [0, N-1] and
- * calculated as follows:
+ * The data of a voxel are stored in a list of 4 single-precision floating-point
+ * numbers ka_min, ka_max, ks_min, ks_max. For a given voxel, the data
+ * corresponding to the operation 'O' on the coefficient 'C' are stored at the
+ * index 'id' between [0, N-1] and calculated as follows:
*
* id = C * RNATM_SVX_OPS_COUNT__ + O
*/
/* Total number of floating-point numbers per voxel */
-#define NFLOATS_PER_VOXEL (RNATM_RADCOEFS_COUNT__ * RNATM_SVX_OPS_COUNT__)
+#define NFLOATS_PER_VOXEL 4
/* Calculate the data index of voxel */
static FINLINE size_t
@@ -50,7 +44,7 @@ voxel_idata
(const enum rnatm_radcoef radcoef,
const enum rnatm_svx_op op)
{
- ASSERT((unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
+ ASSERT(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Ks);
ASSERT((unsigned)op < RNATM_SVX_OPS_COUNT__);
return radcoef*RNATM_SVX_OPS_COUNT__ + op;
}
@@ -58,20 +52,10 @@ voxel_idata
static FINLINE void
voxel_clear(float* voxel)
{
- int radcoef, op;
-
- /* Make degenerated the radcoef range */
- FOR_EACH(radcoef, 0, RNATM_RADCOEFS_COUNT__) {
- FOR_EACH(op, 0, RNATM_SVX_OPS_COUNT__) {
- float val = 0;
- switch(op) {
- case RNATM_SVX_OP_MIN: val = FLT_MAX; break;
- case RNATM_SVX_OP_MAX: val =-FLT_MAX; break;
- default: FATAL("Unreachable code\n"); break;
- }
- voxel[voxel_idata(radcoef, op)] = val;
- }
- }
+ voxel[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)] = FLT_MAX;
+ voxel[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)] =-FLT_MAX;
+ voxel[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)] = FLT_MAX;
+ voxel[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] =-FLT_MAX;
}
#endif /* RNATM_VOXEL_H */
diff --git a/src/rnatm_write_vtk_octree.c b/src/rnatm_write_vtk_octree.c
@@ -139,6 +139,7 @@ register_leaf
/* Register leaf data */
kext_min = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN);
kext_max = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX);
+ ASSERT(kext_min <= kext_max);
CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
}
