commit 942ee2563dab26ab1d95762ef50ea6cdeb2c435f
parent aa183e6471ab84e12f5239eb358692fb37af6047
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Oct 2025 15:55:30 +0200
Add an internal structure for external source properties.
External source properties depend on time (position and power). But
since radiative transfer is instantaneous, it is sufficient to retrieve
them once per radiative path sampling. This saves computation time, as
calculating position and power can be very time-consuming.
The source_sample function now takes these source properties as
parameters in order to calculate a source sample, i.e., a direction
toward the source. So, by construction, a sample can be generated only
if the source properties have been retrieved.
The Monte Carlo algorithm that manages net flow is updated to support
these changes.
Diffstat:
3 files changed, 89 insertions(+), 80 deletions(-)
diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h
@@ -115,6 +115,7 @@ static res_T
XD(compute_incident_diffuse_flux)
(struct sdis_scene* scn,
struct ssp_rng* rng,
+ const struct source_props* props,
const double in_pos[DIM], /* position */
const double in_N[DIM], /* Surface normal. (Away from the surface) */
const double time,
@@ -128,7 +129,7 @@ XD(compute_incident_diffuse_flux)
double N[3] = {0}; /* Surface normal. Always 3D */
size_t nbounces = 0; /* For debug */
res_T res = RES_OK;
- ASSERT(in_pos && in_N && in_hit && diffuse_flux);
+ ASSERT(props && in_pos && in_N && in_hit && diffuse_flux);
/* Local copy of input argument */
dX(set)(pos, in_pos);
@@ -142,7 +143,7 @@ XD(compute_incident_diffuse_flux)
for(;;) {
/* External sources */
- struct source_sample src_sample = SOURCE_SAMPLE_NULL;
+ struct source_sample samp = SOURCE_SAMPLE_NULL;
/* Interface */
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
@@ -201,12 +202,11 @@ XD(compute_incident_diffuse_flux)
/* Calculate the direct contribution if the rebound is specular */
if(bounce.cpnt == BRDF_SPECULAR) {
- res = source_trace_to(scn->source, pos, bounce.dir, time, &src_sample);
+ res = source_trace_to(scn->source, props, pos, bounce.dir, &samp);
if(res != RES_OK) goto error;
- if(!SOURCE_SAMPLE_NONE(&src_sample)) {
- const double Ld = XD(direct_contribution)
- (scn, &src_sample, pos, enc_id, &hit);
+ if(!SOURCE_SAMPLE_NONE(&samp)) {
+ double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit);
L = Ld; /* [W/m^2/sr] */
}
@@ -217,9 +217,9 @@ XD(compute_incident_diffuse_flux)
/* Sample an external source to handle its direct contribution at the
* bounce position */
- res = source_sample(scn->source, rng, pos, time, &src_sample);
+ res = source_sample(scn->source, props, rng, pos, &samp);
CHK(res == RES_OK);
- cos_theta = d3_dot(src_sample.dir, N);
+ cos_theta = d3_dot(samp.dir, N);
/* The source is behind the surface */
if(cos_theta <= 0) {
@@ -227,9 +227,8 @@ XD(compute_incident_diffuse_flux)
/* The source is above the surface */
} else {
- const double Ld = XD(direct_contribution)
- (scn, &src_sample, pos, enc_id, &hit);
- L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2/sr] */
+ double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit);
+ L = Ld * cos_theta / (PI * samp.pdf); /* [W/m^2/sr] */
}
}
diffuse_flux->reflected += L; /* [W/m^2/sr] */
@@ -257,7 +256,8 @@ XD(handle_external_net_flux)
struct sdis_green_external_flux_terms green =
SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL;
- /* Sampling external sources */
+ /* External source */
+ struct source_props src_props = SOURCE_PROPS_NULL;
struct source_sample src_sample = SOURCE_SAMPLE_NULL;
/* External flux */
@@ -312,9 +312,11 @@ XD(handle_external_net_flux)
if(emissivity == 0) goto exit;
- /* Sample the external source */
- res = source_sample
- (scn->source, rng, frag.P, frag.time, &src_sample);
+ res = source_get_props(scn->source, frag.time, &src_props);
+ if(res != RES_OK) goto error;
+
+ /* Sample a direction toward the source to add its direct contribution */
+ res = source_sample(scn->source, &src_props, rng, frag.P, &src_sample);
if(res != RES_OK) goto error;
/* Setup the normal to ensure that it points toward the fluid medium */
@@ -331,8 +333,8 @@ XD(handle_external_net_flux)
}
/* Calculate the incident diffuse flux [W/m^2] */
- res = XD(compute_incident_diffuse_flux)(scn, rng, frag.P, N, frag.time,
- enc_ids[frag.side], args->XD(hit), &incident_flux_diffuse);
+ res = XD(compute_incident_diffuse_flux)(scn, rng, &src_props, frag.P, N,
+ frag.time, enc_ids[frag.side], args->XD(hit), &incident_flux_diffuse);
if(res != RES_OK) goto error;
/* Calculate the incident flux without the part scattered by the environment.
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -163,43 +163,25 @@ sdis_source_get_id(const struct sdis_source* source)
res_T
source_sample
(const struct sdis_source* src,
+ const struct source_props* props,
struct ssp_rng* rng,
const double pos[3],
- const double time,
struct source_sample* sample)
{
- double src_pos[3]; /* [m] */
double main_dir[3];
double half_angle; /* [radians] */
double cos_half_angle; /* [radians] */
double dst; /* [m] */
- double radius; /* Source radius [m] */
- double power; /* Source power [W] */
- double area; /* Source area [m^2] */
res_T res = RES_OK;
ASSERT(src && rng && pos && sample);
- /* Retrieve current source position, radius and power */
- src->spherical.position(time, src_pos, src->data);
- power = src->spherical.power(time, src->data);
- radius = src->spherical.radius;
-
- if(power < 0) {
- log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
- FUNC_NAME, power);
- res = RES_BAD_ARG;
- goto error;
- }
-
- area = 4*PI*radius*radius; /* [m^2] */
-
/* compute the direction of `pos' toward the center of the source */
- d3_sub(main_dir, src_pos, pos);
+ d3_sub(main_dir, props->pos, pos);
/* Normalize the direction and keep the distance from `pos' to the center of
* the source */
dst = d3_normalize(main_dir, main_dir);
- if(dst <= radius) {
+ if(dst <= props->radius) {
log_err(src->dev,
"%s: the position from which the external source is sampled "
"is included in the source:\n"
@@ -208,34 +190,32 @@ source_sample
"\tposition = %g, %g, %g\n"
"\ttime = %g\n"
"\tdistance from position to source = %g\n",
- FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst);
res = RES_BAD_ARG;
goto error;
}
/* Point source */
- if(area == 0) {
+ if(props->area == 0) {
d3_set(sample->dir, main_dir);
sample->pdf = 1;
sample->dst = dst;
- sample->power = power; /* [W] */
sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */
- sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+ sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
/* Spherical source */
} else {
/* Sample the source according to its solid angle,
* i.e. 2*PI*(1 - cos(half_angle)) */
- half_angle = asin(radius/dst);
+ half_angle = asin(props->radius/dst);
cos_half_angle = cos(half_angle);
ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */
(rng, main_dir, cos_half_angle, sample->dir, &sample->pdf);
/* Set other sample variables */
- sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->power = power; /* [W] */
- sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
- sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+ sample->dst = dst - props->radius; /* From pos to source boundaries [m] */
+ sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */
+ sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
}
exit:
@@ -247,47 +227,31 @@ error:
res_T
source_trace_to
(const struct sdis_source* src,
+ const struct source_props* props,
const double pos[3], /* Ray origin */
const double dir[3], /* Ray direction */
- const double time, /* Time at which ray is traced */
struct source_sample* sample)
{
- double src_pos[3]; /* [m] */
double main_dir[3];
- double radius; /* [m] */
- double power; /* [W] */
double dst; /* Distance from pos to the source center [m] */
double half_angle; /* [radian] */
res_T res = RES_OK;
- ASSERT(src && pos && dir && sample);
+ ASSERT(src && props && pos && dir && sample);
ASSERT(d3_is_normalized(dir));
- radius = src->spherical.radius;
-
/* Point sources cannot be targeted */
- if(radius == 0) {
+ if(props->radius == 0) {
*sample = SOURCE_SAMPLE_NULL;
goto exit;
}
- /* Retrieve current source position and power */
- src->spherical.position(time, src_pos, src->data);
- power = src->spherical.power(time, src->data);
-
- if(power < 0) {
- log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
- FUNC_NAME, power);
- res = RES_BAD_ARG;
- goto error;
- }
-
/* compute the direction of `pos' toward the center of the source */
- d3_sub(main_dir, src_pos, pos);
+ d3_sub(main_dir, props->pos, pos);
/* Normalize the direction and keep the distance from `pos' to the center of
* the source */
dst = d3_normalize(main_dir, main_dir);
- if(dst <= radius) {
+ if(dst <= props->radius) {
log_err(src->dev,
"%s: the position from which the external source is targeted "
"is included in the source:\n"
@@ -296,13 +260,13 @@ source_trace_to
"\tposition = %g, %g, %g\n"
"\ttime = %g\n"
"\tdistance from position to source = %g\n",
- FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst);
res = RES_BAD_ARG;
goto error;
}
/* Compute the half angle of the source as seen from pos */
- half_angle = asin(radius/dst);
+ half_angle = asin(props->radius/dst);
/* The source is missed */
if(d3_dot(dir, main_dir) < cos(half_angle)) {
@@ -310,14 +274,11 @@ source_trace_to
/* The source is intersected */
} else {
- const double area = 4*PI*radius*radius; /* [m^2] */
-
d3_set(sample->dir, dir);
sample->pdf = 1;
- sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->power = power; /* [W] */
- sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
- sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+ sample->dst = dst - props->radius; /* From pos to source boundaries [m] */
+ sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */
+ sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
}
exit:
@@ -327,6 +288,36 @@ error:
goto exit;
}
+res_T
+source_get_props
+ (const struct sdis_source* src,
+ const double time, /* [s] */
+ struct source_props* props)
+{
+ res_T res = RES_OK;
+ ASSERT(src && props);
+
+ /* Retrieve the source properties */
+ src->spherical.position(time, props->pos, src->data);
+ props->power = src->spherical.power(time, src->data);
+ props->radius = src->spherical.radius;
+
+ if(props->power < 0) {
+ log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ FUNC_NAME, props->power);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ props->area = 4*PI*props->radius*props->radius; /* [m^2] */
+ props->time = time; /* [s] */
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
double /* [W] */
source_get_power(const struct sdis_source* src, const double time /* [s] */)
{
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -22,11 +22,21 @@
struct sdis_source;
struct ssp_rng;
+struct source_props {
+ double pos[3]; /* [m/fp_to_meter] */
+ double radius; /* [m/fp_to_meter] */
+ double power; /* [W] */
+ double area; /* [m^2/fp_to_meter] */
+ double time; /* [s] */
+};
+#define SOURCE_PROPS_NULL__ {0}
+static const struct source_props SOURCE_PROPS_NULL = SOURCE_PROPS_NULL__;
+
struct source_sample {
double dir[3]; /* Direction _to_ the source */
double pdf; /* pdf of sampled direction */
double dst; /* Distance to the source [m] */
- double power; /* [W] */
+
double radiance; /* [W/m^2/sr] */
/* Radiance relative to power, i.e. the source power is assumed to be equal to
@@ -36,18 +46,24 @@ struct source_sample {
* green function */
double radiance_term; /* [W/m^2/sr] */
};
-#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0, 0, 0}
+#define SOURCE_SAMPLE_NULL__ {0}
static const struct source_sample SOURCE_SAMPLE_NULL = SOURCE_SAMPLE_NULL__;
/* Helper macro used to define whether a sample is valid or not */
#define SOURCE_SAMPLE_NONE(Sample) ((Sample)->pdf == 0)
extern LOCAL_SYM res_T
+source_get_props
+ (const struct sdis_source* source,
+ const double time, /* Time at which props are retrieved [s] */
+ struct source_props* props);
+
+extern LOCAL_SYM res_T
source_sample
(const struct sdis_source* source,
+ const struct source_props* props,
struct ssp_rng* rng,
const double pos[3], /* Position from which the source is sampled */
- const double time, /* Time at which the source is sampled */
struct source_sample* sample);
/* Trace a ray toward the source. The returned sample has a pdf of 1 or 0
@@ -56,9 +72,9 @@ source_sample
extern LOCAL_SYM res_T
source_trace_to
(const struct sdis_source* source,
+ const struct source_props* props,
const double pos[3], /* Ray origin */
const double dir[3], /* Ray direction */
- const double time, /* Time at which ray is traced */
struct source_sample* sample); /* pdf == 0 if no source is reached */
extern LOCAL_SYM double /* [W] */
