stardis-solver

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

commit edc3c424bdffc248fc2adbe7a794fdc84d7225c3
parent 743a809ffd4b2b60b877131342ce7930f6939fa8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 22 Dec 2023 12:17:09 +0100

Updating solver management of external sources

The interface functors used to manage external sources have been
removed. They are replaced by the optional source submitted as an input
argument to scene creation. The Monte Carlo algorithm calculating the
external flux of an interface is updated accordingly. Note that
interfaces now have a variable to enable or disable external flux
calculation per interface, regardless of the presence of an external
source.

In contrast to the previous proposal, the Monte Carlo algorithm handles
only one external source. Although it may seem less generic, it is
intended to be more realistic and simpler. Dealing with multiple sources
would require the implementation of an importance sampling strategy. The
current proposal does not (falsely) suggest that the solver could handle
several external sources in a single step. Further development would be
required.

Nevertheless, it should be noted that the algorithm remains independent
of the type of source as long as we are able to sample it and test its
intersection by a ray.

Diffstat:
Msrc/sdis.h | 50++++++++++----------------------------------------
Msrc/sdis_heat_path_boundary_Xd_handle_external_net_flux.h | 82++++++++++++++++++++++++++++++++++++++-----------------------------------------
Msrc/sdis_interface_c.h | 21+++------------------
Msrc/sdis_scene.c | 1+
Msrc/sdis_scene_Xd.h | 5+++++
Msrc/sdis_scene_c.h | 2++
Msrc/test_sdis_utils.h | 3+--
7 files changed, 61 insertions(+), 103 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -201,37 +201,6 @@ typedef double (const struct sdis_interface_fragment* frag, /* Interface position */ struct sdis_data* data); /* User data */ -/* A sample of an external source */ -struct sdis_external_sources_sample { - double dir[3]; /* Direction _to_ the source */ - double pdf; /* Pdf of sampled direction */ - double distance; /* Distance to the source [m] */ - double radiance; /* Constant source radiance [W/m^2/sr] */ -}; -#define SDIS_EXTERNAL_SOURCES_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0} -static const struct sdis_external_sources_sample -SDIS_EXTERNAL_SOURCES_SAMPLE_NULL = SDIS_EXTERNAL_SOURCES_SAMPLE_NULL__; - -/* Is the sample valid */ -#define SDIS_EXTERNAL_SOURCES_SAMPLE_NONE(Sample) ((Sample)->pdf != 0) - -/* Functor for sampling external sources. The returned sample is used to - * evaluate an external flux at the interface. */ -typedef res_T -(*sdis_sample_external_sources_T) - (const struct sdis_interface_fragment* frag, /* Interface position */ - struct ssp_rng* rng, /* Random Number Generator to use */ - struct sdis_external_sources_sample* sample, /* Returned sample */ - struct sdis_data* data); /* User data */ - -/* Functor returning the external sources hit by a ray */ -typedef res_T -(*sdis_trace_external_sources_T) - (const struct sdis_interface_fragment* frag, /* Interface position */ - const double dir[3], /* Ray direction */ - struct sdis_external_sources_sample* sample, /* Hit source */ - struct sdis_data* data); /* User data */ - /* Define the physical properties of a solid */ struct sdis_solid_shader { /* Properties */ @@ -291,15 +260,11 @@ struct sdis_interface_side_shader { /* Reference temperature used in Picard 1 */ sdis_interface_getter_T reference_temperature; - /* Manage external sources, i.e. sample them to evaluate the corresponding - * external flux. Can be NULL <=> no external source */ - sdis_sample_external_sources_T sample_external_sources; - - /* Does the ray target an external source? Can only be NULL if - * sample_external_sources is also NULL */ - sdis_trace_external_sources_T trace_external_sources; + /* Define whether external sources interact with the interface, i.e. whether + * external fluxes should be processed or not */ + int handle_external_flux; }; -#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL, NULL, NULL, NULL } +#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL, NULL, 1 } static const struct sdis_interface_side_shader SDIS_INTERFACE_SIDE_SHADER_NULL = SDIS_INTERFACE_SIDE_SHADER_NULL__; @@ -474,6 +439,10 @@ struct sdis_scene_create_args { /* Min/max temperature used to linearise the radiative temperature */ double t_range[2]; + + /* External source. Can be NULL <=> no external flux will be calculated on + * scene interfaces */ + struct sdis_source* source; }; #define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \ @@ -485,7 +454,8 @@ struct sdis_scene_create_args { 0, /* #vertices */ \ 1.0, /* #Floating point to meter scale factor */ \ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\ - {0.0, -1.0} /* Temperature range */ \ + {0.0, -1.0}, /* Temperature range */ \ + NULL /* source */ \ } static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT = SDIS_SCENE_CREATE_ARGS_DEFAULT__; 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 @@ -17,6 +17,7 @@ #include "sdis_interface_c.h" #include "sdis_log.h" #include "sdis_scene_c.h" +#include "sdis_source_c.h" #include <rsys/cstr.h> /* res_to_cstr */ @@ -103,26 +104,26 @@ sample_brdf ******************************************************************************/ static INLINE res_T XD(check_handle_external_net_flux_args) - (const struct sdis_device* dev, + (const struct sdis_scene* scn, const char* func_name, const struct XD(handle_external_net_flux_args)* args) { - sdis_sample_external_sources_T functor = NULL; + int net_flux = 0; res_T res = RES_OK; /* Handle bugs */ - ASSERT(dev && func_name && args); + ASSERT(scn && 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); + net_flux = interface_side_is_external_flux_handled(args->interf, args->frag); + net_flux = net_flux && (scn->source != NULL); - if(functor && args->picard_order != 0) { + if(net_flux && args->picard_order != 0) { res = RES_BAD_ARG; - log_err(dev, + log_err(scn->dev, "%s: Impossible to process external fluxes when Picard order is not " "equal to 1; Picard order is currently set to %lu.\n", func_name, (unsigned long)args->picard_order); @@ -164,7 +165,7 @@ XD(trace_ray) static INLINE double /* [W/m^2/sr] */ XD(direct_contribution) (const struct sdis_scene* scn, - struct sdis_external_sources_sample* sample, + struct source_sample* sample, const double pos[DIM], const struct sXd(hit)* hit_from) { @@ -172,7 +173,7 @@ XD(direct_contribution) ASSERT(scn && sample && pos && hit_from); /* Is the source hidden */ - XD(trace_ray)(scn, pos, sample->dir, sample->distance, hit_from, &hit); + XD(trace_ray)(scn, pos, sample->dir, sample->dst, hit_from, &hit); if(SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */ return sample->radiance; /* [W/m^2/sr] */ @@ -239,7 +240,7 @@ XD(compute_incident_diffuse_flux) const struct sXd(hit)* in_hit) /* Current intersection */ { struct sXd(hit) hit = SXD_HIT_NULL; - double pos[DIM] = {0}; + double pos[3] = {0}; /* In 3D for ray tracing ray to the source */ 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] */ @@ -256,10 +257,7 @@ XD(compute_incident_diffuse_flux) 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; + struct source_sample src_sample = SOURCE_SAMPLE_NULL; /* Interface */ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; @@ -289,10 +287,6 @@ XD(compute_incident_diffuse_flux) 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; @@ -304,32 +298,33 @@ XD(compute_incident_diffuse_flux) /* 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 */ + res = source_trace_to(scn->source, pos, brdf_sample.dir, time, &src_sample); + CHK(res == RES_OK); - if(!SDIS_EXTERNAL_SOURCES_SAMPLE_NONE(&extsrc_sample)) { - const double Ld = XD(direct_contribution)(scn, &extsrc_sample, pos, &hit); + if(!SOURCE_SAMPLE_NONE(&src_sample)) { + const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); L = Ld; /* [W/m^2] */ } /* Calculate the direct contribution of the rebound is diffuse */ } else { + double cos_theta = 0; 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 */ + res = source_sample(scn->source, rng, pos, time, &src_sample); + CHK(res == RES_OK); + cos_theta = d3_dot(src_sample.dir, N); /* The source is behind the surface */ - if(d3_dot(extsrc_sample.dir, N) <= 0) { + if(cos_theta <= 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] */ + const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); + L = Ld * cos_theta/PI * src_sample.pdf; /* [W/m^2] */ } } incident_diffuse_flux += L; @@ -350,9 +345,7 @@ XD(handle_external_net_flux) struct XD(temperature)* T) { /* Sampling external sources */ - struct sdis_external_sources_sample extsrc_sample = - SDIS_EXTERNAL_SOURCES_SAMPLE_NULL; - sdis_sample_external_sources_T sample_sources = NULL; + struct source_sample src_sample = SOURCE_SAMPLE_NULL; /* External flux */ double incident_flux = 0; /* [W/m^2] */ @@ -367,31 +360,34 @@ XD(handle_external_net_flux) double emissivity = 0; /* Emissivity */ double Ld = 0; /* Incident radiance [W/m^2/sr] */ double cos_theta = 0; + int handle_flux = 0; res_T res = RES_OK; ASSERT(scn && args && T); - res = XD(check_handle_external_net_flux_args)(scn->dev, FUNC_NAME, args); + res = XD(check_handle_external_net_flux_args)(scn, FUNC_NAME, args); if(res != RES_OK) goto error; - /* Retrieve the functor to sample external sources */ - sample_sources = interface_side_get_external_sources_sampling_functor - (args->interf, args->frag); - /* No external sources <=> no external fluxes. Nothing to do */ - if(!sample_sources) goto exit; + handle_flux = interface_side_is_external_flux_handled(args->interf, args->frag); + handle_flux = net_flux && (scn->source != NULL); + if(handle_flux) goto exit; - /* Sample an external sources */ - res = sample_sources(args->frag, rng, &extsrc_sample, args->interf->data); + /* Sample the external source */ + res = source_sample + (scn->source, rng, args->frag->P, args->frag->time, &src_sample); if(res != RES_OK) goto error; /* Local path data */ dX(set)(N, args->frag->Ng); if(args->frag->side == SDIS_BACK) dX(minus)(N, N); - /* 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 direct flux if the external source is above the + * interface side */ + cos_theta = d3_dot(N, src_sample.dir); + if(cos_theta > 0) { + Ld = XD(direct_contribution)(scn, &src_sample, args->frag->P, args->hit); + incident_flux_direct = cos_theta * Ld / src_sample.pdf; /* [W/m^2] */ + } /* Calculate the incident diffuse flux [W/m^2] */ incident_flux_diffuse = XD(compute_incident_diffuse_flux) diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -185,8 +185,8 @@ interface_side_get_reference_temperature ? shader->reference_temperature(frag, interf->data) : -1; } -static INLINE sdis_sample_external_sources_T -interface_side_get_external_sources_sampling_functor +static INLINE int +interface_side_is_external_flux_handled (const struct sdis_interface* interf, const struct sdis_interface_fragment* frag) { @@ -197,22 +197,7 @@ interface_side_get_external_sources_sampling_functor case SDIS_BACK: shader = &interf->shader.back; break; default: FATAL("Unreachable code\n"); break; } - return shader->sample_external_sources; -} - -static INLINE sdis_trace_external_sources_T -interface_side_get_external_sources_tracing_functor - (const struct sdis_interface* interf, - const struct sdis_interface_fragment* frag) -{ - const struct sdis_interface_side_shader* shader; - ASSERT(interf && frag); - switch(frag->side) { - case SDIS_BACK: shader = &interf->shader.back; break; - case SDIS_FRONT: shader = &interf->shader.front; break; - default: FATAL("Unreachable code\n"); break; - } - return shader->trace_external_sources; + return shader->handle_external_flux; } /******************************************************************************* diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -108,6 +108,7 @@ scene_release(ref_T * ref) if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view)); if(scn->senc2d_scn) SENC2D(scene_ref_put(scn->senc2d_scn)); if(scn->senc3d_scn) SENC3D(scene_ref_put(scn->senc3d_scn)); + if(scn->source) SDIS(source_ref_put(scn->source)); MEM_RM(dev->allocator, scn); SDIS(device_ref_put(dev)); } diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -929,6 +929,11 @@ XD(scene_create) htable_enclosure_init(dev->allocator, &scn->enclosures); htable_d_init(dev->allocator, &scn->tmp_hc_ub); + if(args->source) { + SDIS(source_ref_get(args->source)); + scn->source = args->source; + } + res = XD(run_analyze) (scn, args->nprimitives, diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -218,6 +218,8 @@ struct sdis_scene { double tmin; /* Minimum temperature of the system (In Kelvin) */ double tmax; /* Maximum temperature of the system (In Kelvin) */ + struct sdis_source* source; /* External source. May be NULL */ + ref_T ref; struct sdis_device* dev; }; diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -193,8 +193,7 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = { dummy_interface_getter, /* Emissivity */ \ dummy_interface_getter, /* Specular fraction */ \ dummy_interface_getter, /* Reference temperature */ \ - NULL, /* Sample external sources */ \ - NULL, /* Trace externa sources */ \ + 1 /* Handle external flux */ \ } static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = { dummy_interface_getter, /* Convection coef */