stardis-solver

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

commit 85cad347d896dc98aed22abb9c2a96a8536fb807
parent 9d0c19a29254dac9e0d3adef7d607b2f8693a97e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 23 May 2018 12:52:05 +0200

Output more informations on error during a diffusive random walk

Diffstat:
Msrc/sdis_scene.c | 25++++++++++++++++++++-----
Msrc/sdis_scene_c.h | 15+++++++++++++++
Msrc/sdis_solve.c | 4++--
Msrc/sdis_solve_Xd.h | 51+++++++++++++++++++++++++++++++++++++++------------
4 files changed, 76 insertions(+), 19 deletions(-)

diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -24,9 +24,6 @@ #include <rsys/double3.h> #include <rsys/mem_allocator.h> -#include <star/s2d.h> -#include <star/s3d.h> - #include <limits.h> /* Context used to wrap the user geometry to Star-XD. */ @@ -455,6 +452,7 @@ static INLINE res_T scene_get_medium_2d (const struct sdis_scene* scn, const double pos[2], + struct get_medium_info* info, /* May be NULL */ const struct sdis_medium** out_medium) { const struct sdis_medium* medium = NULL; @@ -502,6 +500,14 @@ scene_get_medium_2d interf = scene_get_interface(scn, hit.prim.prim_id); medium = interface_get_medium (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK); + + /* Register the get_medium_info */ + if(info) { + f2_set(info->pos_tgt, attr.value); + f2_set(info->ray_org, P); + f2_set(info->ray_dir, dir); + info->hit_2d = hit; + } break; } } @@ -519,6 +525,7 @@ static INLINE res_T scene_get_medium_3d (const struct sdis_scene* scn, const double pos[3], + struct get_medium_info* info, const struct sdis_medium** out_medium) { const struct sdis_medium* medium = NULL; @@ -566,6 +573,13 @@ scene_get_medium_3d interf = scene_get_interface(scn, hit.prim.prim_id); medium = interface_get_medium (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK); + + if(info) { + f3_set(info->pos_tgt, attr.value); + f3_set(info->ray_org, P); + f3_set(info->ray_dir, dir); + info->hit_3d = hit; + } break; } } @@ -775,10 +789,11 @@ res_T scene_get_medium (const struct sdis_scene* scn, const double pos[], + struct get_medium_info* info, const struct sdis_medium** out_medium) { return scene_is_2d(scn) - ? scene_get_medium_2d(scn, pos, out_medium) - : scene_get_medium_3d(scn, pos, out_medium); + ? scene_get_medium_2d(scn, pos, info, out_medium) + : scene_get_medium_3d(scn, pos, info, out_medium); } diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -19,6 +19,20 @@ #include <rsys/dynamic_array.h> #include <rsys/ref_count.h> +#include <star/s2d.h> +#include <star/s3d.h> + +struct get_medium_info { + /* Targeted position */ + float pos_tgt[3]; + /* Ray trace to the targeted position in order to define the current medium */ + float ray_org[3]; + float ray_dir[3]; + /* Hit encouters along the ray and used to define the current medium */ + struct s2d_hit hit_2d; + struct s3d_hit hit_3d; +}; + static INLINE void interface_init (struct mem_allocator* allocator, @@ -62,6 +76,7 @@ extern LOCAL_SYM res_T scene_get_medium (const struct sdis_scene* scene, const double position[], + struct get_medium_info* info, /* May be NULL */ const struct sdis_medium** medium); static FINLINE int diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -209,7 +209,7 @@ sdis_solve_probe if(res != RES_OK) goto error; /* Retrieve the medium in which the submitted position lies */ - res = scene_get_medium(scn, position, &medium); + res = scene_get_medium(scn, position, NULL, &medium); if(res != RES_OK) goto error; /* Here we go! Launch the Monte Carlo estimation */ @@ -448,7 +448,7 @@ sdis_solve_camera } /* Retrieve the medium in which the submitted position lies */ - res = scene_get_medium(scn, cam->position, &medium); + res = scene_get_medium(scn, cam->position, NULL, &medium); if(res != RES_OK) goto error; if(medium->type != SDIS_FLUID) { diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -28,7 +28,7 @@ /* Emperical scale factor to apply to the upper bound of the ray range in order * to handle numerical imprecisions */ -#define RAY_RANGE_MAX_SCALE 1.0001f +#define RAY_RANGE_MAX_SCALE 1.01f #define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */ @@ -663,7 +663,7 @@ XD(solid_temperature) (void)ctx; /* Check the random walk consistency */ - CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK); + CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK); if(mdm != rwalk->mdm) { log_err(scn->dev, "%s: invalid solid random walk. " "Unexpected medium at {%g, %g, %g}.\n", @@ -674,6 +674,7 @@ XD(solid_temperature) dX(set)(position_start, rwalk->vtx.P); do { /* Solid random walk */ + struct get_medium_info info; struct sXd(hit) hit0, hit1; double lambda; /* Thermal conductivity */ double rho; /* Volumic mass */ @@ -760,7 +761,7 @@ XD(solid_temperature) /* Fetch the current medium */ if(SXD_HIT_NONE(&rwalk->hit)) { - CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK); + CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK); } else { const struct sdis_interface* interf; interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); @@ -771,15 +772,28 @@ XD(solid_temperature) if(mdm != rwalk->mdm) { log_err(scn->dev, "%s: inconsistent medium during the solid random walk.\n", FUNC_NAME); - if(DIM == 2) { - log_err(scn->dev, - " start position: %g %g; current position: %g %g\n", - SPLIT2(position_start), SPLIT2(rwalk->vtx.P)); - } else { - log_err(scn->dev, - " start position: %g %g %g; current position: %g %g %g\n", - SPLIT3(position_start), SPLIT3(rwalk->vtx.P)); +#if DIM == 2 + #define VEC_STR "%g %g" + #define VEC_SPLIT SPLIT2 +#else + #define VEC_STR "%g %g %g" + #define VEC_SPLIT SPLIT3 +#endif + log_err(scn->dev, + " start position: " VEC_STR "; current position: " VEC_STR "\n", + VEC_SPLIT(position_start), VEC_SPLIT(rwalk->vtx.P)); + if(SXD_HIT_NONE(&rwalk->hit)) { + float hit_pos[DIM]; + fX(mulf)(hit_pos, info.ray_dir, info.XD(hit).distance); + fX(add)(hit_pos, info.ray_org, hit_pos); + log_err(scn->dev, " ray org: " VEC_STR "; ray dir: " VEC_STR "\n", + VEC_SPLIT(info.ray_org), VEC_SPLIT(info.ray_dir)); + log_err(scn->dev, " targeted point: " VEC_STR "\n", + VEC_SPLIT(info.pos_tgt)); + log_err(scn->dev, " hit pos: " VEC_STR "\n", VEC_SPLIT(hit_pos)); } +#undef VEC_STR +#undef VEC_SPLIT return RES_BAD_OP; } @@ -804,15 +818,28 @@ XD(compute_temperature) struct XD(temperature)* stack = NULL; size_t istack = 0; #endif + const size_t max_fails = 10; res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); do { + /* Save the current random walk state */ + const struct XD(rwalk) rwalk_bkp = *rwalk; + const struct XD(temperature) T_bkp = *T; + + size_t nfails = 0; + #ifndef NDEBUG sa_push(stack, *T); ++istack; #endif - res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T); + + /* Reject the current step if a BAD_OP occurs and retry up to "max_fails" + * times */ + do { + res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T); + if(res == RES_BAD_OP) { *rwalk = rwalk_bkp; *T = T_bkp; } + } while(res == RES_BAD_OP && ++nfails < max_fails); if(res != RES_OK) goto error; } while(!T->done);