commit 6c7aa5ceb2fcb86d382e4e17e561ea28ad231b5c
parent e931a643301a4874ec132a8fe0b5bec2a46a5d72
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 25 Aug 2022 18:32:26 +0200
Consider aerosols when building octrees
Diffstat:
| M | src/rnatm_octree.c | | | 187 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------- |
1 file changed, 170 insertions(+), 17 deletions(-)
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -24,6 +24,7 @@
#include "rnatm_log.h"
#include "rnatm_voxel_partition.h"
+#include <star/sars.h>
#include <star/sck.h>
#include <star/suvm.h>
#include <star/svx.h>
@@ -32,6 +33,7 @@
#include <rsys/condition.h>
#include <rsys/cstr.h>
#include <rsys/double3.h>
+#include <rsys/float4.h>
#include <rsys/math.h>
#include <rsys/morton.h>
#include <rsys/mutex.h>
@@ -40,6 +42,7 @@
#include <math.h> /* lround */
#include <omp.h>
+#define GAS SIZE_MAX
#define VOXELIZE_MSG "voxelize atmosphere: %3d%%\r"
/* Structure used to synchronise the voxelization and the build threads */
@@ -138,6 +141,7 @@ struct voxelize_args {
const struct accel_struct* accel_structs;
size_t batch_size; /* #spectral items */
+
};
#define VOXELIZE_ARGS_NULL__ { \
NULL, /* Tetra ids */ \
@@ -448,7 +452,7 @@ compute_grid_definition(struct rnatm* atm, const struct rnatm_create_args* args)
}
static INLINE void
-setup_tetra_radcoefs
+setup_tetra_radcoefs_gas
(struct rnatm* atm,
const struct suvm_primitive* tetra,
const size_t iband,
@@ -490,6 +494,110 @@ setup_tetra_radcoefs
radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
}
+static void
+setup_tetra_radcoefs_aerosol
+ (struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t iband,
+ const size_t iquad_pt,
+ const struct aerosol* aerosol,
+ struct radcoefs* radcoefs)
+{
+ struct sck_band gas_band;
+ struct sars_band ars_band;
+ double gas_spectral_range[2]; /* In cm^-1 */
+ size_t ars_ibands[2];
+ ASSERT(atm && aerosol && tetra && radcoefs);
+
+ /* Look for the aerosol bands covered by the gas band */
+ SCK(get_band(atm->gas.ck, iband, &gas_band));
+ gas_spectral_range[0] = gas_band.lower;
+ gas_spectral_range[1] = gas_band.upper;
+ SARS(find_overlaped_bands(aerosol->sars, gas_spectral_range, ars_ibands));
+
+ /* TODO ask VE what to do if:
+ * - the spectral meshes have holes
+ * - no aerosol band is overlaid by the gas band
+ * - the aerosol bands are entirely included in the gas band
+ * - gas/an aerosol band is degenerated (i.e. lower == upper) */
+
+ /* Aerosol spectral data is not overlaid by the gas band */
+ if(ars_ibands[0] > ars_ibands[1]) FATAL("Not implemented yet!\n");
+
+ SARS(get_band(aerosol->sars, ars_ibands[0], &ars_band));
+
+ /* The gas band is included in an aerosol band: use gas spectral data */
+ if(ars_ibands[0] == ars_ibands[1]
+ && ars_band.lower <= gas_band.lower
+ && ars_band.upper >= gas_band.upper) {
+ setup_tetra_radcoefs_gas(atm, tetra, iband, iquad_pt, radcoefs);
+
+ /* The gas band overlaid N aerosol bands (N >= 1) */
+ } else {
+ float tau_ka[4] = {0, 0, 0, 0};
+ float tau_ks[4] = {0, 0, 0, 0};
+ float ka[4], ks[4], kext[4];
+ float rcp_gas_band_len;
+ size_t iars_band;
+
+ FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) {
+ double lambda_min;
+ double lambda_max;
+ float lambda_len;
+
+ SARS(get_band(aerosol->sars, iars_band, &ars_band));
+ lambda_min = MMIN(gas_band.lower, ars_band.lower);
+ lambda_max = MMAX(gas_band.upper, ars_band.upper);
+ lambda_len = (float)(lambda_max - lambda_min);
+
+ tau_ks[0] += ars_band.ks_list[tetra->indices[0]] * lambda_len;
+ tau_ks[1] += ars_band.ks_list[tetra->indices[1]] * lambda_len;
+ tau_ks[2] += ars_band.ks_list[tetra->indices[2]] * lambda_len;
+ tau_ks[3] += ars_band.ks_list[tetra->indices[3]] * lambda_len;
+
+ tau_ka[0] += ars_band.ka_list[tetra->indices[0]] * lambda_len;
+ tau_ka[1] += ars_band.ka_list[tetra->indices[1]] * lambda_len;
+ tau_ka[2] += ars_band.ka_list[tetra->indices[2]] * lambda_len;
+ tau_ka[3] += ars_band.ka_list[tetra->indices[3]] * lambda_len;
+ }
+
+ /* Compute the radiative coefficients of the tetrahedron */
+ ASSERT(gas_band.upper != gas_band.lower); /* FIXME is it possible? */
+ rcp_gas_band_len = 1.f / (float)(gas_band.upper - gas_band.lower);
+ f4_mulf(ks, tau_ks, rcp_gas_band_len);
+ f4_mulf(ka, tau_ka, rcp_gas_band_len);
+ f4_add(kext, ka, ks);
+
+ /* Compute the radiative coefficients range of the tetrahedron */
+ radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
+ radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
+ radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
+ radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
+ radcoefs->kext_min = MMIN(MMIN(kext[0], kext[1]), MMIN(kext[2], kext[3]));
+ radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
+ }
+}
+
+static INLINE void
+setup_tetra_radcoefs
+ (struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t iband,
+ const size_t iquad_pt,
+ const size_t cpnt,
+ struct radcoefs* radcoefs)
+{
+ ASSERT(atm && tetra && radcoefs);
+ ASSERT(cpnt == GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
+
+ if(cpnt == GAS) {
+ setup_tetra_radcoefs_gas(atm, tetra, iband, iquad_pt, radcoefs);
+ } else {
+ const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+cpnt;
+ setup_tetra_radcoefs_aerosol(atm, tetra, iband, iquad_pt, aerosol, radcoefs);
+ }
+}
+
static INLINE void
update_voxel
(struct rnatm* atm,
@@ -530,10 +638,22 @@ update_voxel
}
static res_T
-voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
+voxelize_tetra
+ (struct rnatm* atm,
+ const struct voxelize_args* args,
+ const size_t cpnt)
{
+ const struct suvm_volume* volume = NULL;
size_t i;
ASSERT(atm && check_voxelize_args(args) == RES_OK);
+ ASSERT(cpnt == GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
+
+ /* Fetch the component volumetric mesh */
+ if(cpnt == GAS) {
+ volume = atm->gas.volume;
+ } else {
+ volume = darray_aerosol_cdata_get(&atm->aerosols)[cpnt].volume;
+ }
partition_clear_voxels(args->part);
@@ -553,7 +673,7 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
enum suvm_intersection_type intersect;
/* Recover the tetrahedron and setup its polyhedron */
- SUVM(volume_get_primitive(atm->gas.volume, itetra, &tetra));
+ SUVM(volume_get_primitive(volume, itetra, &tetra));
SUVM(primitive_setup_polyhedron(&tetra, &poly));
ASSERT(poly.lower[0] <= args->part_upp[0]);
ASSERT(poly.lower[1] <= args->part_upp[1]);
@@ -586,7 +706,7 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
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);
+ setup_tetra_radcoefs(atm, &tetra, iband, iquad_pt, cpnt, radcoefs);
}
/* Iterate voxels intersected by the AABB of the polyedron */
@@ -622,31 +742,64 @@ voxelize_gas(struct rnatm* atm, const struct voxelize_args* args)
}
static res_T
-voxelize_partition(struct rnatm* atm, const struct voxelize_args* args)
+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 && check_voxelize_args(args) == RES_OK);
+ ASSERT(atm && args);
+ ASSERT(cpnt == GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
+
+ /* Get the volumetric mesh of the component */
+ if(cpnt == 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
- (atm->gas.volume,
- args->part_low,
- args->part_upp,
- register_tetra,
- args->tetra_ids);
+ (volume, args->part_low, args->part_upp, register_tetra, args->tetra_ids);
if(res != RES_OK) goto error;
- /* The partition is not covered by any tetrahedron */
- if(darray_size_t_size_get(args->tetra_ids) == 0) {
- partition_empty(args->part);
- goto exit;
+ if(darray_size_t_size_get(args->tetra_ids) != 0) {
+ res = voxelize_tetra(atm, args, cpnt);
+ if(res != RES_OK) goto error;
}
- res = voxelize_gas(atm, args);
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+voxelize_partition(struct rnatm* atm, const struct voxelize_args* args)
+{
+ size_t iaerosol = 0;
+ size_t naerosols = 0;
+ int part_is_empty = 1;
+ res_T res = RES_OK;
+ ASSERT(atm && check_voxelize_args(args) == RES_OK);
+
+ /* Voxelize the gas */
+ res = voxelize_partition_component(atm, args, GAS);
if(res != RES_OK) goto error;
+ part_is_empty = part_is_empty && !darray_size_t_size_get(args->tetra_ids);
+
+ /* Voxelize the aerosols */
+ naerosols = darray_aerosol_size_get(&atm->aerosols);
+ FOR_EACH(iaerosol, 0, naerosols) {
+ res = voxelize_partition_component(atm, args, iaerosol);
+ if(res != RES_OK) goto error;
+ part_is_empty = part_is_empty && !darray_size_t_size_get(args->tetra_ids);
+ }
- /* TODO voxelise aerosols */
+ /* The partition is not overlaped by any tetrahedron */
+ if(part_is_empty) partition_empty(args->part);
exit:
return res;
