commit a96a5c89ca005fef90b4fb987a4d790abaea113d
parent cd9255ceecefacde912d4ba328b693dd7214dbb8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 26 Mar 2021 09:57:59 +0100
Add Henyey-Greenstein as an workaround of the RDG-FA phase function
Diffstat:
3 files changed, 65 insertions(+), 35 deletions(-)
diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c
@@ -38,22 +38,21 @@
* Helper functions
******************************************************************************/
static void
-release_rdgfa(struct htrdr_combustion* cmd)
+release_phase_functions
+ (struct htrdr_combustion* cmd,
+ struct ssf_phase* phases[])
{
size_t i;
ASSERT(cmd);
- if(!cmd->rdgfa_phase_functions) return; /* Nothing to release */
-
+ if(!phases) return; /* Nothing to release */
FOR_EACH(i, 0, htrdr_get_threads_count(cmd->htrdr)) {
- if(cmd->rdgfa_phase_functions[i]) {
- SSF(phase_ref_put(cmd->rdgfa_phase_functions[i]));
+ if(phases[i]) {
+ SSF(phase_ref_put(phases[i]));
}
}
-
- MEM_RM(htrdr_get_allocator(cmd->htrdr), cmd->rdgfa_phase_functions);
- cmd->rdgfa_phase_functions = NULL;
+ MEM_RM(htrdr_get_allocator(cmd->htrdr), phases);
}
static res_T
@@ -166,46 +165,47 @@ setup_laser
}
static res_T
-setup_rdgfa
+setup_phase_functions
(struct htrdr_combustion* cmd,
- const struct htrdr_combustion_args* args)
+ const struct ssf_phase_type* phase_type,
+ struct ssf_phase** out_phases[])
{
struct mem_allocator* allocator = NULL;
+ struct ssf_phase** phases = NULL;
size_t nthreads;
size_t i;
res_T res = RES_OK;
- ASSERT(cmd);
- (void)args; /* Not use */
+ ASSERT(cmd && phase_type && out_phases);
nthreads = htrdr_get_threads_count(cmd->htrdr);
allocator = htrdr_get_allocator(cmd->htrdr);
- /* Allocate the list of per thread RDG-FA */
- cmd->rdgfa_phase_functions = MEM_CALLOC
- (allocator, nthreads, sizeof(*cmd->rdgfa_phase_functions));
- if(!cmd->rdgfa_phase_functions) {
+ /* Allocate the list of per thread phase function */
+ phases = MEM_CALLOC(allocator, nthreads, sizeof(*phases));
+ if(!phases) {
htrdr_log_err(cmd->htrdr,
"Could not allocate the per thread RDG-FA phase function.\n");
res = RES_MEM_ERR;
goto error;
}
- /* Create the per thread RDG-FA */
+ /* Create the per thread phase function */
FOR_EACH(i, 0, nthreads) {
- res = ssf_phase_create
- (allocator, &ssf_phase_rdgfa, cmd->rdgfa_phase_functions + i);
+ res = ssf_phase_create(allocator, phase_type, phases+i);
if(res != RES_OK) {
htrdr_log_err(cmd->htrdr,
- "Could not create the RDG-FA phase function for the thread %lu -- %s.\n",
+ "Could not create the phase function for the thread %lu -- %s.\n",
(unsigned long)i, res_to_cstr(res));
goto error;
}
}
exit:
+ *out_phases = phases;
return res;
error:
- release_rdgfa(cmd);
+ release_phase_functions(cmd, phases);
+ phases = NULL;
goto exit;
}
@@ -332,7 +332,8 @@ combustion_release(ref_T* ref)
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);
- release_rdgfa(cmd);
+ release_phase_functions(cmd, cmd->rdgfa_phase_functions);
+ release_phase_functions(cmd, cmd->hg_phase_functions);
str_release(&cmd->output_name);
htrdr = cmd->htrdr;
@@ -378,7 +379,9 @@ htrdr_combustion_create
if(res != RES_OK) goto error;
res = setup_laser(cmd, args);
if(res != RES_OK) goto error;
- res = setup_rdgfa(cmd, args);
+ res = setup_phase_functions(cmd, &ssf_phase_rdgfa, &cmd->rdgfa_phase_functions);
+ if(res != RES_OK) goto error;
+ res = setup_phase_functions(cmd, &ssf_phase_hg, &cmd->hg_phase_functions);
if(res != RES_OK) goto error;
res = setup_buffer(cmd, args);
if(res != RES_OK) goto error;
diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h
@@ -57,6 +57,7 @@ struct htrdr_combustion {
double wavelength; /* Wavelength of the laser in nanometer */
struct ssf_phase** rdgfa_phase_functions; /* Per thread RDG-FA phase func */
+ struct ssf_phase** hg_phase_functions; /* Per thread Henyey-Greenstein func */
struct htrdr_buffer_layout buf_layout;
struct htrdr_buffer* buf; /* NULL on non master processes */
diff --git a/src/combustion/htrdr_combustion_compute_radiance_sw.c b/src/combustion/htrdr_combustion_compute_radiance_sw.c
@@ -30,6 +30,8 @@
#include <rsys/double4.h>
#include <rsys/dynamic_array.h>
+#define USE_HG_PHASE_FUNCTION 1
+
/* Define a position along the ray into the semi-transparent medium */
struct position {
double distance; /* Distance up to the position */
@@ -178,8 +180,12 @@ sample_position_hit_filter
/* Fetch the radiative coefficient at the collision position */
ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords));
- ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
- k = radcoefs[ctx->radcoef];
+ if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) {
+ k = 0;
+ } else {
+ ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
+ k = radcoefs[ctx->radcoef];
+ }
r = ssp_rng_canonical(ctx->rng);
@@ -375,8 +381,12 @@ sample_scattering_limited_hit_filter
/* Fetch the radiative coefficient at the collision position */
ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords));
- ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
- ks = radcoefs[ATRSTM_RADCOEF_Ks];
+ if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) {
+ ks = 0;
+ } else {
+ ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
+ ks = radcoefs[ATRSTM_RADCOEF_Ks];
+ }
r = ssp_rng_canonical(ctx->rng);
@@ -479,13 +489,15 @@ laser_once_scattered
const double pos[3],
const double dir[3])
{
+#ifndef USE_HG_PHASE_FUNCTION
/* RDG-FA phase function */
struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL;
struct atrstm_fetch_rdgfa_args fetch_rdgfa_args =
ATRSTM_FETCH_RDGFA_ARGS_DEFAULT;
struct ssf_phase_rdgfa_setup_args setup_rdgfa_args =
SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT;
- struct ssf_phase* rdgfa = NULL;
+#endif
+ struct ssf_phase* phase = NULL;
/* Scattering position into the laser sheet */
struct position_limited sc_sample = POSITION_LIMITED_NULL;
@@ -558,6 +570,11 @@ laser_once_scattered
Tr_ext_xsc_lse = transmissivity(cmd, rng, ATRSTM_RADCOEF_Kext, xsc, wi, range);
if(Tr_ext_xsc_lse == 0) return 0; /* No laser contribution */
+#ifdef USE_HG_PHASE_FUNCTION
+ phase = cmd->hg_phase_functions[ithread];
+ SSF(phase_hg_setup(phase, 0));
+ (void)wlen;
+#else
/* Retrieve the RDG-FA phase function parameters from the semi transparent
* medium */
fetch_rdgfa_args.wavelength = wlen;
@@ -569,15 +586,16 @@ laser_once_scattered
ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param));
/* Setup the RDG-FA phase function */
- rdgfa = cmd->rdgfa_phase_functions[ithread];
+ phase = cmd->rdgfa_phase_functions[ithread];
setup_rdgfa_args.wavelength = rdgfa_param.wavelength;
setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension;
setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius;
- SSF(phase_rdgfa_setup(rdgfa, &setup_rdgfa_args));
+ SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args));
+#endif
/* Evaluate the phase function at the scattering position */
d3_minus(wo, dir); /* Ensure SSF convention */
- R = ssf_phase_eval(rdgfa, wo, wi); /* In sr^-1 */
+ R = ssf_phase_eval(phase, wo, wi); /* In sr^-1 */
laser_flux_density = htrdr_combustion_laser_get_flux_density(cmd->laser);
@@ -600,13 +618,15 @@ combustion_compute_radiance_sw
const double pos_in[3],
const double dir_in[3])
{
+#ifndef USE_HG_PHASE_FUNCTION
/* RDG-FA phase function */
struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL;
struct atrstm_fetch_rdgfa_args fetch_rdgfa_args =
ATRSTM_FETCH_RDGFA_ARGS_DEFAULT;
struct ssf_phase_rdgfa_setup_args setup_rdgfa_args =
SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT;
- struct ssf_phase* rdgfa = NULL;
+#endif
+ struct ssf_phase* phase = NULL;
/* Transmissivity between the probe position (i.e. 'pos_in') and the current
* scattering position over the reverse scattering path */
@@ -626,7 +646,6 @@ combustion_compute_radiance_sw
d3_set(dir, dir_in);
d2(range, 0, FLT_MAX);
- rdgfa = cmd->rdgfa_phase_functions[ithread];
wlen = htrdr_combustion_laser_get_wavelength(cmd->laser);
Tr_abs = 1;
@@ -661,6 +680,11 @@ combustion_compute_radiance_sw
pos[1] = pos[1] + dir[1] * scattering.distance;
pos[2] = pos[2] + dir[2] * scattering.distance;
+#ifdef USE_HG_PHASE_FUNCTION
+ phase = cmd->hg_phase_functions[ithread];
+ SSF(phase_hg_setup(phase, 0));
+ (void)wlen;
+#else
/* Retrieve the RDG-FA phase function parameters */
fetch_rdgfa_args.wavelength = wlen;
fetch_rdgfa_args.prim = scattering.prim;
@@ -671,14 +695,16 @@ combustion_compute_radiance_sw
ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param));
/* Setup the RDG-FA phase function corresponding to the scattering event */
+ phase = cmd->rdgfa_phase_functions[ithread];
setup_rdgfa_args.wavelength = rdgfa_param.wavelength;
setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension;
setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius;
- SSF(phase_rdgfa_setup(rdgfa, &setup_rdgfa_args));
+ SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args));
+#endif
/* Sample a new optical path direction from the phase function */
d3_minus(wo, dir); /* Ensure SSF convention */
- ssf_phase_sample(rdgfa, rng, wo, wi, NULL);
+ ssf_phase_sample(phase, rng, wo, wi, NULL);
/* Update the optical path direction */
dir[0] = wi[0];
