stardis-solver

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

commit 5570c9f9a9cc0d42eca81ceaf5899db2e58ac918
parent 37063248b849cd84003f899b77875385c981830b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 24 Jun 2020 18:44:56 +0200

Update the API of sdis_solve_probe_boundary[_green_function]

Diffstat:
Msrc/sdis.h | 46++++++++++++++++++++++++++++++----------------
Msrc/sdis_realisation_Xd.h | 2+-
Msrc/sdis_solve.c | 30++++++------------------------
Msrc/sdis_solve_probe_boundary_Xd.h | 61+++++++++++++++++++++++++++++++++----------------------------
Msrc/test_sdis_solve_boundary.c | 104++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
5 files changed, 131 insertions(+), 112 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -297,6 +297,34 @@ struct sdis_solve_probe_args { static const struct sdis_solve_probe_args SDIS_SOLVE_PROBE_ARGS_DEFAULT = SDIS_SOLVE_PROBE_ARGS_DEFAULT__; +struct sdis_solve_probe_boundary_args { + size_t nrealisations; /* #realisations */ + size_t iprim; /* Identifier of the primitive on which the probe lies */ + double uv[2]; /* Parametric coordinates of the probe onto the primitve */ + double time_range[2]; /* Observation time */ + enum sdis_side side; /* Side of iprim on which the probe lies */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + 0, /* Primitive identifier */ \ + {0,0}, /* UV */ \ + {DBL_MAX,DBL_MAX}, \ + SDIS_SIDE_NULL__, \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + SDIS_HEAT_PATH_NONE, \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_probe_boundary_args +SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT = +SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__; + /* Functor used to process the paths registered against the green function */ typedef res_T (*sdis_process_green_path_T) @@ -874,15 +902,7 @@ sdis_solve_probe SDIS_API res_T sdis_solve_probe_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_estimator** estimator); SDIS_API res_T @@ -979,13 +999,7 @@ sdis_solve_probe_green_function SDIS_API res_T sdis_solve_probe_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** green); SDIS_API res_T diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -228,7 +228,7 @@ XD(boundary_realisation) float st[2]; #endif res_T res = RES_OK; - ASSERT(uv && fp_to_meter > 0 && weight && Tref >= 0 && time >= 0); + ASSERT(uv && fp_to_meter > 0 && weight && time >= 0); T.func = XD(boundary_path); rwalk.hit_side = side; diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -248,46 +248,28 @@ sdis_solve_probe_green_function res_T sdis_solve_probe_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_2d(scn, nrealisations, iprim, uv, time_range, - side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_boundary_2d(scn, args, NULL, out_estimator); } else { - return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, time_range, - side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_boundary_3d(scn, args, NULL, out_estimator); } } res_T sdis_solve_probe_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** green) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_2d(scn, nrealisations, iprim, uv, NULL, - side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_probe_boundary_2d(scn, args, green, NULL); } else { - return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, NULL, - side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_probe_boundary_3d(scn, args, green, NULL); } } diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -33,15 +33,7 @@ static res_T XD(solve_probe_boundary) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** out_green, struct sdis_estimator** out_estimator) { @@ -52,14 +44,16 @@ XD(solve_probe_boundary) struct ssp_rng** rngs = NULL; struct accum* acc_temps = NULL; struct accum* acc_times = NULL; + size_t nrealisations = 0; int64_t irealisation = 0; size_t i; int progress = 0; + int register_paths = SDIS_HEAT_PATH_NONE; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv - || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || args->fp_to_meter <= 0 || ((unsigned)args->side >= SDIS_SIDE_NULL__)) { res = RES_BAD_ARG; goto error; } @@ -68,8 +62,10 @@ XD(solve_probe_boundary) goto error; } if(out_estimator) { - if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + if(args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } @@ -82,12 +78,12 @@ XD(solve_probe_boundary) #endif /* Check the primitive identifier */ - if(iprim >= scene_get_primitives_count(scn)) { + if(args->iprim >= scene_get_primitives_count(scn)) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. " "It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)iprim, + (unsigned long)args->iprim, (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; @@ -96,25 +92,25 @@ XD(solve_probe_boundary) /* Check parametric coordinates */ #if SDIS_XD_DIMENSION == 2 { - const double v = CLAMP(1.0 - uv[0], 0, 1); - if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) { + const double v = CLAMP(1.0 - args->uv[0], 0, 1); + if(args->uv[0] < 0 || args->uv[0] > 1 || !eq_eps(args->uv[0]+v, 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates %g." "u + (1-u) must be equal to 1 with u [0, 1].\n", - FUNC_NAME, uv[0]); + FUNC_NAME, args->uv[0]); res = RES_BAD_ARG; goto error; } } #else /* SDIS_XD_DIMENSION == 3 */ { - const double w = CLAMP(1 - uv[0] - uv[1], 0, 1); - if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1 - || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) { + const double w = CLAMP(1 - args->uv[0] - args->uv[1], 0, 1); + if(args->uv[0] < 0 || args->uv[1] < 0 || args->uv[0] > 1 || args->uv[1] > 1 + || !eq_eps(w + args->uv[0] + args->uv[1], 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates [%g, %g]. " "u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n", - FUNC_NAME, uv[0], uv[1]); + FUNC_NAME, args->uv[0], args->uv[1]); res = RES_BAD_ARG; goto error; } @@ -122,9 +118,15 @@ XD(solve_probe_boundary) #endif /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC @@ -160,6 +162,8 @@ XD(solve_probe_boundary) } /* Here we go! Launch the Monte Carlo estimation */ + nrealisations = args->nrealisations; + register_paths = out_estimator ? args->register_paths : SDIS_HEAT_PATH_NONE; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -185,7 +189,7 @@ XD(solve_probe_boundary) time_current(&t0); if(!out_green) { - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); if(register_paths) { heat_path_init(scn->dev->allocator, &heat_path); pheat_path = &heat_path; @@ -202,8 +206,9 @@ XD(solve_probe_boundary) pgreen_path = &green_path; } - res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side, - fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w); + res_simul = XD(boundary_realisation)(scn, rng, args->iprim, args->uv, time, + args->side, args->fp_to_meter, args->ambient_radiative_temperature, + args->reference_temperature, pgreen_path, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -191,16 +191,16 @@ main(int argc, char** argv) struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; + struct sdis_solve_probe_boundary_args solve_args = + SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; struct interf* interf_props = NULL; struct fluid* fluid_param; - double uv[2]; double pos[3]; double time_range[2] = { INF, INF }; double tr[2]; double ref; size_t prims[4]; enum sdis_side sides[4]; - size_t iprim; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -297,39 +297,55 @@ main(int argc, char** argv) #define SOLVE sdis_solve_probe_boundary #define GREEN sdis_solve_probe_boundary_green_function - #define F SDIS_FRONT - uv[0] = 0.3; - uv[1] = 0.3; - iprim = 6; - - BA(SOLVE(NULL, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, 0, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, 12, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, NULL, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, NULL, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, -1, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - tr[0] = 1; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - tr[1] = 0; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - OK(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos)); + solve_args.nrealisations = N; + solve_args.uv[0] = 0.3; + solve_args.uv[1] = 0.3; + solve_args.iprim = 6; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.side = SDIS_FRONT; + + BA(SOLVE(NULL, &solve_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.nrealisations = N; + solve_args.iprim = 12; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.iprim = 6; + solve_args.side = SDIS_SIDE_NULL__; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.side = SDIS_FRONT; + solve_args.time_range[0] = solve_args.time_range[1] = -1; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.time_range[0] = 1; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.time_range[1] = 0; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.time_range[0] = solve_args.time_range[1] = INF; + + OK(SOLVE(box_scn, &solve_args, &estimator)); + OK(sdis_scene_get_boundary_position + (box_scn, solve_args.iprim, solve_args.uv, pos)); printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos)); check_estimator(estimator, N, ref); - BA(GREEN(NULL, N, iprim, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, 0, iprim, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, 12, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, NULL, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, -1, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, F, 0, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, NULL)); + BA(GREEN(NULL, &solve_args, &green)); + BA(GREEN(box_scn, NULL, &green)); + BA(GREEN(box_scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(GREEN(box_scn, &solve_args, &green)); + solve_args.nrealisations = N; + solve_args.iprim = 12; + BA(GREEN(box_scn, &solve_args, &green)); + solve_args.iprim = 6; + solve_args.side = SDIS_SIDE_NULL__; + BA(GREEN(box_scn, &solve_args, &green)); + solve_args.side = SDIS_FRONT; + OK(GREEN(box_scn, &solve_args, &green)); - OK(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, &green)); check_green_function(green); OK(sdis_green_function_solve(green, time_range, &estimator2)); check_estimator(estimator2, N, ref); @@ -339,27 +355,29 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); /* Dump paths */ - OK(SOLVE(box_scn, N_dump, iprim, uv, time_range, F, 1.0, 0, 0, - SDIS_HEAT_PATH_ALL, &estimator)); + solve_args.nrealisations = N_dump; + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + OK(SOLVE(box_scn, &solve_args, &estimator)); dump_heat_paths(fp, estimator); OK(sdis_estimator_ref_put(estimator)); /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = UNKNOWN_TEMPERATURE; - BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); + BA(SOLVE(box_scn, &solve_args, &estimator)); fluid_param->temperature = Tf; - uv[0] = 0.5; - iprim = 3; - BA(SOLVE(square_scn, N, 4, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos)); - printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos)); - check_estimator(estimator, N, ref); + solve_args.nrealisations = N; + solve_args.register_paths = SDIS_HEAT_PATH_NONE; + solve_args.uv[0] = 0.5; + solve_args.iprim = 4; + + BA(SOLVE(square_scn, &solve_args, &estimator)); + solve_args.iprim = 3; + OK(SOLVE(square_scn, &solve_args, &estimator)); - OK(GREEN(square_scn, N, iprim, uv, F, 1, 0, 0, &green)); + OK(GREEN(square_scn, &solve_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_estimator_ref_put(estimator)); @@ -368,7 +386,7 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = UNKNOWN_TEMPERATURE; - BA(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); + BA(SOLVE(square_scn, &solve_args, &estimator)); fluid_param->temperature = Tf; #undef F