commit 653e080817221a13132dc6a7b3f3256b1f95987a
parent fbcc6acaec441675842417f397999cf8f458d3ec
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 20 May 2021 14:43:51 +0200
Add flux map computation to htrdr-combustion
Diffstat:
8 files changed, 260 insertions(+), 35 deletions(-)
diff --git a/doc/htrdr-combustion.1.txt.in b/doc/htrdr-combustion.1.txt.in
@@ -198,6 +198,26 @@ OPTIONS
Path to the *atrtp*(5) file that stores the thermodynamic properties of the
combustion medium.
+*-R* <__rectangle-parameter__:...>::
+ Switch in flux map computation. The rectangular sensor onto which
+ the flux is integrated is defined by the following parameters:
+
+ **pos**=**_x_**,**_y_**,**_z_**;;
+ Position of the center of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_POS@}.
+
+ **tgt**=**_x_**,**_y_**,**_z_**;;
+ Position targeted by the rectangle, i.e. *tgt* - *pos* is the rectangle
+ normal. By default it is set to {@HTRDR_ARGS_DEFAULT_RECTANGLE_TGT@}.
+
+ **up**=**_x_**,**_y_**,**_z_**;;
+ Up vector of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_UP@}.
+
+ **sz**=**_width_**,**_height_**;;
+ Size of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_SZ@}.
+
*-r* _refract_ids_::
Path the the *atrri*(5) file that lists the spectrally varying refractive
indices of the combustion medium.
@@ -296,11 +316,12 @@ structure.
NOTES
-----
-1. Effects of multiple scattering on radiative properties of soot fractal
-aggregates. J. Yon et al, JQSRT 133, 374-381, 2014.
+1. Effects of multiple scattering on radiative properties of soot
+fractal aggregates. J. Yon et al, JQSRT 133, 374-381, 2014 -
+<https://doi.org/10.1016/j.jqsrt.2013.08.022>
2. Null-collision meshless Monte-Carlo - a new reverse Monte-Carlo algorithm
-designed for laser-source emission in absorbing/scattering inhomogeneous media. M.
+designed for laser-source emission in absorbing/scattering inhomogeneous media. M.
Sans et al, JQSRT, 2021 - <https://doi.org/10.1016/j.jqsrt.2021.107725>
3. VTK file format -
@@ -311,7 +332,7 @@ Sans et al, JQSRT, 2021 - <https://doi.org/10.1016/j.jqsrt.2021.107725>
COPYRIGHT
---------
Copyright © 2018, 2019, 2020, 2021 |Meso|Star> <contact@meso-star.com>.
-Copyright © 2018, 2019, 2021 CNRS
+Copyright © 2018, 2019, 2021 CNRS.
Copyright © 2018, 2019 Université Paul Sabatier
<contact-edstar@laplace.univ-tlse.fr>.
diff --git a/src/atmosphere/htrdr_atmosphere_args.c b/src/atmosphere/htrdr_atmosphere_args.c
@@ -65,7 +65,7 @@ print_help(const char* cmd)
printf(
" -m MIE filename of the Mie's data.\n");
printf(
-" -n SKY-NAME name used to identify the sky in the MATERIALS file.\n"
+" -n SKY-NAME name used to identify the sky in the MATERIALS file.\n"
" Its default value is `%s'.\n",
HTRDR_ATMOSPHERE_ARGS_DEFAULT.sky_mtl_name);
printf(
@@ -76,7 +76,7 @@ print_help(const char* cmd)
" written to standard output.\n");
printf(
" -p <rectangle> switch in flux computation by defining the rectangular\n"
-" sensor onto wich the flux is computed. Refer to the\n"
+" sensor onto which the flux is computed. Refer to the\n"
" %s man page for the list of rectangle options.\n", cmd);
printf(
" -R infinitely repeat the ground along the X and Y axis.\n");
diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c
@@ -310,7 +310,7 @@ draw_pixel_flux
time.nweights += 1;
}
- /* Save the per realisation integration time */
+ /* Write the accumulators */
pixel->flux = flux;
pixel->time = time;
}
diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c
@@ -166,6 +166,7 @@ setup_camera
{
double proj_ratio = 0;
ASSERT(cmd && args && args->image.definition[0] && args->image.definition[1]);
+ ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE);
proj_ratio =
(double)args->image.definition[0]
@@ -182,6 +183,40 @@ setup_camera
}
static res_T
+setup_flux_map
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ ASSERT(cmd && args);
+ ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP);
+ return htrdr_rectangle_create
+ (cmd->htrdr,
+ args->flux_map.size,
+ args->flux_map.position,
+ args->flux_map.target,
+ args->flux_map.up,
+ &cmd->flux_map);
+}
+
+static res_T
+setup_sensor
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ res_T res = RES_OK;
+ switch(cmd->output_type) {
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
+ res = setup_flux_map(cmd, args);
+ break;
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
+ res = setup_camera(cmd, args);
+ break;
+ default: /* Nothing to do */ break;
+ }
+ return res;
+}
+
+static res_T
setup_laser
(struct htrdr_combustion* cmd,
const struct htrdr_combustion_args* args)
@@ -250,7 +285,9 @@ setup_buffer
res_T res = RES_OK;
ASSERT(cmd && args);
- if(cmd->output_type != HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE) goto exit;
+ if(cmd->output_type != HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP
+ && cmd->output_type != HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE)
+ goto exit;
combustion_get_pixel_format(cmd, &pixfmt);
@@ -394,6 +431,9 @@ dump_laser_sheet(const struct htrdr_combustion* cmd)
res_T res = RES_OK;
ASSERT(cmd);
+ htrdr_log(cmd->htrdr, "Write laser sheet to '%s'.\n",
+ str_cget(&cmd->output_name));
+
/* Compute the extent of the geometry that will represent the laser sheet */
extent = compute_laser_mesh_extent(cmd);
@@ -462,6 +502,7 @@ combustion_release(ref_T* ref)
if(cmd->mats) htrdr_materials_ref_put(cmd->mats);
if(cmd->medium) ATRSTM(ref_put(cmd->medium));
if(cmd->camera) htrdr_camera_ref_put(cmd->camera);
+ if(cmd->flux_map) htrdr_rectangle_ref_put(cmd->flux_map);
if(cmd->laser) htrdr_combustion_laser_ref_put(cmd->laser);
if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
@@ -509,7 +550,7 @@ htrdr_combustion_create
if(res != RES_OK) goto error;
res = setup_geometry(cmd, args);
if(res != RES_OK) goto error;
- res = setup_camera(cmd, args);
+ res = setup_sensor(cmd, args);
if(res != RES_OK) goto error;
res = setup_laser(cmd, args);
if(res != RES_OK) goto error;
@@ -554,6 +595,9 @@ htrdr_combustion_run(struct htrdr_combustion* cmd)
ASSERT(cmd);
switch(cmd->output_type) {
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
+ res = combustion_draw_map(cmd);
+ break;
case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
res = combustion_draw_map(cmd);
break;
@@ -585,6 +629,15 @@ combustion_get_pixel_format
{
ASSERT(cmd && fmt);
(void)cmd;
- fmt->size = sizeof(struct combustion_pixel);
- fmt->alignment = ALIGNOF(struct combustion_pixel);
+ switch(cmd->output_type) {
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
+ fmt->size = sizeof(struct combustion_pixel_flux);
+ fmt->alignment = ALIGNOF(struct combustion_pixel_flux);
+ break;
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
+ fmt->size = sizeof(struct combustion_pixel_image);
+ fmt->alignment = ALIGNOF(struct combustion_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
}
diff --git a/src/combustion/htrdr_combustion_args.c b/src/combustion/htrdr_combustion_args.c
@@ -73,6 +73,11 @@ print_help(const char* cmd)
" -l <laser> define the geometry of the laser sheet. Refer to the\n"
" man page for the list of laser options.\n");
printf(
+" -R <rectangle> switch in flux computation bu defining the\n"
+" rectangular sensor onto which the flux is computed.\n"
+" Refer to the man page for the list of rectangle\n"
+" options.\n");
+ printf(
" -m TETRAHEDRA path toward the volumetric mesh.\n");
printf(
" -N precompute the tetrahedra normals.\n");
@@ -239,9 +244,10 @@ htrdr_combustion_args_init
*args = HTRDR_COMBUSTION_ARGS_DEFAULT;
- while((opt = getopt(argc, argv, "C:D:d:F:fg:hi:l:m:NO:o:p:r:sT:t:V:vw:")) != -1) {
+ while((opt = getopt(argc, argv, "C:D:d:F:fg:hi:l:m:NO:o:p:R:r:sT:t:V:vw:")) != -1) {
switch(opt) {
case 'C':
+ args->output_type = HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE;
res = htrdr_args_camera_parse(&args->camera, optarg);
break;
case 'D':
@@ -275,6 +281,10 @@ htrdr_combustion_args_init
case 'o': args->path_output = optarg; break;
case 'p': args->path_therm_props = optarg; break;
case 'r': args->path_refract_ids = optarg; break;
+ case 'R':
+ args->output_type = HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP;
+ res = htrdr_args_rectangle_parse(&args->flux_map, optarg);
+ break;
case 's': args->use_simd = 1; break;
case 'T':
res = cstr_to_double(optarg, &args->optical_thickness);
diff --git a/src/combustion/htrdr_combustion_args.h.in b/src/combustion/htrdr_combustion_args.h.in
@@ -24,6 +24,7 @@
#include <limits.h> /* UINT_MAX support */
enum htrdr_combustion_args_output_type {
+ HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP,
HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE,
HTRDR_COMBUSTION_ARGS_OUTPUT_LASER_SHEET,
HTRDR_COMBUSTION_ARGS_OUTPUT_OCTREES,
@@ -62,6 +63,7 @@ struct htrdr_combustion_args {
const char* path_output; /* Name of the output file */
struct htrdr_args_camera camera; /* Pinhole Camera */
+ struct htrdr_args_rectangle flux_map; /* Flux map */
struct htrdr_args_rectangle laser; /* Laser surface emission */
double wavelength; /* Wavelength of the laser in nanometre */
@@ -98,6 +100,7 @@ struct htrdr_combustion_args {
NULL, /* Output path */ \
\
HTRDR_ARGS_CAMERA_DEFAULT__, /* Pinhole camera */ \
+ HTRDR_ARGS_RECTANGLE_DEFAULT__, /* Flux map */ \
\
HTRDR_ARGS_RECTANGLE_DEFAULT__, /* Laser surface emission */ \
@HTRDR_COMBUSTION_ARGS_DEFAULT_WAVELENGTH@, /* Wavelength in nanometre */ \
diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h
@@ -40,16 +40,27 @@ struct htrdr_rectangle;
struct ssf_phase;
struct ssp_rng;
-struct combustion_pixel {
+struct combustion_pixel_flux {
+ struct htrdr_accum flux; /* In W/m^2 */
+ struct htrdr_accum time; /* In microseconds */
+};
+#define COMBUSTION_PIXEL_FLUX_NULL__ { \
+ HTRDR_ACCUM_NULL__, /* Flux */ \
+ HTRDR_ACCUM_NULL__, /* Time */ \
+}
+static const struct combustion_pixel_flux COMBUSTION_PIXEL_FLUX_NULL =
+ COMBUSTION_PIXEL_FLUX_NULL__;
+
+struct combustion_pixel_image {
struct htrdr_estimate radiance; /* In W/m^2/sr */
struct htrdr_accum time; /* In microseconds */
};
-#define COMBUSTION_PIXEL_NULL__ { \
+#define COMBUSTION_PIXEL_IMAGE_NULL__ { \
HTRDR_ESTIMATE_NULL__, /* Radiance */ \
HTRDR_ACCUM_NULL__, /* Time */ \
}
-static const struct combustion_pixel COMBUSTION_PIXEL_NULL =
- COMBUSTION_PIXEL_NULL__;
+static const struct combustion_pixel_image COMBUSTION_PIXEL_IMAGE_NULL =
+ COMBUSTION_PIXEL_IMAGE_NULL__;
struct htrdr_combustion {
struct htrdr_geometry* geom; /* Combustion chamber geometry */
@@ -57,6 +68,7 @@ struct htrdr_combustion {
struct atrstm* medium; /* Semi transparent medium */
struct htrdr_camera* camera; /* Pinhole camera */
+ struct htrdr_rectangle* flux_map; /* Flux map */
struct htrdr_combustion_laser* laser; /* Laser sheet */
double wavelength; /* Wavelength of the laser in nanometer */
diff --git a/src/combustion/htrdr_combustion_draw_map.c b/src/combustion/htrdr_combustion_draw_map.c
@@ -21,6 +21,7 @@
#include "core/htrdr_camera.h"
#include "core/htrdr_draw_map.h"
#include "core/htrdr_log.h"
+#include "core/htrdr_rectangle.h"
#include <star/ssp.h>
@@ -30,7 +31,7 @@
* Helper functions
******************************************************************************/
static void
-draw_pixel
+draw_pixel_image
(struct htrdr* htrdr,
const struct htrdr_draw_pixel_args* args,
void* data)
@@ -38,7 +39,7 @@ draw_pixel
struct htrdr_accum radiance = HTRDR_ACCUM_NULL;
struct htrdr_accum time = HTRDR_ACCUM_NULL;
struct htrdr_combustion* cmd = NULL;
- struct combustion_pixel* pixel = NULL;
+ struct combustion_pixel_image* pixel = NULL;
size_t isamp;
res_T res = RES_OK;
ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
@@ -96,31 +97,131 @@ draw_pixel
}
static void
-dump_pixel
- (const struct combustion_pixel* pix,
- struct htrdr_accum* time_acc, /* May be NULL */
+draw_pixel_flux
+ (struct htrdr* htrdr,
+ const struct htrdr_draw_pixel_args* args,
+ void* data)
+{
+ struct htrdr_accum flux = HTRDR_ACCUM_NULL;
+ struct htrdr_accum time = HTRDR_ACCUM_NULL;
+ struct htrdr_combustion* cmd = NULL;
+ struct combustion_pixel_flux* pixel = NULL;
+ size_t isamp;
+ res_T res = RES_OK;
+ ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
+ (void)htrdr;
+
+ cmd = args->context;
+ pixel = data;
+
+ FOR_EACH(isamp, 0, args->spp) {
+ struct time t0, t1;
+ double pix_samp[2];
+ double ray_org[3];
+ double ray_dir[3];
+ double normal[3];
+ double weight;
+ double usec;
+
+ /* Begin the registration of the time spent in the realisation */
+ time_current(&t0);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = (double)args->pixel_coord[0] + ssp_rng_canonical(args->rng);
+ pix_samp[1] = (double)args->pixel_coord[1] + ssp_rng_canonical(args->rng);
+ pix_samp[0] *= args->pixel_normalized_size[0];
+ pix_samp[1] *= args->pixel_normalized_size[1];
+
+ /* Retrieve the world space position of pix_samp */
+ htrdr_rectangle_sample_pos(cmd->flux_map, pix_samp, ray_org);
+
+ /* Sample a ray direction */
+ htrdr_rectangle_get_normal(cmd->flux_map, normal);
+ ssp_ran_hemisphere_cos(args->rng, normal, ray_dir, NULL);
+
+ /* Backward trace the path */
+ res = combustion_compute_radiance_sw(cmd, args->ithread,
+ args->rng, ray_org, ray_dir, &weight);
+ if(res != RES_OK) continue; /* Reject the path */
+ weight *= PI; /* Transform form W/m^2/sr to W/m^2 */
+
+ /* End the registration of the per realisation time */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update the pixel accumulator */
+ flux.sum_weights += weight;
+ flux.sum_weights_sqr += weight*weight;
+ flux.nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Write the accumulators */
+ pixel->flux = flux;
+ pixel->time = time;
+}
+
+static INLINE void
+dump_accum
+ (const struct htrdr_accum* acc, /* Accum to dump */
+ struct htrdr_accum* out_acc, /* May be NULL */
FILE* stream)
{
- struct htrdr_estimate time = HTRDR_ESTIMATE_NULL;
- ASSERT(pix && stream && pix->time.nweights);
+ ASSERT(acc && stream);
- htrdr_accum_get_estimation(&pix->time, &time);
+ if(acc->nweights == 0) {
+ fprintf(stream, "0 0 ");
+ } else {
+ struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL;
- fprintf(stream, "%g %g 0 0 0 0 %g %g\n",
- pix->radiance.E, pix->radiance.SE, time.E,time.SE);
+ htrdr_accum_get_estimation(acc, &estimate);
+ fprintf(stream, "%g %g ", estimate.E, estimate.SE);
- if(time_acc) {
- time_acc->sum_weights += pix->time.sum_weights;
- time_acc->sum_weights_sqr += pix->time.sum_weights_sqr;
- time_acc->nweights += pix->time.nweights;
+ if(out_acc) {
+ out_acc->sum_weights += acc->sum_weights;
+ out_acc->sum_weights_sqr += acc->sum_weights_sqr;
+ out_acc->nweights += acc->nweights;
+ }
}
}
+static void
+dump_pixel_flux
+ (const struct combustion_pixel_flux* pix,
+ struct htrdr_accum* time_acc, /* May be NULL */
+ struct htrdr_accum* flux_acc, /* May be NULL */
+ FILE* stream)
+{
+ ASSERT(pix && stream);
+ dump_accum(&pix->flux, flux_acc, stream);
+ fprintf(stream, "0 0 0 0 ");
+ dump_accum(&pix->time, time_acc, stream);
+ fprintf(stream, "\n");
+}
+
+static void
+dump_pixel_image
+ (const struct combustion_pixel_image* pix,
+ struct htrdr_accum* time_acc, /* May be NULL */
+ FILE* stream)
+{
+ ASSERT(pix && stream);
+ fprintf(stream, "%g %g ", pix->radiance.E, pix->radiance.SE);
+ fprintf(stream, "0 0 0 0 ");
+ dump_accum(&pix->time, time_acc, stream);
+ fprintf(stream, "\n");
+}
+
static res_T
dump_buffer
(struct htrdr_combustion* cmd,
struct htrdr_buffer* buf,
struct htrdr_accum* time_acc, /* May be NULL */
+ struct htrdr_accum* flux_acc, /* May be NULL */
FILE* stream)
{
struct htrdr_pixel_format pixfmt;
@@ -138,9 +239,17 @@ dump_buffer
FOR_EACH(y, 0, layout.height) {
FOR_EACH(x, 0, layout.width) {
- const struct combustion_pixel* pix = htrdr_buffer_at(buf, x, y);
- ASSERT(IS_ALIGNED(pix, pixfmt.alignment));
- dump_pixel(pix, time_acc, stream);
+ void* pix_raw = htrdr_buffer_at(buf, x, y);
+ ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment));
+ if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP) {
+ const struct combustion_pixel_flux* pix = pix_raw;
+ dump_pixel_flux(pix, time_acc, flux_acc, stream);
+ } else if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE) {
+ const struct combustion_pixel_image* pix = pix_raw;
+ dump_pixel_image(pix, time_acc, stream);
+ } else {
+ FATAL("Unreachable code.\n");
+ }
}
fprintf(stream, "\n");
}
@@ -155,12 +264,22 @@ combustion_draw_map(struct htrdr_combustion* cmd)
{
struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL;
struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
+ struct htrdr_accum flux_acc = HTRDR_ACCUM_NULL;
struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL;
+ struct htrdr_estimate flux = HTRDR_ESTIMATE_NULL;
res_T res = RES_OK;
ASSERT(cmd);
/* Setup the draw map input arguments */
- args.draw_pixel = draw_pixel;
+ switch(cmd->output_type) {
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
+ args.draw_pixel = draw_pixel_image;
+ break;
+ case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
+ args.draw_pixel = draw_pixel_flux;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
args.buffer_layout = cmd->buf_layout;
args.spp = cmd->spp;
args.context = cmd;
@@ -172,7 +291,7 @@ combustion_draw_map(struct htrdr_combustion* cmd)
if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
/* Write buffer to output */
- res = dump_buffer(cmd, cmd->buf, &path_time_acc, cmd->output);
+ res = dump_buffer(cmd, cmd->buf, &path_time_acc, &flux_acc, cmd->output);
if(res != RES_OK) goto error;
htrdr_accum_get_estimation(&path_time_acc, &path_time);
@@ -180,6 +299,13 @@ combustion_draw_map(struct htrdr_combustion* cmd)
"Time per radiative path (in micro seconds): %g +/- %g\n",
path_time.E, path_time.SE);
+ if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP) {
+ htrdr_accum_get_estimation(&flux_acc, &flux);
+ htrdr_log(cmd->htrdr,
+ "Radiative flux density (in W/m^2): %g +/- %g\n",
+ flux.E, flux.SE);
+ }
+
exit:
return res;
error:
