commit dd0cdda823b4cfd81f2fe8704b48ca7034e1ce6d
parent 638f2fff672dc9e39577379ecada3612975a0dc1
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 14 Jan 2021 14:09:48 +0100
Finalize the "build octree" process
Diffstat:
9 files changed, 418 insertions(+), 43 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -23,6 +23,7 @@ option(NO_TEST "Do not build tests" OFF)
################################################################################
# Check dependencies
################################################################################
+find_package(AtrRI 0.0 REQUIRED)
find_package(AtrTP 0.0 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
find_package(RCMake 0.4 REQUIRED)
@@ -36,10 +37,12 @@ include(rcmake)
include(rcmake_runtime)
include_directories(
+ ${AtrRI_INCLUDE_DIR}
${AtrTP_INCLUDE_DIR}
${RSys_INCLUDE_DIR}
${StarTetraHedra_INCLUDE_DIR}
- ${StarUVM_INCLUDE_DIR})
+ ${StarUVM_INCLUDE_DIR}
+ ${StarVX_INCLUDE_DIR})
################################################################################
# Configure and define targets
@@ -53,12 +56,14 @@ set(ATRSTM_FILES_SRC
atrstm.c
atrstm_build_octrees.c
atrstm_log.c
+ atrstm_optprops.c
atrstm_partition.c
atrstm_setup_uvm.c)
set(ATRSTM_FILES_INC
atrstm_build_octrees.h
atrstm_c.h
atrstm_log.h
+ atrstm_optprops.h
atrstm_partition.h)
set(ATRSTM_FILES_INC_API
atrstm.h)
@@ -72,7 +77,7 @@ rcmake_prepend_path(ATRSTM_FILES_INC_API ${ATRSTM_SOURCE_DIR})
rcmake_prepend_path(ATRSTM_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(atrstm SHARED ${ATRSTM_FILES_SRC} ${ATRSTM_FILES_INC} ${ATRSTM_FILES_INC_API})
-target_link_libraries(atrstm AtrTP RSys StarTetraHedra StarUVM StarVX)
+target_link_libraries(atrstm AtrRI AtrTP RSys StarTetraHedra StarUVM StarVX m)
if(CMAKE_COMPILER_IS_GNUCC)
set_target_properties(atrstm PROPERTIES LINK_FLAGS "${OpenMP_C_FLAGS}")
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -17,6 +17,7 @@
#include "atrstm_c.h"
#include "atrstm_log.h"
+#include <astoria/atrri.h>
#include <astoria/atrtp.h>
#include <rsys/mem_allocator.h>
@@ -53,6 +54,7 @@ release_atrstm(ref_T* ref)
ASSERT(ref);
atrstm = CONTAINER_OF(ref, struct atrstm, ref);
if(atrstm->atrtp) ATRTP(ref_put(atrstm->atrtp));
+ if(atrstm->atrri) ATRRI(ref_put(atrstm->atrri));
if(atrstm->suvm) SUVM(device_ref_put(atrstm->suvm));
if(atrstm->svx) SVX(device_ref_put(atrstm->svx));
if(atrstm->logger == &atrstm->logger__) logger_release(&atrstm->logger__);
@@ -100,6 +102,11 @@ atrstm_create
atrstm->allocator = allocator;
atrstm->verbose = args->verbose;
atrstm->nthreads = MMIN(args->nthreads, (unsigned)nthreads_max);
+ atrstm->gyration_radius_prefactor = args->gyration_radius_prefactor;
+ atrstm->fractal_dimension = args->fractal_dimension;
+ atrstm->spectral_type = args->spectral_type;
+ atrstm->wlen_range[0] = args->wlen_range[0];
+ atrstm->wlen_range[1] = args->wlen_range[1];
if(logger) {
atrstm->logger = logger;
@@ -124,6 +131,14 @@ atrstm_create
goto error;
}
+ /* Load the refractive index */
+ res = atrri_create
+ (atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->atrri);
+ if(res != RES_OK) goto error;
+ res = atrri_load(atrstm->atrri, args->atrri_filename);
+ if(res != RES_OK) goto error;
+
+ /* Create the Star-UnstructuredVolumetricMesh device */
res = suvm_device_create
(atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->suvm);
if(res != RES_OK) goto error;
@@ -143,6 +158,7 @@ atrstm_create
/* Create the Star-VoXel device */
res = svx_device_create
(atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->svx);
+ if(res != RES_OK) goto error;
exit:
if(out_atrstm) *out_atrstm = atrstm;
diff --git a/src/atrstm.h b/src/atrstm.h
@@ -47,8 +47,8 @@ enum atrstm_property {
/* List of medium components */
enum atrstm_component {
- ATRSTM_CPNT_GAS,
ATRSTM_CPNT_SOOT,
+ ATRSTM_CPNT_GAS,
ATRSTM_CPNTS_COUNT__
};
@@ -70,11 +70,18 @@ struct atrstm_args {
const char* atrtp_filename; /* Filename of the Thermodynamic properties */
const char* atrri_filename; /* Filename of the refraction index LUT */
const char* name; /* Name of the gas mixture */
+
enum atrstm_spectral_type spectral_type; /* Longwave/shortwave */
double wlen_range[2]; /* Spectral range to handle In nm */
+
+ double gyration_radius_prefactor;
+ double fractal_dimension;
+
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
double optical_thickness; /* Threshold used during octree building */
+
int precompute_normals; /* Pre-compute the tetrahedra normals */
+
unsigned nthreads; /* Hint on the number of threads to use */
int verbose; /* Verbosity level */
};
@@ -87,6 +94,8 @@ struct atrstm_args {
"gas mixture", /* Name */ \
ATRSTM_SPECTRAL_TYPES_COUNT__, /* Spectral type */ \
{DBL_MAX,-DBL_MAX}, /* Spectral integration range */ \
+ 1.70, /* Prefactor */ \
+ 1.75, /* Fractal dimension */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Acceleration grid max definition */ \
1, /* Optical thickness */ \
1, /* Precompute tetrahedra normals */ \
diff --git a/src/atrstm_build_octrees.c b/src/atrstm_build_octrees.c
@@ -13,11 +13,16 @@
* 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 /* nextafterf */
+
#include "atrstm_c.h"
#include "atrstm_build_octrees.h"
#include "atrstm_log.h"
-#include "atrstm_svx.h"
+#include "atrstm_optprops.h"
#include "atrstm_partition.h"
+#include "atrstm_svx.h"
+
+#include <astoria/atrri.h>
#include <star/suvm.h>
#include <star/svx.h>
@@ -26,6 +31,7 @@
#include <rsys/dynamic_array_size_t.h>
#include <rsys/morton.h>
+#include <math.h>
#include <omp.h>
struct build_octree_context {
@@ -78,6 +84,60 @@ register_primitive
CHK(darray_size_t_push_back(prims, &prim->iprim) == RES_OK);
}
+static void
+voxel_commit_optprops_range
+ (float* vox,
+ const enum atrstm_component cpnt,
+ const struct optprops* optprops_min,
+ const struct optprops* optprops_max)
+{
+ double ka[2], ks[2], kext[2];
+ float vox_ka[2], vox_ks[2], vox_kext[2];
+ ASSERT(vox);
+ ASSERT(optprops_min);
+ ASSERT(optprops_max);
+
+ ka[0] = optprops_min->ka;
+ ka[1] = optprops_max->ka;
+ ks[0] = optprops_min->ks;
+ ks[1] = optprops_max->ks;
+ kext[0] = optprops_min->kext;
+ kext[1] = optprops_max->kext;
+
+ /* Ensure that the single precision bounds include their double precision
+ * version. */
+ if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX);
+ if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX);
+ if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX);
+ if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX);
+ if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX);
+ if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
+
+ /* Fetch the range of the voxel optical properties */
+ vox_ka[0] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN);
+ vox_ka[1] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX);
+ vox_ks[0] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN);
+ vox_ks[1] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX);
+ vox_kext[0] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN);
+ vox_kext[1] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX);
+
+ /* Update the range of the voxel optical properties */
+ vox_ka[0] = MMIN(vox_ka[0], (float)ka[0]);
+ vox_ka[1] = MMAX(vox_ka[1], (float)ka[1]);
+ vox_ks[0] = MMIN(vox_ks[0], (float)ks[0]);
+ vox_ks[1] = MMAX(vox_ks[1], (float)ks[1]);
+ vox_kext[0] = MMIN(vox_kext[0], (float)kext[0]);
+ vox_kext[1] = MMAX(vox_kext[1], (float)kext[1]);
+
+ /* Commit the upd to the voxel */
+ vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, vox_ka[0]);
+ vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, vox_ka[1]);
+ vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, vox_ks[0]);
+ vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, vox_ks[1]);
+ vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, vox_kext[0]);
+ vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, vox_kext[1]);
+}
+
static res_T
voxelize_partition
(struct atrstm* atrstm,
@@ -85,16 +145,18 @@ voxelize_partition
const double part_low[3], /* Lower bound of the partition */
const double part_upp[3], /* Upper bound of the partition */
const double vxsz[3], /* Size of a partition voxel */
+ const struct atrri_refractive_index* refract_id, /* Refractive index */
struct part* part)
{
size_t iprim;
res_T res = RES_OK;
- ASSERT(atrstm && prims && part_low && part_upp && vxsz && part);
+ ASSERT(atrstm && prims && part_low && part_upp && vxsz && part && refract_id);
ASSERT(part_low[0] < part_upp[0]);
ASSERT(part_low[1] < part_upp[1]);
ASSERT(part_low[2] < part_upp[2]);
ASSERT(vxsz[0] > 0 && vxsz[1] > 0 && vxsz[2] > 0);
+ /* Clean the partition voxels */
part_clear_voxels(part);
FOR_EACH(iprim, 0, darray_size_t_size_get(prims)) {
@@ -154,23 +216,33 @@ voxelize_partition
intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
if(intersect != SUVM_INTERSECT_NONE) {
+ struct optprops optprops_min;
+ struct optprops optprops_max;
float* vox = part_get_voxel(part, mcode[2]);
- /* TODO Register the optical properties of the primitive against
- * the current voxel. */
- (void)vox;
+
+ res = primitive_compute_optprops_range
+ (atrstm, refract_id, &prim, &optprops_min, &optprops_max);
+ if(UNLIKELY(res != RES_OK)) goto error;
+
+ voxel_commit_optprops_range
+ (vox, ATRSTM_CPNT_SOOT, &optprops_min, &optprops_max);
}
}
}
}
}
+exit:
return res;
+error:
+ goto exit;
}
static res_T
voxelize_volumetric_mesh
(struct atrstm* atrstm,
const struct build_octrees_args* args,
+ const struct atrri_refractive_index* refract_id,
struct pool* pools)
{
struct darray_prims_list prims_list; /* Per thread list of primitives */
@@ -256,8 +328,9 @@ voxelize_volumetric_mesh
if(!part) { /* An error occurs */
ATOMIC_SET(&res, res_local);
continue;
- }
- res = voxelize_partition(atrstm, prims, part_low, part_upp, vxsz, part);
+ }
+ res = voxelize_partition
+ (atrstm, prims, part_low, part_upp, vxsz, refract_id, part);
if(res != RES_OK) {
pool_free_partition(pool, part);
ATOMIC_SET(&res, res_local);
@@ -279,8 +352,8 @@ voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
const struct build_octree_context* ctx = context;
struct pool* pool = NULL;
struct part* part = NULL;
- float* vxl = NULL;
- uint64_t ivxl, ipart, ipool;
+ float* vox = NULL;
+ uint64_t ivox, ipart, ipool;
ASSERT(xyz && dst && ctx);
(void)xyz;
@@ -289,7 +362,7 @@ voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
* morton order and thus the partition morton ID is encoded in the MSB of the
* morton code while the voxels morton ID is stored in its LSB. */
ipart = (mcode >> (LOG2_PARTITION_DEFINITION*3));
- ivxl = (mcode & BIT_U64(LOG2_PARTITION_DEFINITION*3));
+ ivox = (mcode & BIT_U64(LOG2_PARTITION_DEFINITION*3));
/* Compute the pool index containing the partition. Partitions are
* alternatively stored into the per thread pool. Consequentlu the i^th
@@ -305,13 +378,13 @@ voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
if(!part) { /* An error occurs */
memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
} else {
- vxl = part_get_voxel(part, ivxl); /* Fetch the voxel data */
+ vox = part_get_voxel(part, ivox); /* Fetch the voxel data */
/* Copy voxel data */
- memcpy(dst, vxl, NFLOATS_PER_VOXEL * sizeof(float));
+ memcpy(dst, vox, NFLOATS_PER_VOXEL * sizeof(float));
/* Free the partition if the fetch voxel was its last one */
- if(ivxl == PARTITION_NVOXELS - 1) {
+ if(ivox == PARTITION_NVOXELS - 1) {
pool_free_partition(pool, part);
}
}
@@ -328,7 +401,7 @@ voxels_merge_component
float ks[2] = {FLT_MAX,-FLT_MAX};
float kext[2] = {FLT_MAX,-FLT_MAX};
size_t ivox;
- ASSERT(dst && voxs && nvoxs);
+ ASSERT(dst && voxels && nvoxs);
FOR_EACH(ivox, 0, nvoxs) {
const float* vox = voxels[ivox];
@@ -413,14 +486,20 @@ voxels_challenge_merge
/* Compute the Kext range for each component */
voxels_component_compute_kext_range(ATRSTM_CPNT_GAS, voxels, nvoxs, kext_gas);
voxels_component_compute_kext_range(ATRSTM_CPNT_SOOT, voxels, nvoxs, kext_soot);
- ASSERT(kext_gas[1] >= kext_gas[0]);
- ASSERT(kext_soot[1] >= kext_soot[0]);
-
- /* Check if the voxels are candidates to the merge process by evaluating its
- * optical thickness regarding the maximum size of their AABB and a user
- * defined threshold */
- merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold;
- merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold;
+
+ if(kext_gas[0] > kext_gas[1]) { /* Empty voxels */
+ /* Always merge empty voxels */
+ merge_gas = 1;
+ merge_soot = 1;
+ } else {
+ /* Check if the voxels are candidates to the merge process by evaluating its
+ * optical thickness regarding the maximum size of their AABB and a user
+ * defined threshold */
+ ASSERT(kext_gas[0] <= kext_gas[1]);
+ ASSERT(kext_soot[0] <= kext_soot[1]);
+ merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold;
+ merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold;
+ }
return merge_gas && merge_soot;
}
@@ -443,7 +522,8 @@ build_octree
ctx.pools = pools;
ctx.tau_threshold = args->optical_thickness_threshold;
- /* Setup the voxel descriptor */
+ /* Setup the voxel descriptor. TODO in shortwave, the NFLOATS_PER_VOXEL
+ * could be divided by 2 since there is no gas to handle */
vox_desc.get = voxel_get;
vox_desc.merge = voxels_merge;
vox_desc.challenge_merge = voxels_challenge_merge;
@@ -472,11 +552,23 @@ error:
res_T
build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
{
+ struct atrri_refractive_index refract_id = ATRRI_REFRACTIVE_INDEX_NULL;
struct pool* pools = NULL;
size_t i;
+ double wlen;
ATOMIC res = RES_OK;
ASSERT(atrstm && check_build_octrees_args(args));
+ /* Currently only shortwave is handled and in this case, radiative
+ * transfert is monochromatic */
+ ASSERT(atrstm->spectral_type == ATRSTM_SPECTRAL_SW);
+ ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]);
+ wlen = atrstm->wlen_range[0];
+
+ /* Fetch the submitted wavelength */
+ res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &refract_id);
+ if(res != RES_OK) goto error;
+
/* Allocate the per thread pool of partition */
pools = MEM_CALLOC(atrstm->allocator, atrstm->nthreads, sizeof(*pools));
if(!pools) { res = RES_MEM_ERR; goto error; }
@@ -490,7 +582,8 @@ build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
{
#pragma omp section
{
- res_T res_local = voxelize_volumetric_mesh(atrstm, args, pools);
+ const res_T res_local = voxelize_volumetric_mesh
+ (atrstm, args, &refract_id, pools);
if(res_local != RES_OK) {
size_t ipool;
log_err(atrstm, "Error voxelizing the volumetric mesh -- %s\n",
@@ -504,7 +597,7 @@ build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
#pragma omp section
{
- res_T res_local = build_octree(atrstm, args, pools);
+ const res_T res_local = build_octree(atrstm, args, pools);
if(res_local != RES_OK) {
size_t ipool;
log_err(atrstm, "Error buildin the octree -- %s\n",
@@ -516,6 +609,7 @@ build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
}
}
}
+ if(res != RES_OK) goto error;
exit:
if(pools) {
diff --git a/src/atrstm_c.h b/src/atrstm_c.h
@@ -16,17 +16,21 @@
#ifndef ATRSTM_C_H
#define ATRSTM_C_H
+#include "atrstm.h"
+
#include <rsys/logger.h>
#include <rsys/ref_count.h>
#include <rsys/str.h>
/* Forward declaration of external data types */
+struct atrri;
struct atrtp;
struct mem_allocator;
struct suvm;
struct atrstm {
- struct atrtp* atrtp; /* Handle toward the Astoria Thermodynamic properties */
+ struct atrtp* atrtp; /* Handle toward the Astoria Thermodynamic Properties */
+ struct atrri* atrri; /* Handle toward the Astoria Refractive Indices */
/* Unstructured volumetric mesh */
struct suvm_device* suvm;
@@ -39,6 +43,12 @@ struct atrstm {
struct str name; /* Name of the gas mixture */
+ double gyration_radius_prefactor;
+ double fractal_dimension;
+
+ enum atrstm_spectral_type spectral_type;
+ double wlen_range[2]; /* Spectral range in nm */
+
struct mem_allocator* allocator;
struct logger* logger;
struct logger logger__; /* Default logger */
diff --git a/src/atrstm_optprops.c b/src/atrstm_optprops.c
@@ -0,0 +1,74 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#include "atrstm_c.h"
+#include "atrstm_optprops.h"
+
+#include <astoria/atrri.h>
+#include <astoria/atrtp.h>
+
+#include <star/suvm.h>
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+primitive_compute_optprops
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct optprops optprops[])
+{
+ struct optprops_compute_args args = OPTPROPS_COMPUTE_ARGS_NULL;
+ struct atrtp_desc desc = ATRTP_DESC_NULL;
+ size_t inode;
+ res_T res = RES_OK;
+ ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && optprops);
+
+ res = atrtp_get_desc(atrstm->atrtp, &desc);
+ if(res != RES_OK) goto error;
+
+ /* Setup the constant input arguments of the optical properties computation
+ * routine */
+ args.lambda = refract_id->wavelength;
+ args.n = refract_id->n;
+ args.kappa = refract_id->kappa;
+ args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_dimension = atrstm->fractal_dimension;
+
+ FOR_EACH(inode, 0, 4) {
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]);
+
+ /* Setup the per node input arguments of the optical properties computation
+ * routine */
+ args.soot_volumic_fraction =
+ node[ATRTP_SOOT_VOLFRAC];
+ args.soot_primary_particles_count =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ args.soot_primary_particles_diameter =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+
+ /* Compute the node optical properties */
+ optprops_compute(&optprops[inode], &args);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/atrstm_optprops.h b/src/atrstm_optprops.h
@@ -0,0 +1,150 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#ifndef ATRSTM_OPTPROPS_H
+#define ATRSTM_OPTPROPS_H
+
+#include <rsys/math.h>
+#include <rsys/rsys.h>
+
+/* Optical properties */
+struct optprops {
+ double ka; /* Absorption coefficient */
+ double ks; /* Scattering coefficient */
+ double kext; /* Extinction coefficient */
+};
+#define OPTPROPS_NULL__ {0}
+static const struct optprops OPTPROPS_NULL = OPTPROPS_NULL__;
+
+struct optprops_compute_args {
+ double lambda; /* wavelength [nm] */
+ double n; /* Refractive index (real component) */
+ double kappa; /* Refractive index (imaginary component */
+ double gyration_radius_prefactor;
+ double fractal_dimension;
+ double soot_volumic_fraction; /* In m^3(soot) / m^3 */
+ double soot_primary_particles_count;
+ double soot_primary_particles_diameter; /* In nm */
+};
+#define OPTPROPS_COMPUTE_ARGS_NULL__ {0}
+static const struct optprops_compute_args OPTPROPS_COMPUTE_ARGS_NULL =
+ OPTPROPS_COMPUTE_ARGS_NULL__;
+
+/* Forward declarations */
+struct atrri_refractive_index;
+struct atrstm;
+struct suvm_primitive;
+
+extern LOCAL_SYM res_T
+primitive_compute_optprops
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct optprops optprops[]); /* Per primitive node optical properties */
+
+static INLINE res_T
+primitive_compute_optprops_range
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct optprops* optprops_min,
+ struct optprops* optprops_max)
+{
+ struct optprops props[4];
+ res_T res = RES_OK;
+ ASSERT(optprops_min && optprops_max);
+
+ res = primitive_compute_optprops(atrstm, refract_id, prim, props);
+ if(res != RES_OK) return res;
+
+ optprops_min->ka = MMIN
+ (MMIN(props[0].ka, props[1].ka),
+ MMIN(props[2].ka, props[3].ka));
+ optprops_min->ks = MMIN
+ (MMIN(props[0].ks, props[1].ks),
+ MMIN(props[2].ks, props[3].ks));
+ optprops_min->kext = MMIN
+ (MMIN(props[0].kext, props[1].kext),
+ MMIN(props[2].kext, props[3].kext));
+
+ optprops_max->ka = MMAX
+ (MMAX(props[0].ka, props[1].ka),
+ MMAX(props[2].ka, props[3].ka));
+ optprops_max->ks = MMAX
+ (MMAX(props[0].ks, props[1].ks),
+ MMAX(props[2].ks, props[3].ks));
+ optprops_max->kext = MMAX
+ (MMAX(props[0].kext, props[1].kext),
+ MMAX(props[2].kext, props[3].kext));
+
+ return RES_OK;
+}
+
+/* Compute the optical properties of the semi transparent medium as described
+ * in "Effects of multiple scattering on radiative properties of soot fractal
+ * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and
+ * Radiative Transfer Vol 133, pp. 374-381, 2014 */
+static INLINE void
+optprops_compute
+ (struct optprops* optprops,
+ const struct optprops_compute_args* args)
+{
+ /* Input variables */
+ const double lambda = args->lambda; /* [nm] */
+ const double n = args->n;
+ const double kappa = args->kappa;
+ const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
+ const double Np = args->soot_primary_particles_count;
+ const double Dp = args->soot_primary_particles_diameter; /* [nm] */
+ const double kf = args->gyration_radius_prefactor;
+ const double Df = args->fractal_dimension;
+
+ /* Precomputed values */
+ const double n2 = n*n;
+ const double kappa2 = kappa*kappa;
+ const double four_n2_kappa2 = 4*n2*kappa2;
+ const double xp = (PI*Dp) / lambda;
+ const double xp3 = xp*xp*xp;
+ const double k = (2*PI) / lambda;
+ const double k2 = k*k;
+
+ /* E(m), m = n + i*kappa */
+ const double tmp0 = (n2 - kappa2 + 2);
+ const double denom = tmp0*tmp0 + four_n2_kappa2;
+ const double Em = (6*n*kappa) / denom;
+
+ /* F(m), m = n + i*kappa */
+ const double tmp1 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom;
+ const double Fm = tmp1*tmp1 + Em*Em;
+
+ /* Gyration radius */
+ const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
+ const double Rg2 = Rg*Rg;
+ const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
+
+ /* Cross sections, [nm^3/aggrefate] */
+ const double Cabs = Np * (4*PI*xp3)/(k*k) * Em;
+ const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g;
+
+ /* Optical properties */
+ const double Va = (PI * Dp*Dp*Dp) / 6; /* [nm^3] */
+ const double rho = Fv / Va; /* [#aggregates / nm^3] */
+ const double tmp2 = 1.e-9 * rho;
+ optprops->ka = tmp2 * Cabs;
+ optprops->ks = tmp2 * Csca;
+ optprops->kext = optprops->ka + optprops->kext;
+}
+
+#endif /* ATRSTM_OPTPROPS_H */
diff --git a/src/atrstm_partition.h b/src/atrstm_partition.h
@@ -21,7 +21,7 @@
#include <string.h>
/* Definition of a partition along the 4 axis */
-#define LOG2_PARTITION_DEFINITION 5
+#define LOG2_PARTITION_DEFINITION 5
#define PARTITION_DEFINITION BIT(LOG2_PARTITION_DEFINITION) /*32*/
#define PARTITION_NVOXELS \
( PARTITION_DEFINITION \
@@ -59,13 +59,6 @@ struct pool {
ATOMIC error;
};
-static FINLINE void
-part_clear_voxels(struct part* part)
-{
- ASSERT(part);
- memset(part, 0, sizeof(part->voxels));
-}
-
static FINLINE float*
part_get_voxel
(struct part* part,
@@ -76,6 +69,15 @@ part_get_voxel
return part->voxels + ivoxel * NFLOATS_PER_VOXEL;
}
+static FINLINE void
+part_clear_voxels(struct part* part)
+{
+ size_t ivox;
+ FOR_EACH(ivox, 0, PARTITION_NVOXELS) {
+ vox_clear(part_get_voxel(part, ivox));
+ }
+}
+
extern LOCAL_SYM res_T
pool_init
(struct mem_allocator* allocator, /* May be NULL <=> default allocator */
diff --git a/src/atrstm_svx.h b/src/atrstm_svx.h
@@ -46,31 +46,46 @@
static FINLINE float
vox_get
- (const float* data,
+ (const float* vox,
const enum atrstm_component cpnt,
const enum atrstm_property prop,
const enum atrstm_svx_op op)
{
- ASSERT(data);
+ ASSERT(vox);
ASSERT((unsigned)cpnt < ATRSTM_CPNTS_COUNT__);
ASSERT((unsigned)prop < ATRSTM_PROPS_COUNT__);
ASSERT((unsigned)op < ATRSTM_SVX_OPS_COUNT__);
- return data[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op];
+ return vox[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op];
}
static FINLINE void
vox_set
- (float* data,
+ (float* vox,
const enum atrstm_component cpnt,
const enum atrstm_property prop,
const enum atrstm_svx_op op,
const float val)
{
- ASSERT(data);
+ ASSERT(vox);
ASSERT((unsigned)cpnt < ATRSTM_CPNTS_COUNT__);
ASSERT((unsigned)prop < ATRSTM_PROPS_COUNT__);
ASSERT((unsigned)op < ATRSTM_SVX_OPS_COUNT__);
- data[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op] = val;
+ vox[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op] = val;
+}
+
+static FINLINE void
+vox_clear(float* vox)
+{
+ int cpnt;
+ ASSERT(vox);
+
+ FOR_EACH(cpnt, 0, ATRSTM_CPNTS_COUNT__) {
+ int prop;
+ FOR_EACH(prop, 0, ATRSTM_PROPS_COUNT__) {
+ vox_set(vox, cpnt, prop, ATRSTM_SVX_MIN, FLT_MAX);
+ vox_set(vox, cpnt, prop, ATRSTM_SVX_MAX,-FLT_MAX);
+ }
+ }
}
#endif /* ATRSTM_SVX_H */
