commit 87db36ce6a0e77a970e5f0c7a797b1b1f2bd8672
parent 3e65e4e673fd92d40491c5fda7e3ee86a2d85ddb
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 22 Aug 2022 18:01:42 +0200
Preparing to build multiple octrees
Diffstat:
8 files changed, 245 insertions(+), 151 deletions(-)
diff --git a/src/rnatm.c b/src/rnatm.c
@@ -133,6 +133,7 @@ create_rnatm
str_init(atm->allocator, &atm->name);
gas_init(atm->allocator, &atm->gas);
darray_aerosol_init(atm->allocator, &atm->aerosols);
+ darray_accel_struct_init(atm->allocator, &atm->accel_structs);
/* Setup the number of threads */
nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
@@ -187,14 +188,13 @@ release_rnatm(ref_T* ref)
{
struct rnatm* atm = CONTAINER_OF(ref, struct rnatm, ref);
ASSERT(ref);
-
if(atm->logger == &atm->logger__) logger_release(&atm->logger__);
- if(atm->octree) SVX(tree_ref_put(atm->octree));
+ darray_aerosol_release(&atm->aerosols);
+ darray_accel_struct_release(&atm->accel_structs);
if(atm->svx_allocator_is_init) {
ASSERT(MEM_ALLOCATED_SIZE(&atm->svx_allocator) == 0);
mem_shutdown_regular_allocator(&atm->svx_allocator);
}
- darray_aerosol_release(&atm->aerosols);
gas_release(&atm->gas);
str_release(&atm->name);
MEM_RM(atm->allocator, atm);
@@ -428,3 +428,47 @@ aerosol_copy_and_release(struct aerosol* dst, struct aerosol* src)
src->sars = NULL;
return darray_phase_fn_copy_and_release(&dst->phase_fn_lst, &src->phase_fn_lst);
}
+
+res_T
+accel_struct_init
+ (struct mem_allocator* allocator,
+ struct accel_struct* accel_struct)
+{
+ (void)allocator;
+ ASSERT(accel_struct);
+ accel_struct->octree = NULL;
+ accel_struct->iband = 0;
+ accel_struct->iquad_pt = 0;
+ return RES_OK;
+}
+
+void
+accel_struct_release(struct accel_struct* accel_struct)
+{
+ ASSERT(accel_struct);
+ if(accel_struct->octree) SVX(tree_ref_put(accel_struct->octree));
+}
+
+res_T
+accel_struct_copy(struct accel_struct* dst, const struct accel_struct* src)
+{
+ ASSERT(dst && src);
+ if(dst->octree) SVX(tree_ref_put(dst->octree));
+ dst->octree = src->octree;
+ dst->iband = src->iband;
+ dst->iquad_pt = src->iquad_pt;
+ if(dst->octree) SVX(tree_ref_get(dst->octree));
+ return RES_OK;
+}
+
+res_T
+accel_struct_copy_and_release(struct accel_struct* dst, struct accel_struct* src)
+{
+ ASSERT(dst && src);
+ if(dst->octree) SVX(tree_ref_put(dst->octree));
+ dst->octree = src->octree;
+ dst->iband = src->iband;
+ dst->iquad_pt = src->iquad_pt;
+ src->octree = NULL;
+ return RES_OK;
+}
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -104,7 +104,7 @@ struct rnatm_create_args {
\
1, /* Optical thickness */ \
\
- 16384, /* Hint on the grid definition */ \
+ 512, /* Hint on the grid definition */ \
0, /* Precompute tetrahedra normals */ \
\
NULL, /* Logger */ \
@@ -155,6 +155,7 @@ rnatm_get_k_svx_voxel
RNATM_API res_T
rnatm_write_vtk_octree
(const struct rnatm* atm,
+ const size_t iaccel_struct,
FILE* stream);
#endif /* RNATM_H */
diff --git a/src/rnatm_c.h b/src/rnatm_c.h
@@ -140,17 +140,53 @@ aerosol_copy_and_release
#include <rsys/dynamic_array.h>
/*******************************************************************************
+ * Acceleration structure
+ ******************************************************************************/
+struct accel_struct {
+ struct svx_tree* octree;
+ size_t iband; /* Index of the spectral band */
+ size_t iquad_pt; /* Index of the quadrature point */
+};
+
+extern LOCAL_SYM res_T
+accel_struct_init
+ (struct mem_allocator* allocator,
+ struct accel_struct* accel_struct);
+
+extern LOCAL_SYM void
+accel_struct_release
+ (struct accel_struct* accel_struct);
+
+extern LOCAL_SYM res_T
+accel_struct_copy
+ (struct accel_struct* dst,
+ const struct accel_struct* src);
+
+extern LOCAL_SYM res_T
+accel_struct_copy_and_release
+ (struct accel_struct* dst,
+ struct accel_struct* src);
+
+/* Define the dynamic array of acceleration structures */
+#define DARRAY_NAME accel_struct
+#define DARRAY_DATA struct accel_struct
+#define DARRAY_FUNCTOR_INIT accel_struct_init
+#define DARRAY_FUNCTOR_RELEASE accel_struct_release
+#define DARRAY_FUNCTOR_COPY accel_struct_copy_and_release
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE accel_struct_copy_and_release
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
* Atmosphere
******************************************************************************/
struct rnatm {
struct gas gas;
struct darray_aerosol aerosols;
+ struct darray_accel_struct accel_structs;
struct str name;
unsigned grid_definition[3];
- struct svx_tree* octree;
-
unsigned nthreads;
int verbose;
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -38,16 +38,17 @@
#include <math.h> /* lround */
#include <omp.h>
-struct build_octrees_context {
+struct build_octree_context {
struct pool* pool;
struct partition* part; /* Current partition */
+ size_t iitem; /* Index of the item to handle in the partition */
/* Optical thickness threshold criteria for merging voxels */
double tau_threshold;
};
-#define BUILD_OCTREES_CONTEXT_NULL__ {NULL, NULL, 0}
-static const struct build_octrees_context BUILD_OCTREES_CONTEXT_NULL =
- BUILD_OCTREES_CONTEXT_NULL__;
+#define BUILD_OCTREE_CONTEXT_NULL__ {NULL, NULL, 0, 0}
+static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
+ BUILD_OCTREE_CONTEXT_NULL__;
struct radcoefs {
float ka_min, ka_max;
@@ -60,12 +61,6 @@ struct radcoefs {
FLT_MAX, -FLT_MAX \
}
-struct spectral_item {
- size_t iband;
- size_t iquad_pt;
-};
-#define SPECTRAL_ITEM_NULL__ {0, 0}
-
/* Generate the dynamic array of radiative coefficient lists */
#define DARRAY_NAME radcoef_list
#define DARRAY_DATA struct radcoefs*
@@ -87,8 +82,8 @@ struct voxelize_batch_args {
double vxsz[3]; /* Size of a voxel */
- struct spectral_item* items;
- size_t batch_size; /* #spectral items */
+ struct accel_struct* accel_structs;
+ size_t batch_size; /* #acceleration structures */
};
#define VOXELIZE_BATCH_ARGS_NULL__ { \
NULL, /* Per thread tetra list */ \
@@ -124,7 +119,7 @@ struct voxelize_args {
double vxsz[3]; /* Size of a voxel */
- const struct spectral_item* items;
+ const struct accel_struct* accel_structs;
size_t batch_size; /* #spectral items */
};
#define VOXELIZE_ARGS_NULL__ { \
@@ -170,7 +165,7 @@ check_voxelize_batch_args
|| !args->vxsz[0]
|| !args->vxsz[1]
|| !args->vxsz[2]
- || !args->items
+ || !args->accel_structs
|| !args->batch_size) {
return RES_BAD_ARG;
}
@@ -198,7 +193,7 @@ check_voxelize_args(const struct voxelize_args* args)
|| !args->vxsz[0]
|| !args->vxsz[1]
|| !args->vxsz[2]
- || !args->items
+ || !args->accel_structs
|| !args->batch_size) {
return RES_BAD_ARG;
}
@@ -217,7 +212,7 @@ round_pow2(const unsigned val)
}
/* Return the number of quadrature points for the submitted range of bands.
- * boundaries are inclusives */
+ * Lmits are inclusives */
static INLINE size_t
compute_nquad_pts(const struct rnatm* atm, const size_t bands[2])
{
@@ -236,29 +231,30 @@ compute_nquad_pts(const struct rnatm* atm, const size_t bands[2])
}
static res_T
-setup_batches
+setup_accel_structs
(struct rnatm* atm,
- const size_t bands[2], /* Boundaries are inclusives */
- struct spectral_item** out_items,
- size_t* out_nitems)
+ const size_t bands[2]) /* Range of bands to handle. Limits are inclusives */
{
- struct spectral_item* items = NULL;
- size_t iitem;
- size_t nitems;
- size_t nbands;
+ struct accel_struct* accel_structs = NULL;
+ size_t istruct;
+ size_t naccel_structs;
size_t i;
+ size_t nbands;
res_T res = RES_OK;
- ASSERT(atm && out_items && out_nitems && bands && bands[0] <= bands[1]);
+ ASSERT(atm && bands && bands[0] <= bands[1]);
+
+ naccel_structs = compute_nquad_pts(atm, bands);
- nitems = compute_nquad_pts(atm, bands);
- items = MEM_CALLOC(atm->allocator, nitems, sizeof(*items));
- if(items == NULL) {
- log_err(atm, "%s: failed to allocate memory\n", FUNC_NAME);
- res = RES_MEM_ERR;
+ res = darray_accel_struct_resize(&atm->accel_structs, naccel_structs);
+ if(res != RES_OK) {
+ log_err(atm,
+ "%s: failed to allocate the list of acceleration structures -- %s",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
+ accel_structs = darray_accel_struct_data_get(&atm->accel_structs);
- iitem = 0;
+ istruct = 0;
nbands = bands[1] - bands[0] + 1;
FOR_EACH(i, 0, nbands) {
struct sck_band band;
@@ -267,22 +263,16 @@ setup_batches
SCK(get_band(atm->gas.ck, iband, &band));
FOR_EACH(iquad_pt, 0, band.quad_pts_count) {
- ASSERT(iitem < nitems);
- items[iitem].iband = iband;
- items[iitem].iquad_pt = iquad_pt;
+ ASSERT(istruct < naccel_structs);
+ accel_structs[istruct].iband = iband;
+ accel_structs[istruct].iquad_pt = iquad_pt;
}
}
exit:
- *out_items = items;
- *out_nitems = nitems;
return res;
error:
- if(items) {
- MEM_RM(atm->allocator, items);
- items = NULL;
- nitems = 0;
- }
+ darray_accel_struct_clear(&atm->accel_structs);
goto exit;
}
@@ -404,7 +394,8 @@ update_voxel
(struct rnatm* atm,
const struct radcoefs* radcoefs,
struct partition* part,
- const uint64_t vx_mcode)
+ const uint64_t vx_mcode,
+ const uint64_t vx_iitem)
{
float vx_ka_min, vx_ka_max;
float vx_ks_min, vx_ks_max;
@@ -414,7 +405,7 @@ update_voxel
ASSERT(atm && radcoefs && part);
(void)atm;
- vx = partition_get_voxel(part, vx_mcode, 0/*TODO handle batch_size > 1*/);
+ vx = partition_get_voxel(part, vx_mcode, vx_iitem);
/* Update the range of the radiative coefficients of the voxel */
vx_ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
@@ -443,9 +434,6 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
size_t i;
ASSERT(atm && check_voxelize_args(args) == RES_OK);
- /* FIXME one has to support a batch_size > 1 */
- ASSERT(args->batch_size == 1);
-
partition_clear_voxels(args->part);
FOR_EACH(i, 0, darray_size_t_size_get(args->tetra_ids)) {
@@ -459,6 +447,7 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
uint32_t ivx_upp[3];
uint32_t ivx[3];
uint64_t mcode[3]; /* Cache of 3D morton code */
+ size_t iitem;
const size_t itetra = darray_size_t_cdata_get(args->tetra_ids)[i];
enum suvm_intersection_type intersect;
@@ -491,21 +480,13 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
ASSERT(ivx_upp[1] <= args->part_def);
ASSERT(ivx_upp[2] <= args->part_def);
- /* Compute the range of the tetrahedron radiative coefficients
- * FIXME handle batch_size > 1 */
- #if 1
- setup_tetra_radcoefs(atm, &tetra,
- args->items[0].iband,
- args->items[0].iquad_pt,
- &args->per_item_radcoefs[0]);
- #else
- FOR_EACH(item, 0, args->batch_size) {
- const size_t iband = args->items[item].iband;
- const size_t iquad_pt = args->items[item].iquad_pt;
- struct radcoefs* radcoefs = args->per_item_radcoefs + item;
+ /* Compute the range of the tetrahedron radiative coefficients */
+ FOR_EACH(iitem, 0, args->batch_size) {
+ const size_t iband = args->accel_structs[iitem].iband;
+ const size_t iquad_pt = args->accel_structs[iitem].iquad_pt;
+ struct radcoefs* radcoefs = args->per_item_radcoefs + iitem;
setup_tetra_radcoefs(atm, &tetra, iband, iquad_pt, radcoefs);
}
- #endif
/* Iterate voxels intersected by the AABB of the polyedron */
FOR_EACH(ivx[2], ivx_low[2], ivx_upp[2]) {
@@ -526,15 +507,11 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
intersect = suvm_polyhedron_intersect_aabb(&poly, vx_low, vx_upp);
if(intersect == SUVM_INTERSECT_NONE) continue;
- /* Update the intersected voxel FIXME handle batch_size > 1*/
- #if 1
- update_voxel(atm, &args->per_item_radcoefs[0], args->part, mcode[0]);
- #else
- FOR_EACH(item, 0, args->batch_size) {
- struct radcoefs* radcoefs = args->per_item_radcoefs + item;
- update_voxel(atm, radcoefs, part, mcode[0]);
+ /* 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);
}
- #endif
}
}
}
@@ -585,9 +562,6 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
ATOMIC res = RES_OK;
ASSERT(atm && check_voxelize_batch_args(atm, args) == RES_OK);
- /* FIXME one has to support a batch_size > 1 */
- ASSERT(args->batch_size == 1);
-
/* Print status message */
#define LOG_MSG "voxelize atmosphere: %3d%%\r"
log_info(atm, LOG_MSG, 0);
@@ -611,9 +585,7 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
if(ATOMIC_GET(&res) != RES_OK) continue;
- /* Recover the batch partitions.
- * FIXME add support of multiple spectral item
- * per partition to allow several voxels the _same_ partition id */
+ /* Recover the batch partitions */
part = pool_next_partition(args->pool);
if(!part) { ATOMIC_SET(&res, RES_UNKNOWN_ERR); continue; };
@@ -654,7 +626,7 @@ voxelize_batch(struct rnatm* atm, const struct voxelize_batch_args* args)
d3_set(voxelize_args.vxsz, args->vxsz);
voxelize_args.part = part;
voxelize_args.per_item_radcoefs = per_item_radcoefs;
- voxelize_args.items = args->items;
+ voxelize_args.accel_structs = args->accel_structs;
voxelize_args.batch_size = args->batch_size;
res_local = voxelize_partition(atm, &voxelize_args);
if(res_local == RES_OK) {
@@ -691,22 +663,18 @@ voxelize_atmosphere
(struct rnatm* atm,
/* Range of spectral bands to handle. Limits are inclusives */
const size_t bands[2],
- struct pool* pool,
- const size_t batch_size)
+ struct pool* pool)
{
/* Batch of spectral items to process in a single voxelization */
struct voxelize_batch_args batch_args = VOXELIZE_BATCH_ARGS_NULL;
- size_t nbatches;
- size_t ibatch;
+ size_t batch_size = 0;
+ size_t nbatches = 0;
+ size_t ibatch = 0;
/* Working data structures */
struct darray_size_t_list per_thread_tetra_list;
struct darray_radcoef_list per_thread_radcoef_list;
- /* List of spectral items to handle */
- struct spectral_item* items = NULL;
- size_t nitems = 0;
-
/* Voxelization parameters */
double atm_low[3];
double atm_upp[3];
@@ -716,16 +684,15 @@ voxelize_atmosphere
size_t part_def; /* Definition of a partition */
/* Miscellaneous variables */
- size_t i;
+ struct accel_struct* accel_structs = NULL;
+ size_t naccel_structs = 0;
+ size_t i = 0;
res_T res = RES_OK;
ASSERT(atm && bands && pool);
ASSERT(bands[0] <= bands[1]);
ASSERT(bands[1] < sck_get_bands_count(atm->gas.ck));
- /* FIXME handle batch_size > 1 */
- ASSERT(batch_size == 1);
-
darray_size_t_list_init(atm->allocator, &per_thread_tetra_list);
darray_radcoef_list_init(atm->allocator, &per_thread_radcoef_list);
@@ -736,6 +703,7 @@ voxelize_atmosphere
if(res != RES_OK) goto error;
/* Setup the per thread and per item working data structures */
+ batch_size = pool_get_voxel_width(pool);
FOR_EACH(i, 0, atm->nthreads) {
struct radcoefs* per_item_k = NULL;
per_item_k = MEM_CALLOC(atm->allocator, batch_size, sizeof(*per_item_k));
@@ -743,9 +711,6 @@ voxelize_atmosphere
darray_radcoef_list_data_get(&per_thread_radcoef_list)[i] = per_item_k;
}
- res = setup_batches(atm, bands, &items, &nitems);
- if(res != RES_OK) goto error;
-
/* 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];
@@ -765,13 +730,14 @@ voxelize_atmosphere
/* Calculate the number of batches required to build the accelerating
* structures for the spectral bands submitted; currently we are building one
- * octree per quadrature point of each band (i.e. per spectral item). And
- * each octree requires a complete voxelization of the atmospheric meshes
- * which takes time. To amortize this cost without prohitive increase in
- * memory space, one fills up to N octrees from a single voxelization of the
- * mesh. The number of batches corresponds to the total number of
- * voxelization required */
- nbatches = (nitems + (batch_size - 1)/*ceil*/) / batch_size;
+ * octree per quadrature point of each band. And each octree requires a
+ * complete voxelization of the atmospheric meshes which takes time. To
+ * amortize this cost without prohitive increase in memory space, one fills
+ * up to N octrees from a single voxelization of the mesh. The number of
+ * batches corresponds to the total number of voxelization required */
+ accel_structs = darray_accel_struct_data_get(&atm->accel_structs);
+ naccel_structs = darray_accel_struct_size_get(&atm->accel_structs);
+ nbatches = (naccel_structs + (batch_size - 1)/*ceil*/) / batch_size;
/* Print the size of a voxel */
log_info(atm, "voxel size = {%g, %g, %g}\n", SPLIT3(vxsz));
@@ -793,9 +759,9 @@ voxelize_atmosphere
FOR_EACH(ibatch, 0, nbatches) {
size_t item_range[2];
item_range[0] = ibatch * batch_size;
- item_range[1] = MMIN(item_range[0] + batch_size, nitems);
+ item_range[1] = MMIN(item_range[0] + batch_size, naccel_structs);
- batch_args.items = items + item_range[0];
+ batch_args.accel_structs = accel_structs + item_range[0];
batch_args.batch_size = item_range[1] - item_range[0];
res = voxelize_batch(atm, &batch_args);
if(res != RES_OK) goto error;
@@ -806,7 +772,6 @@ exit:
MEM_RM(atm->allocator,
darray_radcoef_list_data_get(&per_thread_radcoef_list)[i]);
}
- MEM_RM(atm->allocator, items);
darray_size_t_list_release(&per_thread_tetra_list);
darray_radcoef_list_release(&per_thread_radcoef_list);
return res;
@@ -817,7 +782,7 @@ error:
static void
vx_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
{
- struct build_octrees_context* ctx = context;
+ struct build_octree_context* ctx = context;
const float* vx = NULL;
uint64_t ivx, ipart;
int log2_part_def;
@@ -845,7 +810,7 @@ vx_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
}
}
- vx = partition_cget_voxel(ctx->part, ivx, 0/*TODO handle batch_size > 1*/);
+ vx = partition_cget_voxel(ctx->part, ivx, ctx->iitem);
memcpy(dst, vx, NFLOATS_PER_VOXEL * sizeof(float));
}
@@ -883,7 +848,7 @@ vx_challenge_merge
const size_t nvxs,
void* context)
{
- const struct build_octrees_context* ctx = context;
+ const struct build_octree_context* ctx = context;
double low[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
double upp[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
double tau;
@@ -927,59 +892,81 @@ static res_T
build_octrees
(struct rnatm* atm,
const struct rnatm_create_args* args,
+ /* Range of spectral bands to handle. Limits are inclusives */
+ const size_t bands[2],
struct pool* pool)
{
- struct build_octrees_context ctx = BUILD_OCTREES_CONTEXT_NULL;
- struct svx_voxel_desc vx_desc = SVX_VOXEL_DESC_NULL;
struct svx_device* svx = NULL;
double low[3], upp[3];
size_t def[3];
- res_T res = RES_OK;
- ASSERT(atm && args && pool);
+ size_t istruct;
+ size_t naccel_structs;
+ ATOMIC res = RES_OK;
+ ASSERT(atm && args && bands && pool);
+ ASSERT(bands[0] <= bands[1]);
+ ASSERT(bands[1] < sck_get_bands_count(atm->gas.ck));
res = svx_device_create(atm->logger, &atm->svx_allocator, atm->verbose, &svx);
if(res != RES_OK) goto error;
- /* Setup the build context */
- ctx.pool = pool;
- ctx.tau_threshold = args->optical_thickness;
- ctx.part = NULL;
-
- /* Setup the voxel descriptor */
- vx_desc.get = vx_get;
- vx_desc.merge = vx_merge;
- vx_desc.challenge_merge = vx_challenge_merge;
- vx_desc.context = &ctx;
- vx_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
-
/* Recover the AABB from the atmosphere. Note that we have already made sure
* when setting up the meshes of the gas and aerosols that the aerosols are
* included in the gas. Therefore, the AABB of the gas is the same as the
* AABB of the atmosphere */
SUVM(volume_get_aabb(atm->gas.volume, low, upp));
- /* Build the octree. TODO currently only one octree is built while we have to
- * build as many octrees as quadrature points */
+ /* Retrieve grid definition */
def[0] = (size_t)atm->grid_definition[0];
def[1] = (size_t)atm->grid_definition[1];
def[2] = (size_t)atm->grid_definition[2];
- res = svx_octree_create(svx, low, upp, def, &vx_desc, &atm->octree);
+
+ naccel_structs = darray_accel_struct_size_get(&atm->accel_structs);
+
+ /* Build the octrees. Each thread consumes an element of a partition. We
+ * therefore set the number of threads to the width of a voxel of the
+ * partition */
+ omp_set_num_threads((int)pool_get_voxel_width(pool));
+ #pragma omp parallel for schedule(static, 1/*chunk size*/)
+ for(istruct = 0; istruct < naccel_structs; ++istruct) {
+ struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
+ struct svx_voxel_desc vx_desc = SVX_VOXEL_DESC_NULL;
+ struct svx_tree* octree = NULL;
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ /* Setup the build context */
+ ctx.pool = pool;
+ ctx.part = NULL;
+ ctx.iitem = istruct % pool_get_voxel_width(pool);
+ ctx.tau_threshold = args->optical_thickness;
+
+ /* Setup the voxel descriptor */
+ vx_desc.get = vx_get;
+ vx_desc.merge = vx_merge;
+ vx_desc.challenge_merge = vx_challenge_merge;
+ vx_desc.context = &ctx;
+ vx_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
+
+ res_local = svx_octree_create(svx, low, upp, def, &vx_desc, &octree);
+ if(ctx.part) partition_free(ctx.part);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; };
+
+ /* Register the built octree */
+ darray_accel_struct_data_get(&atm->accel_structs)[istruct].octree = octree;
+ }
if(res != RES_OK) goto error;
exit:
- if(ctx.part) partition_free(ctx.part);
if(svx) SVX(device_ref_put(svx));
- return res;
+ return (res_T)res;
error:
- if(atm->octree) {
- SVX(tree_ref_put(atm->octree));
- atm->octree = NULL;
- }
+ darray_accel_struct_clear(&atm->accel_structs);
goto exit;
}
static res_T
-create_pool(struct rnatm* atm, struct pool** out_pool)
+create_pool(struct rnatm* atm, struct pool** out_pool, const size_t nbands)
{
/* Empirical constant that defines the number of non empty partitions to
* pre-allocate per thread to avoid contension between the thread building the
@@ -993,6 +980,12 @@ create_pool(struct rnatm* atm, struct pool** out_pool)
* with a sparse grid without significantly increasing memory usage */
const size_t OVERALL_NPARTITIONS_PER_THREAD = NPARTITIONS_PER_THREAD * 64;
+ /* Number of items stored per voxel data. A width greater than 1 makes it
+ * possible to store by partition the radiative coefficients of several
+ * spectral data. The goal is to voxelize the volumetric mesh once for N
+ * spectral data */
+ const size_t VOXEL_WIDTH = MMIN(4, nbands);
+
/* Definition of a partition on the 3 axis */
const size_t PARTITION_DEFINITION = 8;
@@ -1001,11 +994,11 @@ create_pool(struct rnatm* atm, struct pool** out_pool)
res_T res = RES_OK;
ASSERT(atm && out_pool);
- /* Create the partition pool
- * TODO add support of multiple spectral item per pool */
+ /* Create the partition pool */
pool_args.npartitions = atm->nthreads * OVERALL_NPARTITIONS_PER_THREAD;
pool_args.npreallocated_partitions = atm->nthreads * NPARTITIONS_PER_THREAD;
pool_args.partition_definition = PARTITION_DEFINITION;
+ pool_args.voxel_width = VOXEL_WIDTH;
pool_args.allocator = atm->allocator;
res = pool_create(&pool_args, &pool);
if(res != RES_OK) goto error;
@@ -1024,29 +1017,30 @@ static res_T
create_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
{
struct pool* pool = NULL;
+ const size_t bands[2] = {0,0}; /* TODO handle multiple bands */
ATOMIC res = RES_OK;
ASSERT(atm);
- res = create_pool(atm, &pool);
+ res = create_pool(atm, &pool, bands[1] - bands[0] + 1);
+ if(res != RES_OK) goto error;
+ res = setup_accel_structs(atm, bands);
if(res != RES_OK) goto error;
log_info(atm, "partitionning of radiative properties (grid definition: %ux%ux%u)\n",
SPLIT3(atm->grid_definition));
- /* Enable nested threads to allow multi-threading when calculating radiative
- * properties */
+ /* Enable nested threads to allow multi-threading when voxelizing tetrahedra
+ * and building octrees */
omp_set_nested(1);
#pragma omp parallel sections num_threads(2)
{
#pragma omp section
{
- const size_t bands[2] = {0,0}; /* TODO handle multiple bands */
- const res_T res_local = voxelize_atmosphere
- (atm, bands, pool, 1/*FIXME handle batch_size > 1*/);
+ const res_T res_local = voxelize_atmosphere(atm, bands, pool);
if(res_local != RES_OK) {
- log_err(atm, "atmosphere voxelization error -- %s\n",
+ log_err(atm, "error when voxelizing the atmopshere -- %s\n",
res_to_cstr((res_T)res));
pool_invalidate(pool);
ATOMIC_SET(&res, res_local);
@@ -1055,7 +1049,7 @@ create_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
#pragma omp section
{
- const res_T res_local = build_octrees(atm, args, pool);
+ const res_T res_local = build_octrees(atm, args, bands, pool);
if(res_local != RES_OK) {
log_err(atm, "error building octrees -- %s\n", res_to_cstr((res_T)res));
pool_invalidate(pool);
@@ -1107,4 +1101,3 @@ exit:
error:
goto exit;
}
-
diff --git a/src/rnatm_voxel_partition.c b/src/rnatm_voxel_partition.c
@@ -461,6 +461,13 @@ pool_get_partition_definition(const struct pool* pool)
return pool->partition_definition;
}
+size_t
+pool_get_voxel_width(const struct pool* pool)
+{
+ ASSERT(pool);
+ return pool->voxel_width;
+}
+
struct partition*
pool_next_partition(struct pool* pool)
{
diff --git a/src/rnatm_voxel_partition.h b/src/rnatm_voxel_partition.h
@@ -111,6 +111,10 @@ extern LOCAL_SYM size_t
pool_get_partition_definition
(const struct pool* pool);
+extern LOCAL_SYM size_t
+pool_get_voxel_width
+ (const struct pool* pool);
+
/* Returns a free partition. Waits for a free partition to be available.
* Returns NULL if an error occurs */
extern LOCAL_SYM struct partition*
diff --git a/src/rnatm_write_vtk_octree.c b/src/rnatm_write_vtk_octree.c
@@ -138,18 +138,24 @@ register_leaf
* Exported function
******************************************************************************/
res_T
-rnatm_write_vtk_octree(const struct rnatm* atm, FILE* stream)
+rnatm_write_vtk_octree
+ (const struct rnatm* atm,
+ const size_t iaccel_struct,
+ FILE* stream)
{
char buf[128];
struct time t0, t1;
struct octree_data data;
- const double* leaf_data;
+ struct svx_tree* octree = NULL;
+ const double* leaf_data = NULL;
size_t nvertices;
size_t ncells;
size_t i;
res_T res = RES_OK;
- if(!atm || !stream) {
+ if(!atm
+ || !stream
+ || iaccel_struct >= darray_accel_struct_size_get(&atm->accel_structs)) {
res = RES_BAD_ARG;
goto error;
}
@@ -158,16 +164,18 @@ rnatm_write_vtk_octree(const struct rnatm* atm, FILE* stream)
log_info(atm, "write VTK octree\n");
time_current(&t0);
+ octree = darray_accel_struct_cdata_get(&atm->accel_structs)[iaccel_struct].octree;
+
octree_data_init(atm, &data);
/* Register leaf data */
- SVX(tree_for_each_leaf(atm->octree, register_leaf, &data));
+ SVX(tree_for_each_leaf(octree, register_leaf, &data));
nvertices = darray_double_size_get(&data.vertices)/3/*#coords per vertex*/;
ncells = darray_size_t_size_get(&data.cells)/8/*#ids per cell*/;
#ifndef NDEBUG
{
struct svx_tree_desc octree_desc = SVX_TREE_DESC_NULL;
- SVX(tree_get_desc(atm->octree, &octree_desc));
+ SVX(tree_get_desc(octree, &octree_desc));
ASSERT(ncells == octree_desc.nleaves);
}
#endif
diff --git a/src/test_rnatm.c b/src/test_rnatm.c
@@ -254,7 +254,8 @@ write_vtk_octree(struct rnatm* atm, const char* filename)
goto error;
}
- res = rnatm_write_vtk_octree(atm, fp);
+ res = rnatm_write_vtk_octree
+ (atm, 0/*FIXME handle a user defined octree id*/, fp);
if(res != RES_OK) goto error;
exit:
