stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit 5d1344b82cd16e708a559b408ab1f893282959cd
parent b6a5a6f7317f939a70e3a346d8913e489ec87ff7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 17 Apr 2025 18:32:12 +0200

Further improve the robustness of the radiative paths

Filter intersections on primitive edges if the intersected side does not
look towards the current enclosure. The intersection should belong to
another primitive which shares the same edge, but which this time looks
towards the current enclosure.

Diffstat:
Msrc/sdis_heat_path.h | 12++++++++----
Msrc/sdis_heat_path_boundary_Xd_handle_external_net_flux.h | 29+++++++++++++++++++----------
Msrc/sdis_heat_path_boundary_c.h | 4++--
Msrc/sdis_heat_path_radiative_Xd.h | 14+++++++++-----
4 files changed, 38 insertions(+), 21 deletions(-)

diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -329,41 +329,45 @@ radiative_path_3d extern LOCAL_SYM void trace_ray_2d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double pos[2], const double dir[3], /* Always in 3D */ const double distance, + const unsigned enc_id, const struct s2d_hit* hit_from, struct s2d_hit* hit); extern LOCAL_SYM void trace_ray_3d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double pos[3], const double dir[3], /* Always in 3D */ const double distance, + const unsigned enc_id, const struct s3d_hit* hit_from, struct s3d_hit* hit); /* Trace a ray and setup the fragment at the intersection found, if any. */ extern LOCAL_SYM res_T find_next_fragment_2d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double in_pos[2], const double in_dir[3], /* Always in 3D */ const struct s2d_hit* in_hit, const double time, + const unsigned enc_id, struct s2d_hit* out_hit, struct sdis_interface** out_interf, struct sdis_interface_fragment* out_frag); extern LOCAL_SYM res_T find_next_fragment_3d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double in_pos[3], const double in_dir[3], /* Always in 3D */ const struct s3d_hit* in_hit, const double time, + const unsigned enc_id, struct s3d_hit* out_hit, struct sdis_interface** out_interf, struct sdis_interface_fragment* out_frag); 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 @@ -90,16 +90,17 @@ XD(check_handle_external_net_flux_args) static INLINE double /* [W/m^2/sr] */ XD(direct_contribution) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct source_sample* sample, const double pos[DIM], + const unsigned enc_id, /* Current enclosure */ 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->dst, hit_from, &hit); + XD(trace_ray)(scn, pos, sample->dir, sample->dst, enc_id, hit_from, &hit); if(!SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */ /* Note that the value returned is not the source's actual radiance, but the @@ -112,11 +113,12 @@ XD(direct_contribution) static res_T XD(compute_incident_diffuse_flux) - (const struct sdis_scene* scn, + (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 unsigned enc_id, /* Current enclosure */ const struct sXd(hit)* in_hit, /* Current intersection */ struct incident_diffuse_flux* diffuse_flux) /* [W/m^2] */ { @@ -158,7 +160,7 @@ XD(compute_incident_diffuse_flux) d3_minus(wi, dir); /* Always in 3D */ res = XD(find_next_fragment) - (scn, pos, dir, &hit, time, &hit, &interf, &frag); + (scn, pos, dir, &hit, time, enc_id, &hit, &interf, &frag); if(res != RES_OK) goto error; if(SXD_HIT_NONE(&hit)) { @@ -203,7 +205,8 @@ XD(compute_incident_diffuse_flux) if(res != RES_OK) goto error; if(!SOURCE_SAMPLE_NONE(&src_sample)) { - const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); + const double Ld = XD(direct_contribution) + (scn, &src_sample, pos, enc_id, &hit); L = Ld; /* [W/m^2/sr] */ } @@ -224,7 +227,8 @@ XD(compute_incident_diffuse_flux) /* The source is above the surface */ } else { - const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); + 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] */ } } @@ -244,7 +248,7 @@ error: ******************************************************************************/ res_T XD(handle_external_net_flux) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T) @@ -270,6 +274,7 @@ XD(handle_external_net_flux) struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; /* Miscellaneous */ + unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; double sum_h = 0; double emissivity = 0; /* Emissivity */ double Ld = 0; /* Incident radiance [W/m^2/sr] */ @@ -291,6 +296,9 @@ XD(handle_external_net_flux) frag.side = SDIS_BACK; } + /* Retrieve the enclosures */ + scene_get_enclosure_ids(scn, args->XD(hit)->prim.prim_id, enc_ids); + /* No external sources <=> no external fluxes. Nothing to do */ handle_flux = interface_side_is_external_flux_handled(args->interf, &frag); handle_flux = handle_flux && (scn->source != NULL); @@ -317,13 +325,14 @@ XD(handle_external_net_flux) * interface side */ cos_theta = d3_dot(N, src_sample.dir); if(cos_theta > 0) { - Ld = XD(direct_contribution)(scn, &src_sample, frag.P, args->XD(hit)); + Ld = XD(direct_contribution) + (scn, &src_sample, frag.P, enc_ids[frag.side], args->XD(hit)); incident_flux_direct = cos_theta * Ld / src_sample.pdf; /* [W/m^2] */ } /* Calculate the incident diffuse flux [W/m^2] */ - res = XD(compute_incident_diffuse_flux) - (scn, rng, frag.P, N, frag.time, args->XD(hit), &incident_flux_diffuse); + res = XD(compute_incident_diffuse_flux)(scn, rng, 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_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h @@ -175,14 +175,14 @@ HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL__; extern LOCAL_SYM res_T handle_external_net_flux_2d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T); extern LOCAL_SYM res_T handle_external_net_flux_3d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T); diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -236,7 +236,7 @@ XD(trace_radiative_path) d3_minus(wi, dir); res = XD(find_next_fragment)(scn, pos, dir, &rwalk->XD(hit), - rwalk->vtx.time, &rwalk->XD(hit), &interf, &frag); + rwalk->vtx.time, rwalk->enc_id, &rwalk->XD(hit), &interf, &frag); if(res != RES_OK) goto error; /* The path reaches the radiative environment */ @@ -348,10 +348,11 @@ error: void XD(trace_ray) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double pos[DIM], const double dir[3], const double distance, + const unsigned enc_id, const struct sXd(hit)* hit_from, struct sXd(hit)* hit) { @@ -366,7 +367,9 @@ XD(trace_ray) ray_range[0] = 0; ray_range[1] = (float)distance; filter_data.XD(hit) = *hit_from; - filter_data.epsilon = 1.e-4; + filter_data.epsilon = 1.e-6; + filter_data.scn = scn; /* Enable the filtering wrt the enclosure id */ + filter_data.enc_id = enc_id; #if DIM == 2 SXD(scene_view_trace_ray_3d (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit)); @@ -378,11 +381,12 @@ XD(trace_ray) res_T XD(find_next_fragment) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, const double in_pos[DIM], const double in_dir[3], /* Always in 3D */ const struct sXd(hit)* in_hit, const double time, + const unsigned enc_id, struct sXd(hit)* out_hit, struct sdis_interface** out_interf, struct sdis_interface_fragment* out_frag) @@ -415,7 +419,7 @@ XD(find_next_fragment) res = RES_OK; /* Find the following surface along the direction of propagation */ - XD(trace_ray)(scn, rt_pos, in_dir, INF, in_hit, &hit); + XD(trace_ray)(scn, rt_pos, in_dir, INF, enc_id, in_hit, &hit); if(SXD_HIT_NONE(&hit)) break; /* Retrieve the current position and normal */