star-phor

Radiative transfer solver for photoreactors.
git clone https://www.edstar.cnrs.fr/git/star-phor.git
Log | Files | Refs | README | LICENSE

commit fbfd67ead297ae300eccba37297e738a1b98c207
parent bc76ac9eb5f4857790fe903308a467459115e066
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Thu, 27 Nov 2025 19:41:02 +0100

Integrate geometrical optics into the MVREA algorithm

This commit builds on the input file and star-phor-input evolution
(which is now able to handle refractive indices) and integrates it into
the Monte Carlo computation when a photon intersects an interface.

As BRDFs were already supported in the input file, star-phor
prioritizes user-defined BRDFs. Possible outcomes for ray–interface
interactions with a user-defined BRDF are either reflection or
absorption (the latter stopping the path).

If no BRDF is defined for a given interface, star-phor will then try
to retrieve optical indices from each side of the interface and compute
the ray–interface interaction outcome using geometrical optics (Fresnel
and Snell–Descartes laws). This is done using the
star-scattering-functions library, which is added as a new dependency
of star-phor in this commit.

If no refractive index is set for a given volume, default values of 1
for the real part and 0 for the imaginary part are set.

Diffstat:
MMakefile | 2++
MREADME.md | 1+
Mconfig.mk | 5+++--
Msrc/sphor_compute_mvrea.c | 424++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------
4 files changed, 364 insertions(+), 68 deletions(-)

diff --git a/Makefile b/Makefile @@ -70,6 +70,7 @@ $(LIBNAME): $(OBJ_LIB) $(PKG_CONFIG) --atleast-version $(SPHIN_VERSION) sphin $(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp $(PKG_CONFIG) --atleast-version $(SUNIQ_VERSION) suniq + $(PKG_CONFIG) --atleast-version $(SSF_VERSION) ssf echo "config done" > $@ .SUFFIXES: .c .d .o @@ -127,6 +128,7 @@ pkg: -e 's#@SPHIN_VERSION@#$(SPHIN_VERSION)#g'\ -e 's#@SSP_VERSION@#$(SSP_VERSION)#g'\ -e 's#@SUNIQ_VERSION@#$(SUNIQ_VERSION)#g'\ + -e 's#@SSF_VERSION@#$(SSF_VERSION)#g'\ sphor.pc.in > sphor.pc star-phor.1: doc/star-phor.scd diff --git a/README.md b/README.md @@ -12,6 +12,7 @@ Solver for radiative transfer in photoreactors. - star-3d - star-sp - star-uniq +- star-sf - Open MPI ## Installation diff --git a/config.mk b/config.mk @@ -29,10 +29,11 @@ S3D_VERSION=0.10 SPHIN_VERSION=0.0 SSP_VERSION=0.14 SUNIQ_VERSION=0.0 +SSF_VERSION=0.10 -INCS = $$($(PKG_CONFIG) --cflags $(MPI_PC) rsys s3d sphin star-sp suniq)\ +INCS = $$($(PKG_CONFIG) --cflags $(MPI_PC) rsys s3d sphin star-sp suniq ssf)\ -fopenmp -lm -LIBS = $$($(PKG_CONFIG) --libs $(MPI_PC) rsys s3d sphin star-sp suniq)\ +LIBS = $$($(PKG_CONFIG) --libs $(MPI_PC) rsys s3d sphin star-sp suniq ssf)\ -fopenmp -lm ################################################################################ diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -33,6 +33,7 @@ #include <star/sphin.h> #include <star/s3d.h> #include <star/ssp.h> +#include <star/ssf.h> #include <rsys/dynamic_array.h> #include <rsys/hash_table.h> #include <omp.h> @@ -191,25 +192,60 @@ error: } static res_T +prop_rad_compute_ka + (struct sphin_prop_rad* prop_rad, + double wavelength, + double* ka) +{ + struct sphin_scatterer* scatterer = NULL; + struct sphin_spectral_property* abs_cross_sec = NULL; + double concentration = 0; + double sigma_a = 0; + res_T res = RES_OK; + + res = sphin_prop_rad_get_scatterer(prop_rad, &scatterer); + if (RES_OK != res){ goto error; } + + res = sphin_scatterer_get_concentration(scatterer, &concentration); + if (RES_OK != res){ goto error; } + + res = sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec); + if (RES_OK != res){ goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (abs_cross_sec, wavelength, SPHIN_INTERPOLATION_LINEAR, &sigma_a); + if (RES_OK != res){ goto error; } + + *ka = concentration * sigma_a; +exit: + return res; +error: + goto exit; +} + +static res_T volume_compute_total_ka (struct sphor* sphor, struct interface interface, size_t side_id, double wavelength, - double* ka) + double* out_ka) { - double sigma_a = 0; - double concentration = 0; size_t prop_rad_count = 0; size_t j = 0; + double ka = 0; struct sphin_prop_rad* prop_rad = NULL; - struct sphin_scatterer* scatterer = NULL; - struct sphin_spectral_property* abs_cross_sec = NULL; struct sphin_volume* volume = NULL; res_T res = RES_OK; ASSERT(NULL != sphor); - ASSERT(NULL != ka); + ASSERT(NULL != out_ka); + + /* No volume defined in the interface side */ + if (INVALID_ID == interface.volumes[side_id]) { + ka = 0.; + goto exit; + } res = sphin_config_get_volume (sphor->config, interface.volumes[side_id], &volume); @@ -219,29 +255,19 @@ volume_compute_total_ka res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); if (RES_OK != res){ goto error; } - *ka = 0.; - FOR_EACH(j, 0, prop_rad_count) { + double ka_i = 0; res = sphin_volume_get_prop_rad(volume, j, &prop_rad); if (RES_OK != res){ goto error; } - res = sphin_prop_rad_get_scatterer(prop_rad, &scatterer); - if (RES_OK != res){ goto error; } - - res = sphin_scatterer_get_concentration(scatterer, &concentration); + res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); if (RES_OK != res){ goto error; } - res = sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec); - if (RES_OK != res){ goto error; } - - res = sphin_spectral_property_interpolate_at_wavelength - (abs_cross_sec, wavelength, SPHIN_INTERPOLATION_LINEAR, &sigma_a); - if (RES_OK != res){ goto error; } - - *ka += concentration * sigma_a; + ka += ka_i; } exit: + *out_ka = ka; return res; error: goto exit; @@ -253,28 +279,102 @@ sample_abs_free_path_exp struct ssp_rng* rng, struct interface interface, size_t side_id, - double wavelength, + double ka, double* s) { - double ka = 0; /* Absorption coefficient of current volume */ - res_T res = RES_OK; ASSERT(NULL != sphor); ASSERT(NULL != s); + (void)sphor; + /* Get the absorption coefficient of the current volume to compute an * absorption free length. If the intersected interface does not correspond * to any volume in the config, the absorption coefficient is supposed to be * zero and the absorption free length is set to FLT_MAX */ if (INVALID_ID != interface.volumes[side_id]) { - res = volume_compute_total_ka(sphor, interface, side_id, wavelength, &ka); - if (RES_OK != res){ goto error; } /* Sample free path length according to volume absorption coefficent */ *s = ssp_ran_exp(rng, ka);} else { *s = FLT_MAX; } + return RES_OK; +} + +static res_T +interface_side_get_refractive_index + (struct sphor* sphor, + struct interface* interface, + enum sphin_side side_id, + double wavelength, + double* n_real, + double* n_imag) +{ + struct sphin_volume* volume = NULL; + struct sphin_refractive_index* refr_ind = NULL; + struct sphin_spectral_property* spec_prop = NULL; + struct sphin_spectral_property_descriptor spec_prop_desc + = SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != interface); + ASSERT(NULL != n_real); + ASSERT(NULL != n_imag); + + /* If the interface side doesnt has a volume defined, the refractive index + * receives the default value of 1 to the real part and 0 to the imaginary + * part */ + if (INVALID_ID == interface->volumes[side_id]) { + *n_real = 1; + *n_imag = 0; + goto exit; + } + + res = sphin_config_get_volume + (sphor->config, interface->volumes[side_id], &volume); + if (RES_OK != res){ goto error; } + + res = sphin_volume_get_refractive_index(volume, &refr_ind); + if (RES_OK != res){ goto error; } + + /* If the volume is null, the refractive index + * receives the default value of 1 to the real part and 0 to the imaginary + * part */ + if (NULL == refr_ind) { + *n_real = 1; + *n_imag = 0; + goto exit; + } + + /* Real part of the refractive index */ + res = sphin_refractive_index_get_n_real(refr_ind, &spec_prop); + if (RES_OK != res){ goto error; } + + if (NULL != spec_prop) { + res = sphin_spectral_property_get_desc(spec_prop, &spec_prop_desc); + if (RES_OK != res){ goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (spec_prop, wavelength, SPHIN_INTERPOLATION_LINEAR, n_real); + if (RES_OK != res){ goto error; } + } + else { *n_real = 1; } + + /* Imag part of the refractive index */ + res = sphin_refractive_index_get_n_imag(refr_ind, &spec_prop); + if (RES_OK != res){ goto error; } + + if (NULL != spec_prop) { + res = sphin_spectral_property_get_desc(spec_prop, &spec_prop_desc); + if (RES_OK != res){ goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (spec_prop, wavelength, SPHIN_INTERPOLATION_LINEAR, n_imag); + if (RES_OK != res){ goto error; } + } + else { *n_imag = 0; } exit: return res; error: @@ -282,70 +382,253 @@ error: } static res_T -interface_sample_ray_interaction_type +interface_sample_ray_interaction_type_brdf (struct sphor* sphor, struct ssp_rng* rng, struct interface* interface, enum sphin_side side_id, + struct s3d_hit* hit, + double dir[3], + double wavelength, enum ray_interface_interaction_type* interaction_type) { enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__; - double reflectivity = 0; - double s = 0; /* random number */ struct sphin_surface* surface = NULL; struct sphin_brdf* brdf = NULL; + double reflectivity = 0; + double s = 0; /* random number */ size_t i = 0; res_T res = RES_OK; ASSERT(NULL != sphor); + ASSERT(NULL != rng); + ASSERT(NULL != hit); ASSERT(NULL != interface); ASSERT(NULL != interaction_type); - FOR_EACH(i, 0, interface->surface_count[(size_t)side_id]) { + (void)hit; + (void)dir; + (void)wavelength; + + FOR_EACH(i, 0, (size_t)interface->surface_count[(size_t)side_id]) { res = sphin_config_get_surface (sphor->config, interface->surfaces[(size_t)side_id][i], &surface); if (RES_OK != res){ goto error; } res = sphin_surface_get_brdf(surface, &brdf); if (RES_OK != res){ goto error; } - res = sphin_brdf_get_type(brdf, &brdf_type); - if (RES_OK != res){ goto error; } if (NULL != brdf){ break; } } - /* If the interface side does not have a BRDF (i.e., it defines only a - * sphin_volume or a sphin_surface without a BRDF), the photon - * crosses the interface the path continues */ - if((NULL == brdf) || (interface->surface_count[side_id] == 0)){ - *interaction_type = RAY_INTERFACE_INTERACTION_TRANSMISSION; + if (NULL != brdf) { res = RES_BAD_ARG; goto error; } + + res = sphin_brdf_get_type(brdf, &brdf_type); + if (RES_OK != res){ goto error; } + + switch (brdf_type) { + case SPHIN_BRDF_LAMBERT: + res = sphin_brdf_lambertian_get_reflectivity(brdf, &reflectivity); + if (RES_OK != res){ goto error; } + break; + case SPHIN_BRDF_SPECULAR: + res = sphin_brdf_specular_get_reflectivity(brdf, &reflectivity); + if (RES_OK != res){ goto error; } + break; + case SPHIN_BRDF_NONE__: + res = RES_BAD_ARG; + goto error; + default: + FATAL("Unreachable code\n"); + break; } - /* If the interface side is a sphin_surface with a BRDF declared, - * verify if a reflection takes place based on the reflectivity of the - * surface */ - else{ - switch (brdf_type) { - case SPHIN_BRDF_LAMBERT: - res = sphin_brdf_lambertian_get_reflectivity(brdf, &reflectivity); - if (RES_OK != res){ goto error; } + /* Bernoulli test if the photon was absorbed or reflected */ + s = ssp_rng_canonical(rng); + if (s < reflectivity) { + *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; + } + else { *interaction_type = RAY_INTERFACE_INTERACTION_ABSORPTION; } + +exit: + return res; +error: + goto exit; +} + +static res_T +sample_interface_ray_interaction_from_brdf +(struct sphor* sphor, + struct ssp_rng* rng, + struct interface* interface, + enum sphin_side side_id, + struct s3d_hit* hit, + double wavelength, + double dir[3], + enum ray_interface_interaction_type* interaction_type) +{ + struct sphin_surface* surface = NULL; + struct sphin_brdf* brdf = NULL; + size_t i = 0; + double sample[3] = {0}; + res_T res = RES_OK; + + FOR_EACH(i, 0, (size_t)interface->surface_count[(size_t)side_id]) { + res = sphin_config_get_surface + (sphor->config, interface->surfaces[(size_t)side_id][i], &surface); + if (RES_OK != res){ goto error; } + res = sphin_surface_get_brdf(surface, &brdf); + if (RES_OK != res){ goto error; } + if (NULL != brdf){ break; } + } + + res = interface_sample_ray_interaction_type_brdf + (sphor, rng, interface, side_id, hit, dir, wavelength, interaction_type); + if (RES_OK != res){ goto error; } + + switch (*interaction_type) { + case RAY_INTERFACE_INTERACTION_ABSORPTION: break; - case SPHIN_BRDF_SPECULAR: - res = sphin_brdf_specular_get_reflectivity(brdf, &reflectivity); - if (RES_OK != res){ goto error; } + case RAY_INTERFACE_INTERACTION_REFLECTION: + res = ran_brdf_reflection_direction(brdf, rng, &hit->prim, dir,sample); + d3_set(dir, sample); + break; + case RAY_INTERFACE_INTERACTION_TRANSMISSION: break; - case SPHIN_BRDF_NONE__: + case RAY_INTERFACE_INTERACTION_NONE__: res = RES_BAD_ARG; goto error; default: FATAL("Unreachable code\n"); break; } - /* Test if the photon was absorbed or reflected */ - s = ssp_rng_canonical(rng); - if (s < reflectivity) { - *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; - } - else { *interaction_type = RAY_INTERFACE_INTERACTION_ABSORPTION; } + +exit: + return res; +error: + goto exit; +} + +static res_T +sample_interface_ray_interaction_from_refractive_index +(struct sphor* sphor, + struct ssp_rng* rng, + struct interface* interface, + enum sphin_side side_id, + struct s3d_hit* hit, + double wavelength, + double dir[3], + enum ray_interface_interaction_type* interaction_type) +{ + double invert_normal = 1.; + double normal[3] = {0}; + double n_real_i = 0, n_imag_i = 0; + double n_real_t = 0, n_imag_t = 0; + double pdf = 0; + double sample[3] = {0}; + int flag = 0; + size_t side_id_back = 0; + struct s3d_attrib attrib; + struct ssf_bsdf* bsdf = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != rng); + ASSERT(NULL != hit); + ASSERT(NULL != interface); + ASSERT(NULL != dir); + ASSERT(NULL != interaction_type); + + /* Get refractive index of the medium in the incident side */ + res = interface_side_get_refractive_index + (sphor, interface, side_id, wavelength, &n_real_i, &n_imag_i); + if (RES_OK != res){ goto error; } + + /* Get refractive index of the medium in the transmitted side */ + side_id_back = SPHIN_SIDE_FRONT^side_id; + res = interface_side_get_refractive_index + (sphor, interface, side_id_back, wavelength, &n_real_t, &n_imag_t); + if (RES_OK != res){ goto error; } + + /* Get surface normal in the intersection point and ensure its compatible with + * the ssf API */ + res = s3d_primitive_get_attrib + (&hit->prim, S3D_GEOMETRY_NORMAL, hit->uv, &attrib); + if (RES_OK != res){ goto error; } + d3_set_f3(normal, attrib.value); + /* Ensure normal and incoming direction point to the same hemisphere */ + invert_normal = (d3_dot(normal, dir) > 0 ) ?1. : -1.; + d3_normalize(normal, d3_muld(normal, normal, invert_normal)); + + /* Sample an interaction type and a new direction using ssf */ + res = ssf_bsdf_create + (sphor->allocator, &ssf_specular_dielectric_dielectric_interface, &bsdf); + if (RES_OK != res){ goto error; } + + res = ssf_specular_dielectric_dielectric_interface_setup + (bsdf, n_real_i, n_real_t); + if (RES_OK != res){ goto error; } + ssf_bsdf_sample(bsdf, rng, dir, normal, sample, &flag, &pdf); + + if (flag && SSF_TRANSMISSION) { + *interaction_type = RAY_INTERFACE_INTERACTION_TRANSMISSION; } + else if (flag && SSF_REFLECTION) { + *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; } + else { res = RES_BAD_ARG; } + +exit: + if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } + return res; +error: + goto exit; +} + +static res_T +sample_interface_ray_interaction + (struct sphor* sphor, + struct ssp_rng* rng, + struct interface* interface, + enum sphin_side side_id, + struct s3d_hit* hit, + double wavelength, + double dir[3], + enum ray_interface_interaction_type* interaction_type) +{ + struct sphin_surface* surface = NULL; + struct sphin_brdf* brdf = NULL; + size_t i = 0; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != rng); + ASSERT(NULL != hit); + ASSERT(NULL != interface); + + /* Verify if one of the surfaces defined in the intersected interface + * side has a BRDF defined */ + FOR_EACH(i, 0, (size_t)interface->surface_count[(size_t)side_id]) { + res = sphin_config_get_surface + (sphor->config, interface->surfaces[(size_t)side_id][i], &surface); + if (RES_OK != res){ goto error; } + res = sphin_surface_get_brdf(surface, &brdf); + if (RES_OK != res){ goto error; } + if (NULL != brdf){ break; } + } + + if(NULL != brdf) { + /* If the interface side is a sphin_surface with a BRDF declared, + * verify if a reflection takes place based on the reflectivity of the + * surface */ + res = sample_interface_ray_interaction_from_brdf + (sphor, rng, interface, side_id, hit, wavelength, dir, interaction_type); + if (RES_OK != res){ goto error; } } + /* If the interface side does not have a BRDF (i.e., it defines only a + * sphin_volume or a sphin_surface without a BRDF), we use the refractive + * index of each side of the interface to construct the BRDF */ + else{ + res = sample_interface_ray_interaction_from_refractive_index + (sphor, rng, interface, side_id, hit, wavelength, dir, interaction_type); + if (RES_OK != res){ goto error; } + } exit: return res; error: @@ -380,8 +663,6 @@ trace_ray if (S3D_HIT_NONE(&hit)) { goto exit; } - - /* TODO: separated function */ /* Retrieve the intercepted triangle and the corresponding interface in the * sphor structure */ triangle_id = hit.prim.prim_id; @@ -739,9 +1020,9 @@ compute_MVREA_realization double dir[3] = {0}; /* Current direction in path sampling */ double pos[3] = {0}; /* Current position in path sampling */ double response_function = 0; - double s = 0; /* Random number */ - double sample[3] = {0}; double wavelength = 0; + double ka = 0; + double s = 0; /* Random number */ enum ray_interface_interaction_type ray_interface_interaction_type = RAY_INTERFACE_INTERACTION_NONE__; enum sphin_side side_id = SPHIN_SIDE_NONE__; @@ -751,13 +1032,14 @@ compute_MVREA_realization struct s3d_primitive prim = S3D_PRIMITIVE_NULL; struct s3d_attrib attrib; struct source_view* source_view = NULL; - struct sphin_brdf* brdf = NULL; struct sphin_sensor_volume* sensor_volume = NULL; struct sphin_volume* volume = NULL; res_T res = RES_OK; ASSERT(NULL != sphor); ASSERT(NULL != rng); + ASSERT(NULL != accum2id); + ASSERT(NULL != accums); /* Sample a random position from a source in the whole scene */ res = sample_source_position(sphor, rng, &source_view, &prim, pos); @@ -782,11 +1064,14 @@ compute_MVREA_realization * position, the photon is lost and it counts zero in the weight. */ if (S3D_HIT_NONE(&hit)) { break; } - /* TODO: Retrive prop rads of current volume */ + /* Retrieve the ka of the present volume */ + res = volume_compute_total_ka(sphor, interface, side_id, wavelength, &ka); + if (RES_OK != res){ goto error; } + /* Sample the free path length to absorption from an exponential * distribution with rate parameter ka */ res = sample_abs_free_path_exp - (sphor, rng, interface, side_id, wavelength, &s); + (sphor, rng, interface, side_id, ka, &s); if (RES_OK != res){ goto error; } /* If the distance between the ray origin and the intersection is bigger @@ -824,8 +1109,15 @@ compute_MVREA_realization * than the absorption free path length. The photon intersects a primitive * in the scene view. Sample the type of interaction that the photon has * with the interface */ - res = interface_sample_ray_interaction_type - (sphor, rng, &interface, side_id, &ray_interface_interaction_type); + res = sample_interface_ray_interaction + (sphor, + rng, + &interface, + side_id, + &hit, + wavelength, + dir, + &ray_interface_interaction_type); if (RES_OK != res){ goto error; } int stop = 0; @@ -839,13 +1131,13 @@ compute_MVREA_realization * the photon arriving direction, geometry normak and surface brdf * properties */ case RAY_INTERFACE_INTERACTION_REFLECTION: - res = ran_brdf_reflection_direction(brdf, rng, &hit.prim, dir, sample); - if (RES_OK != res){ goto error; } - d3_set(dir, sample); break; /* If the ray traverses the interface, it simply continue in the same * direction as before */ case RAY_INTERFACE_INTERACTION_TRANSMISSION: + /* The position coordinates are the same, but on the other side of the + * interface */ + side_id = SPHIN_SIDE_FRONT ^ side_id; /* Cross the interface */ break; /* Handle the possible exceptions */ case RAY_INTERFACE_INTERACTION_NONE__: