commit 190039c9b64654e19a148f2ca4c964f32f02e9b0
parent cd26e116bc62d133a77007de6345fdba519bcc86
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 7 Nov 2022 15:29:38 +0100
Fix how atmospheric components are voxelized
Until now, a voxel stored the maximum and minimum of the radiative
coefficients of atmospheric components while it had to store the maximum
and minimum of their sum. In this commit, we use a dummy in which
aerosols are voxelized. This partition is then accumulated in the
current partition in which the gas is directly voxelized.
Diffstat:
4 files changed, 163 insertions(+), 64 deletions(-)
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -81,10 +81,16 @@ struct radcoefs {
#define DARRAY_DATA struct radcoefs*
#include <rsys/dynamic_array.h>
+/* Generate the dynmaic array of partitions */
+#define DARRAY_NAME partition
+#define DARRAY_DATA struct partition*
+#include <rsys/dynamic_array.h>
+
struct voxelize_batch_args {
/* Working data structure */
struct darray_size_t_list* per_thread_tetra_list;
struct darray_radcoef_list* per_thread_radcoef_list;
+ struct darray_partition* per_thread_dummy_partition;
struct pool* pool;
size_t part_def; /* Partition definition */
@@ -105,6 +111,7 @@ struct voxelize_batch_args {
#define VOXELIZE_BATCH_ARGS_NULL__ { \
NULL, /* Per thread tetra list */ \
NULL, /* Per thread radiative coefs */ \
+ NULL, /* Per thread dummy partition */ \
NULL, /* Partition pool */ \
\
0, /* Partition definition */ \
@@ -130,6 +137,7 @@ struct voxelize_args {
struct radcoefs* per_item_radcoefs;
struct partition* part;
+ struct partition* part_dummy;
size_t part_def; /* Partition definition */
/* AABB of the partition */
@@ -147,6 +155,7 @@ struct voxelize_args {
NULL, /* Radiative coefs */ \
\
NULL, /* Partition */ \
+ NULL, /* Dummy partition */ \
0, /* Partition definition */ \
\
/* AABB of the partition */ \
@@ -470,6 +479,7 @@ static res_T
voxelize_tetra
(struct rnatm* atm,
const struct voxelize_args* args,
+ struct partition* part, /* Partition to fill */
const size_t cpnt)
{
const struct suvm_volume* volume = NULL;
@@ -558,7 +568,7 @@ voxelize_tetra
/* Update the intersected voxel */
FOR_EACH(iitem, 0, args->batch_size) {
struct radcoefs* radcoefs = args->per_item_radcoefs + iitem;
- update_voxel(atm, radcoefs, args->part, mcode[0], iitem);
+ update_voxel(atm, radcoefs, part, mcode[0], iitem);
}
}
}
@@ -569,65 +579,65 @@ voxelize_tetra
}
static res_T
-voxelize_partition_component
- (struct rnatm* atm,
- const struct voxelize_args* args,
- const size_t cpnt)
-{
- struct suvm_volume* volume = NULL;
- res_T res = RES_OK;
- ASSERT(atm && args);
- ASSERT(cpnt == RNATM_GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
-
- /* Get the volumetric mesh of the component */
- if(cpnt == RNATM_GAS) {
- volume = atm->gas.volume;
- } else {
- volume = darray_aerosol_cdata_get(&atm->aerosols)[cpnt].volume;
- }
-
- /* Find the list of gas tetrahedra that overlap the partition */
- darray_size_t_clear(args->tetra_ids);
- res = suvm_volume_intersect_aabb
- (volume, args->part_low, args->part_upp, register_tetra, args->tetra_ids);
- if(res != RES_OK) goto error;
-
- if(darray_size_t_size_get(args->tetra_ids) != 0) {
- res = voxelize_tetra(atm, args, cpnt);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
voxelize_partition(struct rnatm* atm, const struct voxelize_args* args)
{
- size_t iaerosol = 0;
+ STATIC_ASSERT((RNATM_GAS+1) == 0, Unexpected_constant_value);
+
size_t naerosols = 0;
+ size_t cpnt = 0;
int part_is_empty = 1;
res_T res = RES_OK;
ASSERT(atm && check_voxelize_args(args) == RES_OK);
partition_clear_voxels(args->part);
- /* Voxelize the gas */
- res = voxelize_partition_component(atm, args, RNATM_GAS);
- if(res != RES_OK) goto error;
- part_is_empty = part_is_empty && !darray_size_t_size_get(args->tetra_ids);
-
- /* Voxelize the aerosols */
+ /* Voxelize atmospheric components */
naerosols = darray_aerosol_size_get(&atm->aerosols);
- FOR_EACH(iaerosol, 0, naerosols) {
- res = voxelize_partition_component(atm, args, iaerosol);
+
+ /* CAUTION: in the following loop, it is assumed that the first component is
+ * necessarily the gas which is therefore voxelized directly into the
+ * partition. Aerosols are voxelized in a dummy partition before their
+ * accumulation in the partition */
+ cpnt = RNATM_GAS;
+
+ do {
+ struct suvm_volume* volume = NULL;
+
+ /* Get the volumetric mesh of the component */
+ if(cpnt == RNATM_GAS) {
+ volume = atm->gas.volume;
+ } else {
+ volume = darray_aerosol_cdata_get(&atm->aerosols)[cpnt].volume;
+ }
+
+ /* Find the list of tetrahedra that overlap the partition */
+ darray_size_t_clear(args->tetra_ids);
+ res = suvm_volume_intersect_aabb
+ (volume, args->part_low, args->part_upp, register_tetra, args->tetra_ids);
if(res != RES_OK) goto error;
- part_is_empty = part_is_empty && !darray_size_t_size_get(args->tetra_ids);
- }
- /* The partition is not overlaped by any tetrahedron */
+ /* The partition is not covered by any tetrahedron of the component */
+ if(darray_size_t_size_get(args->tetra_ids) == 0) continue;
+
+ /* The partition is covered by at least one tetrahedron and is therefore no
+ * longer empty */
+ part_is_empty = 0;
+
+ /* Voxelize the gas directly into the partition */
+ if(cpnt == RNATM_GAS) {
+ res = voxelize_tetra(atm, args, args->part, cpnt);
+ if(res != RES_OK) goto error;
+
+ /* Voxelize each aerosol into the dummy partition before its accumulation */
+ } else {
+ partition_clear_voxels(args->part_dummy);
+ res = voxelize_tetra(atm, args, args->part_dummy, cpnt);
+ if(res != RES_OK) goto error;
+ partition_accum(args->part, args->part_dummy);
+ }
+ } while(++cpnt < naerosols);
+
+ /* The partition is not covered by any tetrahedron */
if(part_is_empty) partition_empty(args->part);
exit:
@@ -656,6 +666,7 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
#pragma omp parallel for schedule(static, 1/*chunk size*/)
for(i = 0; i < (int64_t)args->nparts_adjusted; ++i) {
struct voxelize_args voxelize_args = VOXELIZE_ARGS_NULL;
+ struct partition* part_dummy = NULL;
struct darray_size_t* tetra_ids = NULL;
struct partition* part = NULL;
struct radcoefs* per_item_radcoefs = NULL;
@@ -671,7 +682,10 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
/* Recover the batch partitions */
part = pool_next_partition(args->pool);
- if(!part) { ATOMIC_SET(&res, RES_UNKNOWN_ERR); continue; };
+ if(!part) {
+ ATOMIC_SET(&res, RES_UNKNOWN_ERR);
+ continue;
+ };
/* Is the current partition out of bounds of the atmosphere grid */
morton_xyz_decode_u21((uint64_t)partition_get_id(part), part_ids);
@@ -702,6 +716,11 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
per_item_radcoefs = darray_radcoef_list_data_get
(args->per_thread_radcoef_list)[ithread];
+ /* Obtain the dummy partition needed to temporarily store aerosol
+ * voxelization prior to accumulation in the current partition */
+ part_dummy = darray_partition_data_get
+ (args->per_thread_dummy_partition)[ithread];
+
/* Voxelizes the partition and once done, commits */
voxelize_args.tetra_ids = tetra_ids;
voxelize_args.part_def = args->part_def;
@@ -709,6 +728,7 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
d3_set(voxelize_args.part_upp, part_upp);
d3_set(voxelize_args.vxsz, args->vxsz);
voxelize_args.part = part;
+ voxelize_args.part_dummy = part_dummy;
voxelize_args.per_item_radcoefs = per_item_radcoefs;
voxelize_args.accel_structs = args->accel_structs;
voxelize_args.batch_size = args->batch_size;
@@ -753,6 +773,7 @@ voxelize_atmosphere
/* Working data structures */
struct darray_size_t_list per_thread_tetra_list;
struct darray_radcoef_list per_thread_radcoef_list;
+ struct darray_partition per_thread_dummy_partition;
/* Voxelization parameters */
double atm_low[3];
@@ -773,12 +794,15 @@ voxelize_atmosphere
darray_size_t_list_init(atm->allocator, &per_thread_tetra_list);
darray_radcoef_list_init(atm->allocator, &per_thread_radcoef_list);
+ darray_partition_init(atm->allocator, &per_thread_dummy_partition);
/* Allocate the per thread lists */
res = darray_size_t_list_resize(&per_thread_tetra_list, atm->nthreads);
if(res != RES_OK) goto error;
res = darray_radcoef_list_resize(&per_thread_radcoef_list, atm->nthreads);
if(res != RES_OK) goto error;
+ res = darray_partition_resize(&per_thread_dummy_partition, atm->nthreads);
+ if(res != RES_OK) goto error;
/* Setup the per thread and per item working data structures */
batch_size = pool_get_voxel_width(pool);
@@ -789,6 +813,12 @@ voxelize_atmosphere
darray_radcoef_list_data_get(&per_thread_radcoef_list)[i] = per_item_k;
}
+ /* Setup the per thread dummy partition */
+ FOR_EACH(i, 0, atm->nthreads) {
+ struct partition* part = pool_dummy_partition(pool);
+ darray_partition_data_get(&per_thread_dummy_partition)[i] = part;
+ }
+
/* Recover the AABB atmosphere and compute the size of a voxel */
SUVM(volume_get_aabb(atm->gas.volume, atm_low, atm_upp));
vxsz[0] = (atm_upp[0] - atm_low[0]) / (double)atm->grid_definition[0];
@@ -827,6 +857,7 @@ voxelize_atmosphere
/* Setup regular batch arguments */
batch_args.per_thread_tetra_list = &per_thread_tetra_list;
batch_args.per_thread_radcoef_list = &per_thread_radcoef_list;
+ batch_args.per_thread_dummy_partition = &per_thread_dummy_partition;
batch_args.pool = pool;
batch_args.part_def = part_def;
batch_args.nparts[0] = nparts[0];
@@ -874,11 +905,18 @@ voxelize_atmosphere
exit:
FOR_EACH(i, 0, darray_radcoef_list_size_get(&per_thread_radcoef_list)) {
- MEM_RM(atm->allocator,
- darray_radcoef_list_data_get(&per_thread_radcoef_list)[i]);
+ struct radcoefs* per_item_k = NULL;
+ per_item_k = darray_radcoef_list_data_get(&per_thread_radcoef_list)[i];
+ MEM_RM(atm->allocator,per_item_k);
+ }
+ FOR_EACH(i, 0, darray_partition_size_get(&per_thread_dummy_partition)) {
+ struct partition* part = NULL;
+ part = darray_partition_data_get(&per_thread_dummy_partition)[i];
+ partition_free(part);
}
darray_size_t_list_release(&per_thread_tetra_list);
darray_radcoef_list_release(&per_thread_radcoef_list);
+ darray_partition_release(&per_thread_dummy_partition);
return res;
error:
goto exit;
@@ -1236,12 +1274,21 @@ create_pool(struct rnatm* atm, struct pool** out_pool, const size_t nbands)
* octrees from the partitions and the threads that fill these partitions */
const size_t NPARTITIONS_PER_THREAD = 64;
+ /* Number of dummy partitions per thread. These partitions are used to
+ * temporarily store the voxelization of aerosols that are then accumulated
+ * in the voxelized partition. Actually, one such partition per voxelization
+ * thread is required */
+ size_t NPARTITIONS_PER_THREAD_DUMMY = 1;
+
/* Empirical constant that defines the total number of partitions to
- * pre-allocate per thread. This constant is necessarily greater than or equal
- * to NPARTITIONS_PER_THREAD, the difference representing the number of
- * completely empty partitions. Such partitions help reduce thread contention
- * with a sparse grid without significantly increasing memory usage */
- const size_t OVERALL_NPARTITIONS_PER_THREAD = NPARTITIONS_PER_THREAD * 64;
+ * pre-allocate per thread. This constant is necessarily greater than or
+ * equal to (NPARTITIONS_PER_THREAD + NPARTITIONS_PER_THREAD_DUMMY), the
+ * difference representing the number of completely empty partitions. Such
+ * partitions help reduce thread contention with a sparse grid without
+ * significantly increasing memory usage */
+ const size_t OVERALL_NPARTITIONS_PER_THREAD =
+ NPARTITIONS_PER_THREAD * 64
+ + NPARTITIONS_PER_THREAD_DUMMY;
/* Number of items stored per voxel data. A width greater than 1 makes it
* possible to store by partition the radiative coefficients of several
diff --git a/src/rnatm_voxel.h b/src/rnatm_voxel.h
@@ -23,8 +23,8 @@
#include "rnatm.h"
-/*
- * Memory layout of a voxel
+/*
+ * Memory layout of a voxel
* ------------------------
*
* The data of a voxel are stored in a list of 4 single-precision floating-point
@@ -58,4 +58,16 @@ voxel_clear(float* voxel)
voxel[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)] =-FLT_MAX;
}
+static FINLINE void
+voxel_accum(float* dst, const float* src)
+{
+ #define ACCUM(K, Op) \
+ dst[voxel_idata((K), (Op))] += src[voxel_idata((K), (Op))]
+ ACCUM(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN);
+ ACCUM(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX);
+ ACCUM(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN);
+ ACCUM(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX);
+ #undef ACCUM
+}
+
#endif /* RNATM_VOXEL_H */
diff --git a/src/rnatm_voxel_partition.c b/src/rnatm_voxel_partition.c
@@ -374,6 +374,24 @@ partition_clear_voxels(struct partition* partition)
}
}
+void
+partition_accum(struct partition* dst, struct partition* src)
+{
+ size_t ivoxel = 0;
+ ASSERT(dst->pool->partition_nvoxels == src->pool->partition_nvoxels);
+ ASSERT(dst->pool->voxel_width == src->pool->voxel_width);
+ ASSERT(src->tile && dst->tile); /* Partition cannot be empty */
+
+ FOR_EACH(ivoxel, 0, src->pool->partition_nvoxels) {
+ size_t iitem = 0;
+ FOR_EACH(iitem, 0, src->pool->voxel_width) {
+ const float* src_voxel = partition_cget_voxel(src, ivoxel, iitem);
+ float* dst_voxel = partition_get_voxel(dst, ivoxel, iitem);
+ voxel_accum(dst_voxel, src_voxel);
+ }
+ }
+}
+
/*******************************************************************************
* Pool of partitions
******************************************************************************/
@@ -488,6 +506,23 @@ struct partition*
pool_next_partition(struct pool* pool)
{
struct partition* partition;
+ ASSERT(pool);
+
+ partition = pool_dummy_partition(pool);
+ if(!partition) return NULL;
+
+ mutex_lock(pool->mutex);
+ partition->id = pool->next_part_id;
+ pool->next_part_id += 1;
+ mutex_unlock(pool->mutex);
+
+ return partition;
+}
+
+struct partition*
+pool_dummy_partition(struct pool* pool)
+{
+ struct partition* partition;
res_T res = RES_OK;
ASSERT(pool);
@@ -497,12 +532,6 @@ pool_next_partition(struct pool* pool)
partition_free(partition);
return NULL;
}
-
- mutex_lock(pool->mutex);
- partition->id = pool->next_part_id;
- pool->next_part_id += 1;
- mutex_unlock(pool->mutex);
-
return partition;
}
diff --git a/src/rnatm_voxel_partition.h b/src/rnatm_voxel_partition.h
@@ -90,6 +90,11 @@ extern LOCAL_SYM void
partition_clear_voxels
(struct partition* partition);
+extern LOCAL_SYM void
+partition_accum
+ (struct partition* dst,
+ struct partition* src);
+
/******************************************************************************
* Partition pool, i.e. collection of partitions of voxels that are accessible
* concurrently by several threads
@@ -123,6 +128,12 @@ extern LOCAL_SYM struct partition*
pool_next_partition
(struct pool* pool);
+/* Same as pool_next_partition but the returned partition has no ID (it
+ * cannot be commited). Such a partition is used as a temporary variable */
+extern LOCAL_SYM struct partition*
+pool_dummy_partition
+ (struct pool* pool);
+
/* Returns the partition with the identifier 'ipartition'. Waits for the
* partition to be available. Returns NULL if an error occurs */
extern LOCAL_SYM struct partition*
