star-phor

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

commit dc7eb4e681146cc588ae724c857f3ce6fb3fc964
parent 94a5f9e25a3f1859ec2113a70bf5fe17cdcac3f9
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Wed, 28 Jan 2026 15:59:36 +0100

Extract weight update logic into a separate function

Weight updates were previously implemented directly in the realization
function, leading to two main issues that hinder readability: (i) the
update logic spans many lines, obscuring the intent of the surrounding
block and disrupting the flow of the function; (ii) lower-level
implementation details appear at a higher level of abstraction.

Diffstat:
Msrc/sphor_compute_mvrea.c | 238++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------
1 file changed, 168 insertions(+), 70 deletions(-)

diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -643,6 +643,157 @@ error: } static res_T +MVREA_update_volume_weights + (struct sphor* sphor, + struct intersection* intersection, + double wavelength, + struct htable_accum2id* accum2id, + struct darray_accum* accums) +{ + struct accum_key accum_key = ACCUM_KEY_NULL; + size_t prop_rad_count = 0; + size_t volume_id = SIZE_MAX; + size_t i_proprad = 0; + double response_function = 0; + double vol_size = 0; + double tot_pow = 0; + double ka = 0; + double mc_weight = 0; + struct sphin_volume* volume = NULL; + struct sphin_sensor_volume* sensor_volume = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != intersection); + ASSERT(NULL != accum2id); + ASSERT(NULL != accums); + + volume = volume_get_from_primitive(sphor, &intersection->position.primitive); + + /* Get the total three-dimensional space occupied by the volume */ + res = sphin_volume_compute_total_size(volume, &vol_size/*[m^3]*/); + if (RES_OK != res) { goto error; } + res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); + if (RES_OK != res) { goto error; } + + mc_weight = tot_pow / vol_size; + + res = sphin_volume_get_sensor(volume, &sensor_volume); + if (RES_OK != res) { goto error; } + + res = sphin_sensor_volume_get_response_function + (sensor_volume, &response_function); + if (RES_OK != res) { goto error; } + res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); + if (RES_OK != res) { goto error; } + + volume_id = intersection->position.primitive.interface->volumes + [intersection->position.primitive.side]; + + /* Update total accumulator for all volumes in the scene*/ + accum_key.sensor_type = SPHOR_SENSOR_VOLUME; + accum_key.sensor_id = ALL_SENSORS; + accum_key.component_id = ALL_COMPONENTS; + update_accum(accum2id, &accum_key, mc_weight, accums); + + /* Update total accumulator in the volume */ + accum_key.sensor_type = SPHOR_SENSOR_VOLUME; + accum_key.sensor_id = volume_id; + accum_key.component_id = ALL_COMPONENTS; + update_accum(accum2id, &accum_key, mc_weight, accums); + + /* Increment the weight of each chemical species in the volume by + * ( ka_i / ka_tot ) response_function */ + res = volume_compute_total_ka(sphor, volume, wavelength, &ka); + if (RES_OK != res) { goto error; } + FOR_EACH(i_proprad, 0, prop_rad_count) { + double ka_i = 0; + struct sphin_prop_rad* prop_rad = NULL; + res = sphin_volume_get_prop_rad(volume, i_proprad, &prop_rad); + if (RES_OK != res) { goto error; } + res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); + if (RES_OK != res) { goto error; } + accum_key.component_id = i_proprad; + + mc_weight = ka_i / ka * tot_pow / vol_size; + update_accum(accum2id, &accum_key, mc_weight, accums); + } + +exit: + return res; +error: + goto exit; +} + +static res_T +MVREA_update_surface_weights + (struct sphor* sphor, + struct intersection* intersection, + double wavelength, + struct htable_accum2id* accum2id, + struct darray_accum* accums) +{ + struct accum_key accum_key = ACCUM_KEY_NULL; + size_t surface_id = SIZE_MAX; + size_t i = 0; + double response_function = 0; + double tot_pow = 0; + double mc_weight = 0; + double surf_area = 0; + struct sphin_surface* surface = NULL; + struct sphin_sensor_surface* sensor_surface = NULL; + struct primitive prim = PRIMITIVE_NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != intersection); + ASSERT(NULL != accum2id); + ASSERT(NULL != accums); + + (void)wavelength; + + prim = intersection->position.primitive; + + /* Retrieve the sensor and the surface_id corresponding to the sensor */ + FOR_EACH(i, 0, (size_t)prim.interface->surface_count[(size_t)prim.side]) { + surface_id = (size_t)prim.interface->surfaces[prim.side][i]; + res = sphin_config_get_surface(sphor->config, surface_id, &surface); + if (RES_OK != res) { goto error; } + res = sphin_surface_get_sensor(surface, &sensor_surface); + if (RES_OK != res) { goto error; } + + if (NULL != sensor_surface) { break; } + } + + if (NULL == sensor_surface) { res = RES_BAD_ARG; goto error; } + + res = sphin_surface_compute_total_area(surface, &surf_area); + if (RES_OK != res) { goto error; } + res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); + if (RES_OK != res) { goto error; } + res = sphin_sensor_surface_get_response_function + (sensor_surface, &response_function); + if (RES_OK != res) { goto error; } + + mc_weight = response_function * tot_pow / surf_area; + + /* Update total accumulator for all surfaces losses in the scene*/ + accum_key.sensor_type = SPHOR_SENSOR_SURFACE; + accum_key.sensor_id = ALL_SENSORS; + update_accum(accum2id, &accum_key, mc_weight, accums); + + /* Update total accumulator in the surface */ + accum_key.sensor_type = SPHOR_SENSOR_SURFACE; + accum_key.sensor_id = surface_id; + update_accum(accum2id, &accum_key, mc_weight, accums); + +exit: + return res; +error: + goto exit; +} + +static res_T MVREA_write_outputs (struct sphor* sphor, const struct darray_accum* accums, @@ -690,7 +841,6 @@ compute_MVREA_realization /* Miscellaneous */ double dir[3] = {0}; /* Sampled direction during path sampling */ - double response_function = 0; double ka = 0; double free_path = 0; double wavelength = 0; @@ -756,91 +906,40 @@ compute_MVREA_realization } else { /* If the volume is a sensor, compute the weight using the response * function */ - - /* TODO res = update_mc_weights_MVREA_volume(...,mc_weight,...); */ - struct accum_key accum_key = ACCUM_KEY_NULL; - size_t prop_rad_count = 0; - size_t volume_id = SIZE_MAX; - size_t i_proprad = 0; - double vol_size = 0; - double tot_pow = 0; - double mc_weight = 0; - - /* Get the total three-dimensional space occupied by the volume */ - res = sphin_volume_compute_total_size(volume, &vol_size/*[m^3]*/); - if (RES_OK != res) { goto error; } - res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); - if (RES_OK != res) { goto error; } - - mc_weight = tot_pow / vol_size; - res = sphin_sensor_volume_get_response_function - (sensor_volume, &response_function); - if (RES_OK != res) { goto error; } - res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); - if (RES_OK != res) { goto error; } - - volume_id = intersection.position.primitive.interface->volumes - [intersection.position.primitive.side]; - - /* Update total accumulator for all volumes in the scene*/ - accum_key.sensor_type = SPHOR_SENSOR_VOLUME; - accum_key.sensor_id = ALL_SENSORS; - accum_key.component_id = ALL_COMPONENTS; - update_accum(accum2id, &accum_key, mc_weight, accums); - - /* Update total accumulator in the volume */ - accum_key.sensor_type = SPHOR_SENSOR_VOLUME; - accum_key.sensor_id = volume_id; - accum_key.component_id = ALL_COMPONENTS; - update_accum(accum2id, &accum_key, mc_weight, accums); - - /* Increment the weight of each chemical species in the volume by - * ( ka_i / ka_tot ) response_function */ - FOR_EACH(i_proprad, 0, prop_rad_count) { - double ka_i = 0; - struct sphin_prop_rad* prop_rad = NULL; - res = sphin_volume_get_prop_rad(volume, i_proprad, &prop_rad); - if (RES_OK != res) { goto error; } - res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); - if (RES_OK != res) { goto error; } - accum_key.component_id = i_proprad; - - mc_weight = ka_i / ka * tot_pow / vol_size; - update_accum(accum2id, &accum_key, mc_weight, accums); - } + CALL(MVREA_update_volume_weights(sphor, + /* in */ &intersection, wavelength, accum2id, + /* out */ accums)); } - /* Absorption => stop the path */ + /* Absorption by a volume => stop the path */ stop = 1; - } else { /* Else, the distance between the ray origin and the intersection is * smaller 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 as well as the new direction sampled */ CALL(sample_interface_ray_interaction(sphor, - /* in: */ rng, &ray, &intersection, wavelength, - /* out: */ dir, &ray_interface_interaction_type)); + /* in: */ rng, &ray, &intersection, wavelength, + /* out: */ dir, &ray_interface_interaction_type)); if (RAY_INTERFACE_INTERACTION_ABSORPTION == - ray_interface_interaction_type) - { - /* TODO: verify if surface is sensor and update its accum - res = primitive_get_sensor_surface(sphor, - &intersection.position.primitive, &sensor_surface); if (RES_OK != - res){ goto error; } */ + ray_interface_interaction_type) { + CALL(primitive_get_sensor_surface + (sphor, &intersection.position.primitive, &sensor_surface)); + /* If the volume is a sensor, compute the weight using the response * function */ - if (NULL != sensor_surface){ - /* res = update_mc_weights_surface_MVREA(...= */ + if (NULL != sensor_surface) { + CALL(MVREA_update_surface_weights(sphor, + /* in */ &intersection, wavelength, accum2id, + /* out */ accums)); } - /* Absorption => stop the path */ + /* Absorption by a surface => stop the path */ stop = 1; } else if (RAY_INTERFACE_INTERACTION_TRANSMISSION == - ray_interface_interaction_type) - { + ray_interface_interaction_type) { /* Continue the path on the other side of the primitive.*/ intersection.position.primitive.side = !intersection.position.primitive.side; @@ -854,8 +953,7 @@ compute_MVREA_realization CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); } else if (RAY_INTERFACE_INTERACTION_REFLECTION == - ray_interface_interaction_type) - { + ray_interface_interaction_type) { CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); } }