stardis-solver

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

commit 4f7a34d95ab77a42ec48189d3197ebc7b8e5d371
parent 07aaa4ad76556d0477eb54ce99154bfd33f9b498
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 15 Dec 2021 12:37:12 +0100

Add the branch_id variable to the registered heat vertex

This index is the same for all the vertices of a line strip.  It is
equal to the index + 1 of the path from which the branch is spawned. The
main path has an index of 0.

Diffstat:
Msrc/sdis.h | 3++-
Msrc/sdis_heat_path.h | 13+++++++++++++
Msrc/sdis_heat_path_boundary_Xd_c.h | 3++-
Msrc/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h | 33++++++++++++++++++++-------------
Msrc/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h | 63++++++++++++++++++++++++---------------------------------------
Msrc/sdis_heat_path_conductive_Xd.h | 4++--
Msrc/sdis_heat_path_convective_Xd.h | 8++++----
Msrc/sdis_heat_path_radiative_Xd.h | 14++++++++++----
Msrc/sdis_misc.h | 6++++--
Msrc/sdis_realisation.c | 2+-
Msrc/sdis_realisation_Xd.h | 4++--
Msrc/test_sdis_picard.c | 26++++++++++++++++++++++++++
Msrc/test_sdis_utils.c | 15++++++++++++++-
13 files changed, 124 insertions(+), 70 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -259,8 +259,9 @@ struct sdis_heat_vertex { double time; double weight; enum sdis_heat_vertex_type type; + int branch_id; }; -#define SDIS_HEAT_VERTEX_NULL__ {{0,0,0}, 0, 0, SDIS_HEAT_VERTEX_CONDUCTION} +#define SDIS_HEAT_VERTEX_NULL__ {{0,0,0}, 0, 0, SDIS_HEAT_VERTEX_CONDUCTION, 0} static const struct sdis_heat_vertex SDIS_HEAT_VERTEX_NULL = SDIS_HEAT_VERTEX_NULL__; diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -176,6 +176,19 @@ error: goto exit; } +static INLINE void +heat_path_increment_sub_path_branch_id + (struct sdis_heat_path* path, + const size_t ivtx_begin, + const size_t ivtx_end) +{ + size_t ivtx; + FOR_EACH(ivtx, ivtx_begin, ivtx_end) { + struct sdis_heat_vertex* vtx = heat_path_get_vertex(path, ivtx); + vtx->branch_id += 1; + } +} + /* Generate the dynamic array of heat paths */ #define DARRAY_NAME heat_path #define DARRAY_DATA struct sdis_heat_path diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h @@ -855,7 +855,8 @@ XD(solid_reinjection) (args->rwalk_ctx->heat_path, &args->rwalk->vtx, args->T->value, - SDIS_HEAT_VERTEX_CONDUCTION); + SDIS_HEAT_VERTEX_CONDUCTION, + (int)args->rwalk_ctx->nbranchings); if(res != RES_OK) goto error; exit: diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h @@ -152,7 +152,7 @@ XD(solid_fluid_boundary_picard1_path) double delta; /* Orthogonal fitted reinjection dst at the boundary */ double r; - enum sdis_heat_vertex_type current_vertex_type; + struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL; enum sdis_side solid_side = SDIS_SIDE_NULL__; enum sdis_side fluid_side = SDIS_SIDE_NULL__; res_T res = RES_OK; @@ -222,17 +222,18 @@ XD(solid_fluid_boundary_picard1_path) p_conv = h_conv / h_hat; p_cond = h_cond / h_hat; - /* Fetch the current heat vertex type. This will be used below to restart the - * registration of the heat path geometry after a null collision */ - if(ctx->heat_path) { - current_vertex_type = heat_path_get_last_vertex(ctx->heat_path)->type; - } + /* Fetch the last registered heat path vertex */ + if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path); /* Null collision */ for(;;) { double h_radi; /* Radiative coefficient */ double p_radi; /* Radiative proba */ + /* Indices of the registered vertex of the sampled radiative path */ + size_t ihvtx_radi_begin; + size_t ihvtx_radi_end; + r = ssp_rng_canonical(rng); /* Switch in convective path */ @@ -263,6 +264,12 @@ XD(solid_fluid_boundary_picard1_path) /* From there, we know the path is either a radiative path or a * null-collision */ + if(ctx->heat_path) { + /* Fetch the index of the first vertex of the radiative path that is + * going to be traced i.e. the last registered vertex */ + ihvtx_radi_begin = heat_path_get_vertices_count(ctx->heat_path) - 1; + } + /* Sample a radiative path and get the Tref at its end. */ T_s = *T; rwalk_s = *rwalk; @@ -297,14 +304,14 @@ XD(solid_fluid_boundary_picard1_path) } if(ctx->heat_path) { - /* Add a break into the heat path geometry and restart it from the - * current position. The sampled radiative path becomes a branch of the - * current sampled path */ - res = heat_path_add_break(ctx->heat_path); - if(res != RES_OK) goto error; + /* Set the sampled radiative path as a branch of the current path */ + ihvtx_radi_end = heat_path_get_vertices_count(ctx->heat_path); + heat_path_increment_sub_path_branch_id + (ctx->heat_path, ihvtx_radi_begin, ihvtx_radi_end); - res = register_heat_vertex - (ctx->heat_path, &rwalk->vtx, T->value, current_vertex_type); + /* Add a break into the heat path geometry and restart it from the + * position of the input random walk. */ + res = heat_path_restart(ctx->heat_path, &hvtx); if(res != RES_OK) goto error; } } diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h @@ -26,42 +26,6 @@ #include "sdis_Xd_begin.h" /******************************************************************************* - * Non generic helper functions - ******************************************************************************/ -#ifndef SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H -#define SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H - -static INLINE res_T -restart_heat_path - (struct sdis_heat_path* path, - const struct sdis_heat_vertex* vtx) -{ - size_t nverts = 0; - size_t nbreaks = 0; - res_T res = RES_OK; - - if(!path) goto exit; - ASSERT(vtx); - - nbreaks = darray_size_t_size_get(&path->breaks); - nverts = darray_heat_vertex_size_get(&path->vertices); - - res = heat_path_add_break(path); - if(res != RES_OK) goto error; - res = heat_path_add_vertex(path, vtx); - if(res != RES_OK) goto error; - -exit: - return res; -error: - CHK(darray_size_t_resize(&path->breaks, nbreaks) == RES_OK); - CHK(darray_heat_vertex_resize(&path->vertices, nverts) == RES_OK); - goto exit; -} - -#endif /* SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H */ - -/******************************************************************************* * Generic helper functions ******************************************************************************/ static INLINE res_T @@ -96,8 +60,9 @@ XD(sample_path) heat_vtx.time = rwalk.vtx.time; heat_vtx.weight = 0; heat_vtx.type = SDIS_HEAT_VERTEX_RADIATIVE; + heat_vtx.branch_id = (int)ctx->nbranchings + 1; - res = restart_heat_path(ctx->heat_path, &heat_vtx); + res = heat_path_restart(ctx->heat_path, &heat_vtx); if(res != RES_OK) goto error; } @@ -248,6 +213,10 @@ XD(solid_fluid_boundary_picardN_path) double p_radi, p_radi_min, p_radi_max; /* Radiative probas */ double T0, T1, T2, T3, T4, T5; /* Computed temperatures */ + /* Indices of the registered vertex of the sampled radiative path */ + size_t ihvtx_radi_begin; + size_t ihvtx_radi_end; + r = ssp_rng_canonical(rng); /* Switch in convective path */ @@ -275,6 +244,12 @@ XD(solid_fluid_boundary_picardN_path) break; } + if(ctx->heat_path) { + /* Fetch the index of the first vertex of the radiative path that is + * going to be traced i.e. the last registered vertex */ + ihvtx_radi_begin = heat_path_get_vertices_count(ctx->heat_path) - 1; + } + /* Sample a radiative path */ T_s = *T; rwalk_s = *rwalk; @@ -283,6 +258,12 @@ XD(solid_fluid_boundary_picardN_path) res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s); if(res != RES_OK) goto error; + if(ctx->heat_path) { + /* Fetch the index after the last registered vertex of the sampled + * radiative path */ + ihvtx_radi_end = heat_path_get_vertices_count(ctx->heat_path); + } + /* Fetch the last registered heat path vertex of the radiative path */ if(ctx->heat_path) hvtx_s = *heat_path_get_last_vertex(ctx->heat_path); @@ -295,13 +276,17 @@ XD(solid_fluid_boundary_picardN_path) /* Define some helper macros */ #define SWITCH_IN_RADIATIVE { \ *rwalk = rwalk_s; *T = T_s; \ - res = restart_heat_path(ctx->heat_path, &hvtx_s); \ + res = heat_path_restart(ctx->heat_path, &hvtx_s); \ if(res != RES_OK) goto error; \ } (void)0 #define NULL_COLLISION { \ - res = restart_heat_path(ctx->heat_path, &hvtx); \ + res = heat_path_restart(ctx->heat_path, &hvtx); \ if(res != RES_OK) goto error; \ + if(ctx->heat_path) { \ + heat_path_increment_sub_path_branch_id \ + (ctx->heat_path, ihvtx_radi_begin, ihvtx_radi_end); \ + } \ } (void)0 #define COMPUTE_TEMPERATURE(Result, RWalk, Temp) { \ diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -468,8 +468,8 @@ XD(conductive_path) XD(move_pos)(rwalk->vtx.P, dir0, delta); /* Register the new vertex against the heat path */ - res = register_heat_vertex - (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, + SDIS_HEAT_VERTEX_CONDUCTION, (int)ctx->nbranchings); if(res != RES_OK) goto error; ++istep; diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -60,8 +60,8 @@ XD(register_heat_vertex_in_fluid) fX(add)(pos, org, fX(mulf)(dir, dir, dst)); dX_set_fX(vtx.P, pos); - return register_heat_vertex - (ctx->heat_path, &vtx, weight, SDIS_HEAT_VERTEX_CONVECTION); + return register_heat_vertex(ctx->heat_path, &vtx, weight, + SDIS_HEAT_VERTEX_CONVECTION, (int)ctx->nbranchings); } /******************************************************************************* @@ -283,8 +283,8 @@ XD(convective_path) } /* Register the new vertex against the heat path */ - res = register_heat_vertex - (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONVECTION); + res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, + SDIS_HEAT_VERTEX_CONVECTION, (int)ctx->nbranchings); if(res != RES_OK) goto error; /* Setup the fragment of the sampled position into the enclosure. */ diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -42,12 +42,16 @@ XD(trace_radiative_path) * are assumed to be extruded to the infinity along the Z dimension. */ float N[3] = {0, 0, 0}; float dir[3] = {0, 0, 0}; + int branch_id; res_T res = RES_OK; ASSERT(scn && ray_dir && ctx && rwalk && rng && T); f3_set(dir, ray_dir); + /* (int)ctx->nbranchings < 0 <=> Beginning of the realisation */ + branch_id = MMAX((int)ctx->nbranchings, 0); + /* Launch the radiative random walk */ for(;;) { const struct sdis_interface* interf = NULL; @@ -86,12 +90,14 @@ XD(trace_radiative_path) if(ctx->heat_path) { const float empirical_dst = 0.1f; struct sdis_rwalk_vertex vtx; + + vtx = rwalk->vtx; vtx.P[0] += dir[0] * empirical_dst; vtx.P[1] += dir[1] * empirical_dst; vtx.P[2] += dir[2] * empirical_dst; - res = register_heat_vertex - (ctx->heat_path, &vtx, T->value, SDIS_HEAT_VERTEX_RADIATIVE); + res = register_heat_vertex(ctx->heat_path, &vtx, T->value, + SDIS_HEAT_VERTEX_RADIATIVE, branch_id); if(res != RES_OK) goto error; } break; @@ -119,8 +125,8 @@ XD(trace_radiative_path) XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance); /* Register the random walk vertex against the heat path */ - res = register_heat_vertex - (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_RADIATIVE); + res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, + SDIS_HEAT_VERTEX_RADIATIVE, branch_id); if(res != RES_OK) goto error; /* Fetch the new interface and setup the hit fragment */ diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -127,10 +127,11 @@ register_heat_vertex (struct sdis_heat_path* path, const struct sdis_rwalk_vertex* vtx, const double weight, - const enum sdis_heat_vertex_type type) + const enum sdis_heat_vertex_type type, + const int branch_id) { struct sdis_heat_vertex heat_vtx = SDIS_HEAT_VERTEX_NULL; - ASSERT(vtx); + ASSERT(vtx && branch_id >= 0); if(!path) return RES_OK; @@ -140,6 +141,7 @@ register_heat_vertex heat_vtx.time = vtx->time; heat_vtx.weight = weight; heat_vtx.type = type; + heat_vtx.branch_id = branch_id; return heat_path_add_vertex(path, &heat_vtx); } diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c @@ -71,7 +71,7 @@ ray_realisation_3d /* Register the starting position against the heat path */ res = register_heat_vertex - (args->heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE); + (args->heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE, 0); if(res != RES_OK) goto error; res = trace_radiative_path_3d(scn, dir, &ctx, &rwalk, args->rng, &T); diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -194,7 +194,7 @@ XD(probe_realisation) type = args->medium->type == SDIS_SOLID ? SDIS_HEAT_VERTEX_CONDUCTION : SDIS_HEAT_VERTEX_CONVECTION; - res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0, type); + res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0, type, 0); if(res != RES_OK) goto error; if(t0 >= rwalk.vtx.time) { @@ -289,7 +289,7 @@ XD(boundary_realisation) #endif res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0/*weight*/, - SDIS_HEAT_VERTEX_CONDUCTION); + SDIS_HEAT_VERTEX_CONDUCTION, 0/*Branch id*/); if(res != RES_OK) goto error; ctx.green_path = args->green_path; diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c @@ -460,6 +460,26 @@ test_picard } static void +register_heat_paths(struct sdis_scene* scn, const size_t picard_order, FILE* stream) +{ + struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_estimator* estimator = NULL; + CHK(scn && picard_order >= 1 && stream); + + + probe_args.nrealisations = 10; + probe_args.position[0] = 0.05; + probe_args.position[1] = 0; + probe_args.position[2] = 0; + probe_args.picard_order = picard_order; + probe_args.register_paths = SDIS_HEAT_PATH_ALL; + printf("Register %lu heat paths.\n", probe_args.nrealisations); + OK(sdis_solve_probe(scn, &probe_args, &estimator)); + dump_heat_paths(stream, estimator); + OK(sdis_estimator_ref_put(estimator)); +} + +static void create_scene_3d (struct sdis_device* dev, struct sdis_interface* interfaces[INTERFACES_COUNT__], @@ -545,6 +565,7 @@ create_scene_2d int main(int argc, char** argv) { + FILE* stream = NULL; struct mem_allocator allocator; struct sdis_device* dev = NULL; @@ -622,6 +643,8 @@ main(int argc, char** argv) create_scene_2d(dev, interfaces, &scn_2d); create_scene_3d(dev, interfaces, &scn_3d); + CHK((stream = tmpfile()) != NULL); + /* Test picard1 with a constant Tref <=> regular linearisation */ printf("Test Picard1 with a constant Tref of 300 K\n"); ref.T = 314.99999999999989; @@ -691,6 +714,8 @@ main(int argc, char** argv) OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp)); test_picard(scn_2d, 3/*Picard order*/, &ref); test_picard(scn_3d, 3/*Picard order*/, &ref); + register_heat_paths(scn_2d, 3/*Picard order*/, stream); + register_heat_paths(scn_3d, 3/*Picard order*/, stream); printf("\n"); t_range[0] = 280; @@ -746,6 +771,7 @@ main(int argc, char** argv) OK(sdis_medium_ref_put(solid)); OK(sdis_medium_ref_put(dummy)); OK(sdis_device_ref_put(dev)); + CHK(fclose(stream) == 0); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c @@ -17,6 +17,7 @@ #include <rsys/math.h> enum heat_vertex_attrib { + HEAT_VERTEX_BRANCH_ID, HEAT_VERTEX_WEIGHT, HEAT_VERTEX_TIME, HEAT_VERTEX_TYPE @@ -250,11 +251,15 @@ dump_heat_path_vertex_attribs struct sdis_heat_vertex vtx; OK(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); switch(attr) { + case HEAT_VERTEX_BRANCH_ID: + fprintf(stream, "%i\n", vtx.branch_id); + break; case HEAT_VERTEX_WEIGHT: fprintf(stream, "%g\n", vtx.weight); break; case HEAT_VERTEX_TIME: - fprintf(stream, "%g\n", IS_INF(vtx.time) ? FLT_MAX : vtx.time); + fprintf(stream, "%g\n", + IS_INF(vtx.time) || vtx.time > FLT_MAX ? -1 : vtx.time); break; case HEAT_VERTEX_TYPE: switch(vtx.type) { @@ -413,6 +418,14 @@ dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator) dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_TIME); } + /* Write the branch id of the random walk vertices */ + fprintf(stream, "SCALARS BranchID int 1\n"); + fprintf(stream, "LOOKUP_TABLE default\n"); + FOR_EACH(ipath, 0, npaths) { + OK(sdis_estimator_get_path(estimator, ipath, &path)); + dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_BRANCH_ID); + } + /* Write path type */ fprintf(stream, "CELL_DATA %lu\n", (unsigned long)nstrips_overall); fprintf(stream, "SCALARS Path_Type float 1\n");