commit e272f9c46e39be32a57cc1c3dd1587f68d766fbd
parent 487293b5bdf192f326f24dd538d9f00cb6c5b896
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 25 Jul 2022 17:40:36 +0200
Build the acceleration structure
Diffstat:
6 files changed, 260 insertions(+), 19 deletions(-)
diff --git a/README.md b/README.md
@@ -21,8 +21,9 @@ on the [Rad-Net Scattering Function](https://gitlab.com/meso-star/rnsf),
[Star-Aerosol](https://gitlab.com/meso-star/star-aerosol),
[Star-Buffer](https://gitlab.com/meso-star/star-buffer),
[Star-CorrelatedK](https://gitlab.com/meso-star/star-ck),
-[Star-Mesh](https://gitlab.com/meso-star/star-mesh)
-[Star-UnstructuredVolumetricMesh](https://gitlab.com/meso-star/star-uvm)
+[Star-Mesh](https://gitlab.com/meso-star/star-mesh),
+[Star-UnstructuredVolumetricMesh](https://gitlab.com/meso-star/star-uvm),
+[Star-VoXel](https://gitlab.com/meso-star/star-vx)
libraries, and on [OpenMP](https://www.openmp.org) 1.2 the parallelize its
computations. It optionally depends on [scdoc](https://sr.ht/~sircmpwn/scdoc/)
which, if available, is used to generate the man pages.
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -37,6 +37,7 @@ find_package(StarBuffer REQUIRED)
find_package(StarCK REQUIRED)
find_package(StarMesh REQUIRED)
find_package(StarUVM 0.0 REQUIRED)
+find_package(StarVX 0.2 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -49,7 +50,8 @@ include_directories(
${StarBuffer_INCLUDE_DIR}
${StarCK_INCLUDE_DIR}
${StarMesh_INCLUDE_DIR}
- ${StarUVM_INCLUDE_DIR})
+ ${StarUVM_INCLUDE_DIR}
+ ${StarVX_INCLUDE_DIR})
################################################################################
# Configure and define targets
@@ -81,7 +83,8 @@ rcmake_prepend_path(RNATM_FILES_INC_API ${RNATM_SOURCE_DIR})
rcmake_prepend_path(RNATM_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(rnatm SHARED ${RNATM_FILES_SRC} ${RNATM_FILES_INC} ${RNATM_FILES_INC_API})
-target_link_libraries(rnatm RNSF RSys StarAerosol StarBuffer StarCK StarMesh StarUVM m)
+target_link_libraries(rnatm RNSF RSys StarAerosol StarBuffer StarCK StarMesh
+ StarUVM StarVX m)
set_target_properties(rnatm PROPERTIES
COMPILE_FLAGS "${OpenMP_C_FLAGS}"
diff --git a/src/rnatm.c b/src/rnatm.c
@@ -28,6 +28,7 @@
#include <star/sbuf.h>
#include <star/sck.h>
#include <star/suvm.h>
+#include <star/svx.h>
#include <rsys/cstr.h>
#include <rsys/mem_allocator.h>
@@ -87,6 +88,7 @@ check_rnatm_create_args(const struct rnatm_create_args* args)
/* Check miscalleneous arguments */
if(!args->name
+ || args->optical_thickness < 0
|| !args->grid_definition_hint
|| !args->nthreads)
return RES_BAD_ARG;
@@ -113,7 +115,7 @@ create_rnatm
if(!atm) {
if(args->verbose) {
#define ERR_STR \
- "Could not allocate the device of the Rad-Net ATMosphere library"
+ "could not allocate the device of the Rad-Net ATMosphere library"
if(args->logger) {
logger_print(args->logger, LOG_ERROR, ERR_STR);
} else {
@@ -142,7 +144,7 @@ create_rnatm
if(res != RES_OK) {
if(args->verbose) {
fprintf(stderr, MSG_ERROR_PREFIX
- "Could not setup the default logger of the "
+ "could not setup the default logger of the "
"Rad-Net ATMopshere library\n");
}
goto error;
@@ -151,17 +153,26 @@ create_rnatm
res = darray_aerosol_resize(&atm->aerosols, args->naerosols);
if(res != RES_OK) {
- log_err(atm, "Could not allocate aerosol list -- %s\n", res_to_cstr(res));
+ log_err(atm, "could not allocate aerosol list -- %s\n", res_to_cstr(res));
goto error;
}
res = str_set(&atm->name, args->name);
if(res != RES_OK) {
- log_err(atm, "Could not setup the atmosphere name to `%s' -- %s\n",
+ log_err(atm, "could not setup the atmosphere name to `%s' -- %s\n",
args->name, res_to_cstr(res));
goto error;
}
+ res = mem_init_regular_allocator(&atm->svx_allocator);
+ if(res != RES_OK) {
+ log_err(atm,
+ "unable to initialize the allocator used to manage the memory of the "
+ "Star-VoXel library -- %s\n", res_to_cstr(res));
+ goto error;
+ }
+ atm->svx_allocator_is_init = 1;
+
exit:
if(out_atm) *out_atm = atm;
return res;
@@ -177,6 +188,11 @@ release_rnatm(ref_T* ref)
ASSERT(ref);
if(atm->logger == &atm->logger__) logger_release(&atm->logger__);
+ if(atm->octree) SVX(tree_ref_put(atm->octree));
+ 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);
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -84,6 +84,8 @@ struct rnatm_create_args {
size_t naerosols;
char* name; /* Name of the atmosphere */
+ double optical_thickness; /* Threshold used during octree building */
+
unsigned grid_definition_hint; /* Hint on the grid definition */
int precompute_normals; /* Pre-compute the tetrahedra normals */
@@ -99,6 +101,8 @@ struct rnatm_create_args {
0, /* Number of aerosols */ \
"atmosphere", /* Name */ \
\
+ 1, /* Optical thickness */ \
+ \
16384, /* Hint on the grid definition */ \
0, /* Precompute tetrahedra normals */ \
\
diff --git a/src/rnatm_c.h b/src/rnatm_c.h
@@ -149,12 +149,16 @@ struct rnatm {
unsigned grid_definition[3];
+ struct svx_tree* octree;
+
unsigned nthreads;
int verbose;
struct logger* logger;
struct logger logger__;
struct mem_allocator* allocator;
+ struct mem_allocator svx_allocator;
+ int svx_allocator_is_init;
ref_T ref;
};
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -26,7 +26,9 @@
#include <star/sck.h>
#include <star/suvm.h>
+#include <star/svx.h>
+#include <rsys/clock_time.h>
#include <rsys/cstr.h>
#include <rsys/math.h>
#include <rsys/morton.h>
@@ -35,6 +37,18 @@
#include <math.h> /* lround */
#include <omp.h>
+struct build_octrees_context {
+ struct rnatm* atm;
+ struct pool* pool;
+ struct partition* part; /* Current partition */
+
+ /* Optical thickness threshold criteria for merging voxels */
+ double tau_threshold;
+};
+#define BUILD_OCTREES_CONTEXT_NULL__ { NULL, NULL, NULL, 0 }
+static const struct build_octrees_context BUILD_OCTREES_CONTEXT_NULL =
+ BUILD_OCTREES_CONTEXT_NULL__;
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -145,12 +159,14 @@ update_voxel
ASSERT(tetra->nvertices == 4);
nbands = sck_get_bands_count(atm->gas.ck);
+ nbands = 1; /* FIXME For tests only */
#define REDUX_MIN(V4) MMIN(MMIN((V4)[0], (V4)[1]), MMIN((V4)[2], V4[3]))
#define REDUX_MAX(V4) MMAX(MMAX((V4)[0], (V4)[1]), MMAX((V4)[2], V4[3]))
FOR_EACH(iband, 0, nbands) {
struct sck_band band;
+ size_t nquad_pts;
float tetra_ks[4];
float tetra_ks_min, tetra_ks_max;
@@ -163,7 +179,10 @@ update_voxel
tetra_ks_min = REDUX_MIN(tetra_ks);
tetra_ks_max = REDUX_MAX(tetra_ks);
- FOR_EACH(iquad_pt, 0, band.quad_pts_count) {
+ nquad_pts = band.quad_pts_count;
+ nquad_pts = 1; /* FIXME For tests only */
+
+ FOR_EACH(iquad_pt, 0, nquad_pts) {
struct sck_quad_pt quad_pt;
float tetra_ka[4];
float tetra_kext[4];
@@ -381,7 +400,7 @@ voxelize_atmosphere(struct rnatm* atm, struct pool* pool)
nparts_adjusted = nparts_adjusted * nparts_adjusted * nparts_adjusted;
/* Print status message */
- #define LOG_MSG "Voxelization of atmosphere meshes: %3d%%\r"
+ #define LOG_MSG "voxelize atmosphere: %3d%%\r"
log_info(atm, LOG_MSG, 0);
/* Iterates over the partitions of the grid according to their Morton order
@@ -459,8 +478,182 @@ error:
goto exit;
}
+static void
+vx_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
+{
+ struct build_octrees_context* ctx = context;
+ float* vx = NULL;
+ uint64_t ivx, ipart;
+ float kext_min, kext_max;
+ int log2_part_def;
+ ASSERT(xyz && dst && ctx);
+ (void)xyz;
+
+ /* Retrieve the partition's morton id the voxel's morton id in the partition
+ * from the morton code of the voxel in the whole octree. Note that, like
+ * voxels, partitions are sorted in morton order. Therefore, the partition's
+ * morton id is encoded in the most significant bits of the voxel's morton
+ * code, while the voxel's morton id in the partition is stored in its least
+ * significant bits */
+ log2_part_def = log2i((int)pool_get_partition_definition(ctx->pool));
+ ipart = (mcode >> (log2_part_def*3));
+ ivx = (mcode & (BIT_U64(log2_part_def*3)-1));
+
+ /* Recover the partition storing the voxel */
+ if(ctx->part == NULL || partition_get_id(ctx->part) != ipart) {
+ if(ctx->part) pool_free_partition(ctx->pool, ctx->part);
+ ctx->part = pool_fetch_partition(ctx->pool, ipart);
+
+ if(ctx->part == NULL) { /* An error occurs */
+ memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
+ return;
+ }
+ }
+
+ vx = partition_get_voxel(ctx->part, 0, 0, ivx);
+ kext_min = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)];
+ kext_max = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)];
+
+ if(kext_min > kext_max) {
+ /* The voxel is not covered by any tetrahedron. Its coefficients are zero */
+ memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
+ } else {
+ memcpy(dst, vx, NFLOATS_PER_VOXEL * sizeof(float));
+ }
+}
+
+static void
+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);
+ (void)ctx;
+
+ FOR_EACH(ivx, 0, nvxs) {
+ const float* vx = vxs[ivx];
+ ka_min = MMIN(ka_min, vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)]);
+ 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
+vx_challenge_merge
+ (const struct svx_voxel vxs[],
+ const size_t nvxs,
+ void* context)
+{
+ const struct build_octrees_context* ctx = context;
+ double low[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
+ double upp[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
+ double tau;
+ float kext_min = FLT_MAX, kext_max = -FLT_MAX;
+ double sz_max;
+ size_t ivx;
+ ASSERT(vxs && nvxs && context);
+
+ /* 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)]);
+ }
+
+ /* Always merge empty voxels */
+ if(kext_min > kext_max) return 1;
+
+ /* Compute the AABB of the submitted voxels */
+ FOR_EACH(ivx, 0, nvxs) {
+ low[0] = MMIN(vxs[ivx].lower[0], low[0]);
+ low[1] = MMIN(vxs[ivx].lower[1], low[1]);
+ low[2] = MMIN(vxs[ivx].lower[2], low[2]);
+ upp[0] = MMAX(vxs[ivx].upper[0], upp[0]);
+ upp[1] = MMAX(vxs[ivx].upper[1], upp[1]);
+ upp[2] = MMAX(vxs[ivx].upper[2], upp[2]);
+ }
+
+ /* Compute the size of the largest dimension of the ABB */
+ sz_max = upp[0] - low[0];
+ sz_max = MMAX(upp[1] - low[1], sz_max);
+ sz_max = MMAX(upp[2] - low[2], sz_max);
+
+ /* Check if voxels can be merged by comparing the optical thickness along the
+ * maximum size of the merged AABB against a user-defined threshold */
+ tau = (kext_max - kext_min) * sz_max;
+ return tau < ctx->tau_threshold;
+}
+
+static res_T
+build_octrees
+ (struct rnatm* atm,
+ const struct rnatm_create_args* args,
+ 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);
+
+ res = svx_device_create(atm->logger, &atm->svx_allocator, atm->verbose, &svx);
+ if(res != RES_OK) goto error;
+
+ /* Setup the build context */
+ ctx.atm = atm;
+ 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 */
+ 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);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(ctx.part) pool_free_partition(pool, ctx.part);
+ if(svx) SVX(device_ref_put(svx));
+ return res;
+error:
+ if(atm->octree) {
+ SVX(tree_ref_put(atm->octree));
+ atm->octree = NULL;
+ }
+ goto exit;
+}
+
static res_T
-build_octrees(struct rnatm* atm)
+create_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
{
/* Empirical constant that defines the number of voxel partitions to
* pre-allocate per thread to avoid contension between the thread building the
@@ -478,16 +671,16 @@ build_octrees(struct rnatm* atm)
pool_args.nbands = sck_get_bands_count(atm->gas.ck);
pool_args.nquad_pts = band_get_quad_pt_count;
pool_args.context = atm->gas.ck;
- pool_args.partition_definition = 32;
+ pool_args.partition_definition = 4;
pool_args.allocator = atm->allocator;
res = pool_create(&pool_args, &pool);
if(res != RES_OK) {
- log_err(atm, "Failed to create the voxel partition pool -- %s\n",
+ log_err(atm, "failed to create the voxel partition pool -- %s\n",
res_to_cstr((res_T)res));
goto error;
}
- log_info(atm, "Partitions radiative properties (grid definition: %ux%ux%u)\n",
+ 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
@@ -500,8 +693,8 @@ build_octrees(struct rnatm* atm)
{
const res_T res_local = voxelize_atmosphere(atm, pool);
if(res_local != RES_OK) {
- log_err(atm, "Atmosphere voxelization error -- %s\n",
- res_to_cstr(res_local));
+ log_err(atm, "atmosphere voxelization error -- %s\n",
+ res_to_cstr((res_T)res));
pool_invalidate(pool);
ATOMIC_SET(&res, res_local);
}
@@ -509,8 +702,12 @@ build_octrees(struct rnatm* atm)
#pragma omp section
{
- /* TODO Build the octreese */
- FATAL("Not fully implemented yet\n");
+ const res_T res_local = build_octrees(atm, args, pool);
+ if(res_local != RES_OK) {
+ log_err(atm, "error building octrees -- %s\n", res_to_cstr((res_T)res));
+ pool_invalidate(pool);
+ ATOMIC_SET(&res, res_local);
+ }
}
}
if(res != RES_OK) goto error;
@@ -528,14 +725,30 @@ error:
res_T
setup_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
{
+ char buf[128];
+ struct time t0, t1;
+ size_t sz;
res_T res = RES_OK;
ASSERT(atm && args);
+ /* Start time recording */
+ time_current(&t0);
+
res = compute_grid_definition(atm, args);
if(res != RES_OK) goto error;
- res = build_octrees(atm);
+ res = create_octrees(atm, args);
if(res != RES_OK) goto error;
+ /* Log elapsed time */
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(atm, "acceleration structures built in %s\n", buf);
+
+ /* Log memory space used by Star-VX
+ * TODO make human readable the printed size */
+ sz = MEM_ALLOCATED_SIZE(&atm->svx_allocator);
+ log_info(atm, "Star-VoXel memory space: %lu\n", (unsigned long)sz);
+
exit:
return res;
error:
