stardis-solver

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

commit 1750006d7e20c0984cec509403f8d501d70725c7
parent 141e6beb06d328bd0d8bfd984dd9127628141174
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 25 Jun 2020 16:47:10 +0200

Update the API of sdis_solve_probe_boundary_flux

Diffstat:
Msrc/sdis.h | 32+++++++++++++++++++++++++-------
Msrc/sdis_solve.c | 14+++-----------
Msrc/sdis_solve_probe_boundary_Xd.h | 65++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/test_sdis_solve_boundary_flux.c | 72+++++++++++++++++++++++++++++++++++++++++++-----------------------------
4 files changed, 107 insertions(+), 76 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -430,6 +430,30 @@ struct sdis_solve_medium_args { static const struct sdis_solve_medium_args SDIS_SOLVE_MEDIUM_ARGS_DEFAULT = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT__; +struct sdis_solve_probe_boundary_flux_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 */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + 0, /* Primitive identifier */ \ + {0,0}, /* UV */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Refernce temperature */ \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_probe_boundary_flux_args +SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT = + SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__; + BEGIN_DECLS /******************************************************************************* @@ -984,13 +1008,7 @@ sdis_solve_boundary SDIS_API res_T sdis_solve_probe_boundary_flux (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 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_flux_args* args, struct sdis_estimator** estimator); SDIS_API res_T diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -304,22 +304,14 @@ sdis_solve_boundary_green_function res_T sdis_solve_probe_boundary_flux (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 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_flux_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_flux_2d(scn, nrealisations, iprim, uv, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_probe_boundary_flux_2d(scn, args, out_estimator); } else { - return solve_probe_boundary_flux_3d(scn, nrealisations, iprim, uv, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_probe_boundary_flux_3d(scn, args, out_estimator); } } diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -331,13 +331,7 @@ error: static res_T XD(solve_probe_boundary_flux) (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 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_flux_args* args, struct sdis_estimator** out_estimator) { struct sdis_estimator* estimator = NULL; @@ -352,17 +346,17 @@ XD(solve_probe_boundary_flux) struct accum* acc_fl = NULL; /* Per thread flux accumulator */ struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */ struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */ + size_t nrealisations = 0; int64_t irealisation = 0; size_t i; int progress = 0; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv - || !time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1]) - || fp_to_meter <= 0 || Tref < 0 - || !out_estimator) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || 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]) + || args->fp_to_meter <= 0 || !out_estimator) { res = RES_BAD_ARG; goto error; } @@ -374,12 +368,12 @@ XD(solve_probe_boundary_flux) #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; @@ -387,29 +381,33 @@ XD(solve_probe_boundary_flux) /* Check parametric coordinates */ if(scene_is_2d(scn)) { - 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 { - 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; } } /* Check medium is fluid on one side and solid on the other */ - interf = scene_get_interface(scn, (unsigned)iprim); + interf = scene_get_interface(scn, (unsigned)args->iprim); fmd = interface_get_medium(interf, SDIS_FRONT); bmd = interface_get_medium(interf, SDIS_BACK); if(!fmd || !bmd @@ -422,9 +420,15 @@ XD(solve_probe_boundary_flux) fluid_side = (fmd->type == SDIS_FLUID) ? SDIS_FRONT : SDIS_BACK; /* 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 @@ -449,7 +453,7 @@ XD(solve_probe_boundary_flux) /* Prebuild the interface fragment */ res = XD(build_interface_fragment) - (&frag, scn, (unsigned)iprim, uv, fluid_side); + (&frag, scn, (unsigned)args->iprim, args->uv, fluid_side); if(res != RES_OK) goto error; /* Create the estimator */ @@ -457,6 +461,7 @@ XD(solve_probe_boundary_flux) if(res != RES_OK) goto error; /* Here we go! Launch the Monte Carlo estimation */ + nrealisations = args->nrealisations; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -471,6 +476,8 @@ XD(solve_probe_boundary_flux) double time, epsilon, hc, hr; int flux_mask = 0; double T_brf[3] = { 0, 0, 0 }; + const double Tref = args->reference_temperature; + const double Tarad = args->ambient_radiative_temperature; size_t n; int pcent; res_T res_simul = RES_OK; @@ -480,7 +487,7 @@ XD(solve_probe_boundary_flux) /* Begin time registration */ time_current(&t0); - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); /* Compute hr and hc */ frag.time = time; @@ -492,8 +499,8 @@ XD(solve_probe_boundary_flux) flux_mask = 0; if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE; if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE; - res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time, - solid_side, fp_to_meter, Tarad, Tref, flux_mask, T_brf); + res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv, time, + solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c @@ -234,16 +234,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_flux_args solve_args = + SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; struct interf* interf_props = NULL; struct fluid* fluid_param; enum sdis_estimator_type type; - double uv[2]; double pos[3]; - double time_range[2] = { INF, INF }; + double time_range[2] = {INF, INF}; double tr[2]; double analyticT, analyticCF, analyticRF, analyticTF; size_t prims[2]; - size_t iprim; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -345,46 +345,60 @@ main(int argc, char** argv) analyticTF = analyticCF + analyticRF; #define SOLVE sdis_solve_probe_boundary_flux - uv[0] = 0.3; - uv[1] = 0.3; - iprim = 6; - - BA(SOLVE(NULL, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, 0, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, 12, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, NULL, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, NULL, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[0] = 1; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[1] = 0; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[1] = INF; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); + solve_args.nrealisations = N; + solve_args.iprim = 6; + solve_args.uv[0] = 0.3; + solve_args.uv[1] = 0.3; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.ambient_radiative_temperature = Trad; + solve_args.reference_temperature = Tref; + 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.uv[0] = solve_args.uv[1] = 1; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.uv[0] = solve_args.uv[1] = 0.3; + 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[1] = INF; + BA(SOLVE(box_scn, &solve_args, &estimator)); + solve_args.time_range[0] = INF; + OK(SOLVE(box_scn, &solve_args, &estimator)); - OK(SOLVE(box_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); OK(sdis_estimator_get_type(estimator, &type)); CHK(type == SDIS_ESTIMATOR_FLUX); - OK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos)); + OK(sdis_scene_get_boundary_position + (box_scn, solve_args.iprim, solve_args.uv, pos)); printf("Boundary values of the box at (%g %g %g) = ", SPLIT3(pos)); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); - uv[0] = 0.5; - iprim = 3; - BA(SOLVE(square_scn, N, 4, uv, time_range, 1.0, Trad, Tref, &estimator)); - OK(SOLVE(square_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - OK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos)); + 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(sdis_scene_get_boundary_position + (square_scn, solve_args.iprim, solve_args.uv, pos)); printf("Boundary values of the square at (%g %g) = ", SPLIT2(pos)); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); #undef F #undef SOLVE - + #define SOLVE sdis_solve_boundary_flux prims[0] = 6; prims[1] = 7;