commit 0af29fc052ff4d520686da8e7dd91dae4bddd07a
parent d2defc3237bf54c2b3e9d732d2102e18154a17e2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 20 Dec 2023 19:17:39 +0100
Finalize the management of the external flux
Add the calculation of its diffuse part. Even if the sources compile,
this implementation remains a draft that should be highly buggy since it
hasn't been tested. It's still just a milestone to save the current
state. But it's all there, from data management and code structure to
implementation of the Monte Carlo algorithm.
Diffstat:
2 files changed, 330 insertions(+), 33 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
@@ -23,6 +23,82 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Non generic data types and function
+ ******************************************************************************/
+#ifndef SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H
+#define SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H
+
+enum brdf_component {
+ BRDF_SPECULAR,
+ BRDF_DIFFUSE,
+ BRDF_NONE
+};
+
+struct brdf_sample {
+ double dir[3];
+ double pdf;
+ enum brdf_component cpnt;
+};
+#define BRDF_SAMPLE_NULL__ {{0}, 0, BRDF_NONE}
+static const struct brdf_sample BRDF_SAMPLE_NULL = BRDF_SAMPLE_NULL__;
+
+struct brdf {
+ double emissivity;
+ double specular_fraction;
+};
+#define BRDF_NULL__ {0, 0}
+static const struct brdf BRDF_NULL = BRDF_NULL__;
+
+/* Reflect the V wrt the normal N. By convention V points outward the surface.
+ * In fact, this function is a double-precision version of the reflect_3d
+ * function. TODO Clean this "repeat" */
+static FINLINE double*
+reflect(double res[3], const double V[3], const double N[3])
+{
+ double tmp[3];
+ double cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
+ cos_V_N = d3_dot(V, N);
+ d3_muld(tmp, N, 2*cos_V_N);
+ d3_sub(res, tmp, V);
+ return res;
+}
+
+static void
+sample_brdf
+ (const struct brdf* brdf,
+ struct ssp_rng* rng,
+ const double wi[3], /* Incident direction. Point away from the surface */
+ const double N[3], /* Surface normal */
+ struct brdf_sample* sample)
+{
+ double r = 0; /* Random number */
+
+ /* Preconditions */
+ ASSERT(brdf && rng && wi && N && sample);
+ ASSERT(d3_is_normalized(wi) && d3_is_normalized(N));
+ ASSERT(d3_dot(wi, N) > 0);
+
+ r = ssp_rng_canonical(rng);
+
+ /* Sample the specular part */
+ if(r < brdf->specular_fraction) {
+ reflect(sample->dir, wi, N);
+ sample->pdf = 1;
+ sample->cpnt = BRDF_SPECULAR;
+
+ /* Sample the diffuse part */
+ } else {
+ ssp_ran_hemisphere_cos(rng, N, sample->dir, NULL);
+ sample->pdf = 1.0/PI;
+ sample->cpnt = BRDF_DIFFUSE;
+ }
+}
+
+#endif /* SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H */
+
+/*******************************************************************************
* Generic helper functions
******************************************************************************/
static INLINE res_T
@@ -38,6 +114,8 @@ XD(check_handle_external_net_flux_args)
ASSERT(dev && func_name && args);
ASSERT(args->interf && args->frag);
ASSERT(!SXD_HIT_NONE(args->hit));
+ ASSERT(args->h_cond >= 0 && args->h_cond && args->h_radi >= 0);
+ ASSERT(args->h_cond + args->h_cond + args->h_radi > 0);
functor = interface_side_get_external_sources_sampling_functor
(args->interf, args->frag);
@@ -54,6 +132,213 @@ XD(check_handle_external_net_flux_args)
return RES_OK;
}
+static INLINE void
+XD(trace_ray)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ const double dir[3],
+ const double distance,
+ const struct sXd(hit)* hit_from,
+ struct sXd(hit)* hit)
+{
+ struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
+ float ray_org[DIM] = {0};
+ float ray_dir[3] = {0};
+ float ray_range[2] = {0};
+ ASSERT(scn && pos && dir && distance >= 0 && hit_from && hit);
+
+ fX_set_dX(ray_org, pos);
+ f3_set_d3(ray_dir, dir);
+ ray_range[0] = 0;
+ ray_range[1] = (float)distance;
+ filter_data.XD(hit) = *hit_from;
+#if DIM == 2
+ SXD(scene_view_trace_ray_3d
+ (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
+#else
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
+#endif
+}
+
+static INLINE double /* [W/m^2/sr] */
+XD(direct_contribution)
+ (const struct sdis_scene* scn,
+ struct sdis_external_sources_sample* sample,
+ const double pos[DIM],
+ const struct sXd(hit)* hit_from)
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ ASSERT(scn && sample && pos && hit_from);
+
+ /* Is the source hidden */
+ XD(trace_ray)(scn, pos, sample->dir, sample->distance, hit_from, &hit);
+ if(SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */
+
+ return sample->radiance; /* [W/m^2/sr] */
+}
+
+static INLINE void
+XD(setup_fragment)
+ (struct sdis_interface_fragment* frag,
+ const double pos[DIM],
+ const double dir[DIM], /* Direction _toward_ the hit position */
+ const double time, /* Current time */
+ const double N[DIM],/* Surface normal */
+ const struct sXd(hit)* hit)
+{
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ enum sdis_side side = SDIS_SIDE_NULL__;
+ ASSERT(frag && pos && dir && N);
+ ASSERT(dX(is_normalized)(N));
+
+ /* Setup the interface fragment at the intersection position */
+ dX(set)(vtx.P, pos);
+ vtx.time = time;
+ side = dX(dot)(dir, N) < 0 ? SDIS_FRONT : SDIS_BACK;
+ XD(setup_interface_fragment)(frag, &vtx, hit, side);
+}
+
+static INLINE res_T
+XD(setup_brdf)
+ (struct sdis_device* dev,
+ struct brdf* brdf,
+ const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ double epsilon = 0;
+ double alpha = 0;
+ res_T res = RES_OK;
+ ASSERT(brdf && frag);
+
+ epsilon = interface_side_get_emissivity(interf, frag);
+ res = interface_side_check_emissivity(dev, epsilon, frag->P, frag->time);
+ if(res != RES_OK) goto error;
+
+ alpha = interface_side_get_specular_fraction(interf, frag);
+ res = interface_side_check_specular_fraction(dev, alpha, frag->P, frag->time);
+ if(res != RES_OK) goto error;
+
+ brdf->emissivity = epsilon;
+ brdf->specular_fraction = alpha;
+
+exit:
+ return res;
+error:
+ *brdf = BRDF_NULL;
+ goto exit;
+}
+
+static double /* [W/m^2] */
+XD(compute_incident_diffuse_flux)
+ (const struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const double in_pos[DIM], /* position */
+ const double in_N[DIM], /* Surface normal. (Away from the surface) */
+ const double time,
+ const struct sXd(hit)* in_hit) /* Current intersection */
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ double pos[DIM] = {0};
+ double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/
+ double N[3] = {0}; /* Surface normal. Always 3D */
+ double incident_diffuse_flux = 0; /* [W/m^2] */
+ res_T res = RES_OK;
+ ASSERT(in_pos && in_N && in_hit);
+
+ /* Local copy of input argument */
+ dX(set)(pos, in_pos);
+ dX(set)(N, in_N);
+ hit = *in_hit;
+
+ /* Sample a diffusive direction in 3D */
+ ssp_ran_hemisphere_cos(rng, N, dir, NULL);
+
+ for(;;) {
+ /* External sources */
+ struct sdis_external_sources_sample extsrc_sample =
+ SDIS_EXTERNAL_SOURCES_SAMPLE_NULL;
+ sdis_sample_external_sources_T sample_sources = NULL;
+ sdis_trace_external_sources_T trace_sources = NULL;
+
+ /* Interface */
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sdis_interface* interf = NULL;
+
+ /* BRDF */
+ struct brdf brdf = BRDF_NULL;
+ struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL;
+
+ /* Miscellaneous */
+ double L = 0; /* incident direct flux to bounce position */
+ double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */
+ double vec[DIM] = {0}; /* Temporary variable */
+
+ dX(minus)(wi, dir);
+
+ /* Find the following surface along the direction of propagation */
+ XD(trace_ray)(scn, pos, dir, INF, &hit, &hit);
+ if(SXD_HIT_NONE(&hit)) break; /* No surface */
+
+ /* Retrieve the current position and normal */
+ dX(add)(pos, pos, dX(muld)(vec, dir, hit.distance));
+ dX_set_fX(N, hit.normal);
+ dX(normalize(N, N));
+
+ /* Retrieve the current interface properties */
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ XD(setup_fragment)(&frag, pos, dir, time, N, &hit);
+ XD(setup_brdf)(scn->dev, &brdf, interf, &frag);
+ sample_sources = interface_side_get_external_sources_sampling_functor
+ (interf, &frag);
+ trace_sources = interface_side_get_external_sources_tracing_functor
+ (interf, &frag);
+
+ /* Check if path is absorbed */
+ if(ssp_rng_canonical(rng) < brdf.emissivity) break;
+
+ /* Sample rebound direction */
+ if(frag.side == SDIS_BACK) dX(minus)(N, N); /* Revert normal if necessary */
+ sample_brdf(&brdf, rng, wi, N, &brdf_sample);
+ dX(set)(dir, brdf_sample.dir);
+
+ /* Calculate the direct contribution if the rebound is specular */
+ if(brdf_sample.cpnt == BRDF_SPECULAR) {
+ res = trace_sources(&frag, brdf_sample.dir, &extsrc_sample, interf->data);
+ CHK(res == RES_OK); /* TODO handle the error */
+
+ if(!SDIS_EXTERNAL_SOURCES_SAMPLE_NONE(&extsrc_sample)) {
+ const double Ld = XD(direct_contribution)(scn, &extsrc_sample, pos, &hit);
+ L = Ld; /* [W/m^2] */
+ }
+
+ /* Calculate the direct contribution of the rebound is diffuse */
+ } else {
+ ASSERT(brdf_sample.cpnt == BRDF_DIFFUSE);
+
+ /* Sample an external source to handle its direct contribution at the
+ * bounce position */
+ sample_sources(&frag, rng, &extsrc_sample, interf->data);
+ CHK(res == RES_OK); /* TODO handle the error */
+
+ /* The source is behind the surface */
+ if(d3_dot(extsrc_sample.dir, N) <= 0) {
+ L = 0; /* [W/m^2] */
+
+ /* The source is above the surface */
+ } else {
+ const double Ld = XD(direct_contribution)(scn, &extsrc_sample, pos, &hit);
+ const double cos_theta = d3_dot(extsrc_sample.dir, N);
+ L = Ld * cos_theta/PI * extsrc_sample.pdf; /* [W/m^2] */
+ }
+ }
+ incident_diffuse_flux += L;
+ }
+
+ incident_diffuse_flux *= PI;
+ return incident_diffuse_flux;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -65,20 +350,23 @@ XD(handle_external_net_flux)
struct XD(temperature)* T)
{
/* Sampling external sources */
- struct sdis_external_sources_sample sample = SDIS_EXTERNAL_SOURCES_SAMPLE_NULL;
+ struct sdis_external_sources_sample extsrc_sample =
+ SDIS_EXTERNAL_SOURCES_SAMPLE_NULL;
sdis_sample_external_sources_T sample_sources = NULL;
- /* Ray tracing */
- struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
- struct sXd(hit) hit = SXD_HIT_NULL;
- float ray_org[DIM] = {0};
- float ray_dir[DIM] = {0};
- float ray_range[2] = {0};
-
/* External flux */
- double incident_direct_flux = 0; /* [W/m^2/sr] TODO ? */
+ double incident_flux = 0; /* [W/m^2] */
+ double incident_flux_diffuse = 0; /* [W/m^2] */
+ double incident_flux_direct = 0; /* [W/m^2] */
+ double net_flux = 0; /* [W/m^2] */
+
+ /* Sampled path */
+ double N[3] = {0}; /* Normal. Always in 3D */
/* Miscellaneous */
+ double emissivity = 0; /* Emissivity */
+ double Ld = 0; /* Incident radiance [W/m^2/sr] */
+ double cos_theta = 0;
res_T res = RES_OK;
ASSERT(scn && args && T);
@@ -93,31 +381,34 @@ XD(handle_external_net_flux)
if(!sample_sources) goto exit;
/* Sample an external sources */
- res = sample_sources(args->frag, rng, &sample, args->interf->data);
- if(res != RES_OK) {
- log_err(scn->dev,
- "%s: error when sampling external sources -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
+ res = sample_sources(args->frag, rng, &extsrc_sample, args->interf->data);
+ if(res != RES_OK) goto error;
- /* Check whether the sampled external source is occluded or not */
- filter_data.XD(hit) = *args->hit;
- fX_set_dX(ray_org, args->frag->P);
- fX_set_dX(ray_dir, sample.dir);
- ray_range[0] = 0;
- ray_range[1] = (float)sample.distance;
- SXD(scene_view_trace_ray
- (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, &hit));
+ /* Local path data */
+ dX(set)(N, args->frag->Ng);
+ if(args->frag->side == SDIS_BACK) dX(minus)(N, N);
- /* If the source is not occluded, compute its direct contribution */
- if(!SXD_HIT_NONE(&hit)) {
- const double cos_theta = dX(dot)(args->frag->Ng, sample.dir);
- incident_direct_flux = fabs(cos_theta) * sample.radiance / sample.pdf;
- }
+ /* Calculate the incident direct flux */
+ Ld = XD(direct_contribution)(scn, &extsrc_sample, args->frag->P, args->hit);
+ cos_theta = d3_dot(N, extsrc_sample.dir);
+ incident_flux_direct = cos_theta * Ld / extsrc_sample.pdf; /* [W/m^2] */
+
+ /* Calculate the incident diffuse flux [W/m^2] */
+ incident_flux_diffuse = XD(compute_incident_diffuse_flux)
+ (scn, rng, args->frag->P, N, args->frag->time, args->hit);
+
+ /* Calculate the global incident flux */
+ incident_flux = incident_flux_direct + incident_flux_diffuse; /* [W/m^2] */
+
+ /* Calculate the net flux */
+ emissivity = interface_side_get_emissivity(args->interf, args->frag);
+ res = interface_side_check_emissivity
+ (scn->dev, emissivity, args->frag->P, args->frag->time);
+ if(res != RES_OK) goto error;
+ net_flux = incident_flux * emissivity; /* [W/m^2] */
- /* TODO handle diffuse contribution */
- /* TODO update the weight */
+ /* Update the Monte Carlo weight */
+ T->value += net_flux / (args->h_radi + args->h_conv + args->h_cond);
exit:
return res;
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -180,6 +180,9 @@ struct handle_external_net_flux_args_2d {
struct sdis_heat_path* heat_path; /* Save paths */
size_t picard_order;
+ double h_cond; /* Convective coefficient, i.e. lambda/delta */
+ double h_conv; /* Condutive coefficient */
+ double h_radi; /* Radiative coefficient */
};
struct handle_external_net_flux_args_3d {
@@ -191,10 +194,13 @@ struct handle_external_net_flux_args_3d {
struct sdis_heat_path* heat_path; /* Save paths */
size_t picard_order;
+ double h_cond; /* Convective coefficient, i.e. lambda/delta */
+ double h_conv; /* Condutive coefficient */
+ double h_radi; /* Radiative coefficient */
};
-#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___2d {NULL,NULL,NULL,NULL,NULL,0}
-#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___3d {NULL,NULL,NULL,NULL,NULL,0}
+#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___2d {NULL,NULL,NULL,NULL,NULL,0,0,0,0}
+#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___3d {NULL,NULL,NULL,NULL,NULL,0,0,0,0}
static const struct handle_external_net_flux_args_2d
HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL_2d = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___2d;
static const struct handle_external_net_flux_args_3d
