commit aff0cf67939db10a2cb5261dc1eadf82c0b425a6
parent fd89ed39c48a441f035fedf1754007b38d4c1cd3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 22 Jul 2022 12:29:54 +0200
Start building the octrees
Diffstat:
8 files changed, 224 insertions(+), 8 deletions(-)
diff --git a/README.md b/README.md
@@ -21,9 +21,10 @@ 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) and
+[Star-Mesh](https://gitlab.com/meso-star/star-mesh)
[Star-UnstructuredVolumetricMesh](https://gitlab.com/meso-star/star-uvm)
-libraries. It optionally depends on [scdoc](https://sr.ht/~sircmpwn/scdoc/)
+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.
First ensure that CMake is installed on your system. Then install the RCMake
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -28,6 +28,7 @@ option(NO_TEST "Do not build tests" OFF)
################################################################################
# Check dependencies
################################################################################
+find_package(OpenMP 1.2 REQUIRED)
find_package(RCMake 0.4 REQUIRED)
find_package(RNSF REQUIRED)
find_package(RSys 0.9 REQUIRED)
@@ -62,6 +63,7 @@ set(RNATM_FILES_SRC
rnatm.c
rnatm_log.c
rnatm_mesh.c
+ rnatm_octree.c
rnatm_properties.c
rnatm_voxel_partition.c)
set(RNATM_FILES_INC
@@ -79,9 +81,11 @@ 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)
+target_link_libraries(rnatm RNSF RSys StarAerosol StarBuffer StarCK StarMesh StarUVM m)
set_target_properties(rnatm PROPERTIES
+ COMPILE_FLAGS "${OpenMP_C_FLAGS}"
+ LINK_FLAGS "${OpenMP_C_FLAGS}"
DEFINE_SYMBOL RNATM_SHARED_BUILD
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/src/rnatm.c b/src/rnatm.c
@@ -32,6 +32,8 @@
#include <rsys/cstr.h>
#include <rsys/mem_allocator.h>
+#include <omp.h>
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -83,6 +85,12 @@ check_rnatm_create_args(const struct rnatm_create_args* args)
if(res != RES_OK) return res;
}
+ /* Check miscalleneous arguments */
+ if(!args->name
+ || !args->grid_definition_hint
+ || !args->nthreads)
+ return RES_BAD_ARG;
+
return RES_OK;
}
@@ -93,6 +101,7 @@ create_rnatm
{
struct rnatm* atm = NULL;
struct mem_allocator* allocator = NULL;
+ int nthreads_max = 0;
res_T res = RES_OK;
if(!out_atm) { res = RES_BAD_ARG; goto error;}
@@ -115,7 +124,6 @@ create_rnatm
res = RES_MEM_ERR;
goto error;
}
-
ref_init(&atm->ref);
atm->allocator = allocator;
atm->verbose = args->verbose;
@@ -123,6 +131,10 @@ create_rnatm
gas_init(atm->allocator, &atm->gas);
darray_aerosol_init(atm->allocator, &atm->aerosols);
+ /* Setup the number of threads */
+ nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
+ atm->nthreads = MMIN((unsigned)nthreads_max, args->nthreads);
+
if(args->logger) {
atm->logger = args->logger;
} else {
@@ -131,7 +143,7 @@ create_rnatm
if(args->verbose) {
fprintf(stderr, MSG_ERROR_PREFIX
"Could not setup the default logger of the "
- "Rad-Net ATMopshere library.\n");
+ "Rad-Net ATMopshere library\n");
}
goto error;
}
@@ -189,8 +201,8 @@ rnatm_create
if(res != RES_OK) goto error;
res = setup_properties(atm, args);
if(res != RES_OK) goto error;
- /*res = setup_octrees(atm, args);
- if(res != RES_OK) goto error;*/
+ res = setup_octrees(atm, args);
+ if(res != RES_OK) goto error;
exit:
if(out_atm) *out_atm = atm;
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -83,10 +83,13 @@ struct rnatm_create_args {
size_t naerosols;
const char* name; /* Name of the atmosphere */
+ unsigned grid_definition_hint; /* Hint on the grid definition */
int precompute_normals; /* Pre-compute the tetrahedra normals */
struct logger* logger; /* NULL <=> use default logger */
struct mem_allocator* allocator; /* NULL <=> use default allocator */
+
+ unsigned nthreads; /* Hint on the number of threads to use */
int verbose; /* Verbosity level */
};
#define RNATM_CREATE_ARGS_DEFAULT__ { \
@@ -95,10 +98,13 @@ struct rnatm_create_args {
0, /* Number of aerosols */ \
"atmosphere", /* Name */ \
\
+ 65536, /* Hint on the grid definition */ \
0, /* Precompute tetrahedra normals */ \
\
NULL, /* Logger */ \
NULL, /* Allocator */ \
+ \
+ (unsigned)~0, /* #threads */ \
0 /* Verbosity level */ \
}
static const struct rnatm_create_args RNDGR_CREATE_ARGS_DEFAULT =
diff --git a/src/rnatm_c.h b/src/rnatm_c.h
@@ -138,7 +138,11 @@ struct rnatm {
struct darray_aerosol aerosols;
struct str name;
+ unsigned grid_definition[3];
+
+ unsigned nthreads;
int verbose;
+
struct logger* logger;
struct logger logger__;
struct mem_allocator* allocator;
@@ -155,4 +159,9 @@ setup_properties
(struct rnatm* atm,
const struct rnatm_create_args* args);
+extern LOCAL_SYM res_T
+setup_octrees
+ (struct rnatm* atm,
+ const struct rnatm_create_args* args);
+
#endif /* RNATM_C_H */
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -0,0 +1,180 @@
+/* Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Institut de Physique du Globe de Paris
+ * Copyright (C) 2022 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université de Reims Champagne-Ardenne
+ * Copyright (C) 2022 Université de Versaille Saint-Quentin
+ * Copyright (C) 2022 Université Paul Sabatier (contact@laplace.univ-tlse.fr)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* lround support */
+
+#include "rnatm_c.h"
+#include "rnatm_log.h"
+#include "rnatm_voxel_partition.h"
+
+#include <star/sck.h>
+#include <star/suvm.h>
+
+#include <rsys/cstr.h>
+#include <rsys/math.h>
+#include <rsys/rsys.h>
+
+#include <math.h> /* lround */
+#include <omp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static size_t
+band_get_quad_pt_count(const size_t iband, void* ctx)
+{
+ struct sck_band band;
+ struct sck* sck = ctx;
+ ASSERT(ctx);
+
+ SCK(get_band(sck, iband, &band));
+ return band.quad_pts_count;
+}
+
+static FINLINE unsigned
+round_pow2(const unsigned val)
+{
+ const unsigned next_pow2 = (unsigned)round_up_pow2(val);
+ if(next_pow2 - val <= next_pow2/4) {
+ return next_pow2;
+ } else {
+ return next_pow2/2;
+ }
+}
+
+static res_T
+compute_grid_definition(struct rnatm* atm, const struct rnatm_create_args* args)
+{
+ double low[3];
+ double upp[3];
+ double sz[3];
+ double sz_max;
+ double vxsz;
+ unsigned def[3];
+ int iaxis_max;
+ int iaxis_remain[2];
+ int i;
+ ASSERT(atm && args);
+
+ /* 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));
+ sz[0] = upp[0] - low[0];
+ sz[1] = upp[1] - low[1];
+ sz[2] = upp[2] - low[2];
+
+ /* Find the axis along which the atmosphere's AABB is greatest */
+ sz_max = -DBL_MAX;
+ FOR_EACH(i, 0, 3) {
+ if(sz[i] > sz_max) { iaxis_max = i; sz_max = sz[i]; }
+ }
+
+ /* Define the other axes */
+ iaxis_remain[0] = (iaxis_max + 1) % 3;
+ iaxis_remain[1] = (iaxis_max + 2) % 3;
+
+ /* Fix the definition along the maximum axis and calculate the size of a
+ * voxel along this axis */
+ def[iaxis_max] = round_pow2(args->grid_definition_hint);
+ vxsz = sz[iaxis_max] / def[iaxis_max];
+
+ /* Calculate the definitions along the remaining 2 axes. First calculates
+ * them assuming the voxels are cubic, then rounds them to their nearest
+ * power of two */
+ def[iaxis_remain[0]] = (unsigned)lround(sz[iaxis_remain[0]]/vxsz);
+ def[iaxis_remain[1]] = (unsigned)lround(sz[iaxis_remain[1]]/vxsz);
+ def[iaxis_remain[0]] = round_pow2(def[iaxis_remain[0]]);
+ def[iaxis_remain[1]] = round_pow2(def[iaxis_remain[1]]);
+
+ /* Setup grid definition */
+ atm->grid_definition[0] = def[0];
+ atm->grid_definition[1] = def[1];
+ atm->grid_definition[2] = def[2];
+
+ return RES_OK;
+}
+
+static res_T
+build_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
+ * octrees from the partitions and the threads that fill these partitions */
+ const size_t NPARTITIONS_PER_THREAD = 16;
+
+ struct pool_create_args pool_args = POOL_CREATE_ARGS_DEFAULT;
+ struct pool* pool = NULL;
+
+ res_T res = RES_OK;
+ ASSERT(atm && args);
+
+ /* Create the vortex partition pool */
+ pool_args.npartitions = NPARTITIONS_PER_THREAD * atm->nthreads;
+ 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.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",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ log_info(atm, "Partitions radiative properties (grid definition: %ux%ux%u)\n",
+ SPLIT3(atm->grid_definition));
+
+ /* Enable nested threads to allow multi-threading when calculating radiative
+ * properties */
+ omp_set_nested(1);
+
+ /* TODO build the octrees */
+ (void)args;
+
+exit:
+ if(pool) pool_ref_put(pool);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+setup_octrees(struct rnatm* atm, const struct rnatm_create_args* args)
+{
+ res_T res = RES_OK;
+ ASSERT(atm && args);
+
+ res = compute_grid_definition(atm, args);
+ if(res != RES_OK) goto error;
+ res = build_octrees(atm, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -407,6 +407,7 @@ setup_properties(struct rnatm* atm, const struct rnatm_create_args* args)
FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) {
struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i;
+
res = setup_aerosol_properties(atm, aerosol, args->aerosols+i);
if(res != RES_OK) goto error;
}
diff --git a/src/rnatm_voxel_partition.h b/src/rnatm_voxel_partition.h
@@ -26,6 +26,8 @@
#include <rsys/list.h>
#include <rsys/rsys.h>
+/* TODO remove this */
+#if 0
/* Definition of a partition along the 3 axes */
#define PARTITION_LOG2_DEFINITION 5 /* 5 */
#define PARTITION_DEFINITION BIT(PARTITION_LOG2_DEFINITION) /* 2^5 = 32 */
@@ -33,6 +35,7 @@
( PARTITION_DEFINITION \
* PARTITION_DEFINITION \
* PARTITION_DEFINITION)
+#endif
struct pool_create_args {
size_t npartitions; /* Number of partitions to preallocate */
@@ -44,7 +47,7 @@ struct pool_create_args {
struct mem_allocator* allocator; /* NULL <=> default allocator */
};
-#define POOL_CREATE_ARGS_DEFAULT__ {0, 0, NULL, NULL, 0, NULL}
+#define POOL_CREATE_ARGS_DEFAULT__ {0, 0, NULL, NULL, 32, NULL}
static const struct pool_create_args POOL_CREATE_ARGS_DEFAULT =
POOL_CREATE_ARGS_DEFAULT__;
