stardis-solver

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

commit 820cbbcafe0c3d922a7cf7ac8663e222228ddbea
parent 1750006d7e20c0984cec509403f8d501d70725c7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 25 Jun 2020 17:13:34 +0200

Update the API of sdis_solve_boundary_flux

Diffstat:
Msrc/sdis.h | 42++++++++++++++++++++++++++++++------------
Msrc/sdis_solve.c | 14+++-----------
Msrc/sdis_solve_boundary_Xd.h | 65++++++++++++++++++++++++++++++++++++++---------------------------
Msrc/test_sdis_solve_boundary_flux.c | 140++++++++++++++++++++++++++++++++++++++++++++-----------------------------------
4 files changed, 149 insertions(+), 112 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -382,8 +382,8 @@ SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT = struct sdis_solve_boundary_args { size_t nrealisations; /* #realisations */ - size_t* primitives; /* List of boundary primitives to handle */ - enum sdis_side* sides; /* Per primitive side to consider */ + const size_t* primitives; /* List of boundary primitives to handle */ + const enum sdis_side* sides; /* Per primitive side to consider */ size_t nprimitives; /* #primitives */ double time_range[2]; /* Observation time */ double fp_to_meter; /* Scale from floating point units to meters */ @@ -400,7 +400,7 @@ struct sdis_solve_boundary_args { {DBL_MAX,DBL_MAX}, /* Time range */ \ 1, /* FP to meter */ \ -1, /* Ambient radiative temperature */ \ - -1, /* Refernce temperature */ \ + -1, /* Reference temperature */ \ SDIS_HEAT_PATH_NONE, \ NULL /* RNG state */ \ } @@ -423,7 +423,7 @@ struct sdis_solve_medium_args { {DBL_MAX,DBL_MAX}, /* Time range */ \ 1, /* FP to meter */ \ -1, /* Ambient radiative temperature */ \ - -1, /* Refernce temperature */ \ + -1, /* Reference temperature */ \ SDIS_HEAT_PATH_NONE, \ NULL /* RNG state */ \ } @@ -447,13 +447,37 @@ struct sdis_solve_probe_boundary_flux_args { {DBL_MAX,DBL_MAX}, /* Time range */ \ 1, /* FP to meter */ \ -1, /* Ambient radiative temperature */ \ - -1, /* Refernce temperature */ \ + -1, /* Reference 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__; +struct sdis_solve_boundary_flux_args { + size_t nrealisations; /* #realisations */ + const size_t* primitives; /* List of boundary primitives to handle */ + size_t nprimitives; /* #primitives */ + 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_BOUNDARY_FLUX_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + NULL, /* List or primitive ids */ \ + 0, /* #primitives */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_boundary_flux_args +SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT = + SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT__; + BEGIN_DECLS /******************************************************************************* @@ -1014,13 +1038,7 @@ sdis_solve_probe_boundary_flux SDIS_API res_T sdis_solve_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - 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_boundary_flux_args* args, struct sdis_estimator** estimator); SDIS_API res_T diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -318,22 +318,14 @@ sdis_solve_probe_boundary_flux res_T sdis_solve_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - 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_boundary_flux_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_boundary_flux_2d(scn, nrealisations, primitives, nprimitives, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_boundary_flux_2d(scn, args, out_estimator); } else { - return solve_boundary_flux_3d(scn, nrealisations, primitives, nprimitives, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_boundary_flux_3d(scn, args, out_estimator); } } diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -430,13 +430,7 @@ error: static res_T XD(solve_boundary_flux) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - 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_boundary_flux_args* args, struct sdis_estimator** out_estimator) { struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL); @@ -452,17 +446,24 @@ XD(solve_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; size_t i; size_t view_nprims; - int64_t irealisation; int progress = 0; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives - || !time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1]) - || !nprimitives || fp_to_meter < 0 || Tref < 0 + if(!scn + || !args + || !args->nrealisations + || args->nrealisations > INT64_MAX + || !args->primitives + || 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->nprimitives + || args->fp_to_meter < 0 || !out_estimator) { res = RES_BAD_ARG; goto error; @@ -475,12 +476,12 @@ XD(solve_boundary_flux) #endif SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims)); - FOR_EACH(i, 0, nprimitives) { - if(primitives[i] >= view_nprims) { + FOR_EACH(i, 0, args->nprimitives) { + if(args->primitives[i] >= view_nprims) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)primitives[i], + (unsigned long)args->primitives[i], (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; @@ -497,19 +498,20 @@ XD(solve_boundary_flux) /* Initialise the boundary shape with the triangles/segments of the * submitted primitives */ - ctx.primitives = primitives; + ctx.primitives = args->primitives; ctx.view = scn->sXd(view); vdata.get = XD(boundary_get_position); #if SDIS_XD_DIMENSION == 2 vdata.usage = S2D_POSITION; vdata.type = S2D_FLOAT2; - res = s2d_line_segments_setup_indexed_vertices(shape, (unsigned)nprimitives, - XD(boundary_get_indices), (unsigned)(nprimitives*2), &vdata, 1, &ctx); + res = s2d_line_segments_setup_indexed_vertices(shape, + (unsigned)args->nprimitives, XD(boundary_get_indices), + (unsigned)(args->nprimitives*2), &vdata, 1, &ctx); #else /* DIM == 3 */ vdata.usage = S3D_POSITION; vdata.type = S3D_FLOAT3; - res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprimitives, - XD(boundary_get_indices), (unsigned)(nprimitives*3), &vdata, 1, &ctx); + res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)args->nprimitives, + XD(boundary_get_indices), (unsigned)(args->nprimitives*3), &vdata, 1, &ctx); #endif if(res != RES_OK) goto error; @@ -522,9 +524,15 @@ XD(solve_boundary_flux) if(res != RES_OK) goto error; /* 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(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); @@ -550,6 +558,7 @@ XD(solve_boundary_flux) res = estimator_create(scn->dev, SDIS_ESTIMATOR_FLUX, &estimator); if(res != RES_OK) goto error; + 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) { @@ -567,6 +576,8 @@ XD(solve_boundary_flux) const struct sdis_medium *fmd, *bmd; enum sdis_side solid_side, fluid_side; double T_brf[3] = { 0, 0, 0 }; + const double Tref = args->reference_temperature; + const double Tarad = args->ambient_radiative_temperature; double epsilon, hc, hr; size_t iprim; double uv[DIM - 1]; @@ -583,7 +594,7 @@ XD(solve_boundary_flux) /* Begin time registration */ time_current(&t0); - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); /* Sample a position onto the boundary */ #if SDIS_XD_DIMENSION == 2 @@ -605,8 +616,8 @@ XD(solve_boundary_flux) if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } /* Map from boundary scene to sdis scene */ - ASSERT(prim.prim_id < nprimitives); - iprim = primitives[prim.prim_id]; + ASSERT(prim.prim_id < args->nprimitives); + iprim = args->primitives[prim.prim_id]; interf = scene_get_interface(scn, (unsigned)iprim); fmd = interface_get_medium(interf, SDIS_FRONT); @@ -638,7 +649,7 @@ XD(solve_boundary_flux) 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); + 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,14 +234,14 @@ 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 = + struct sdis_solve_probe_boundary_flux_args probe_args = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; + struct sdis_solve_boundary_flux_args bound_args = + SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT; struct interf* interf_props = NULL; struct fluid* fluid_param; enum sdis_estimator_type type; double pos[3]; - double time_range[2] = {INF, INF}; - double tr[2]; double analyticT, analyticCF, analyticRF, analyticTF; size_t prims[2]; (void)argc, (void)argv; @@ -345,53 +345,53 @@ main(int argc, char** argv) analyticTF = analyticCF + analyticRF; #define SOLVE sdis_solve_probe_boundary_flux - 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)); + probe_args.nrealisations = N; + probe_args.iprim = 6; + probe_args.uv[0] = 0.3; + probe_args.uv[1] = 0.3; + probe_args.time_range[0] = INF; + probe_args.time_range[1] = INF; + probe_args.ambient_radiative_temperature = Trad; + probe_args.reference_temperature = Tref; + BA(SOLVE(NULL, &probe_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)); + BA(SOLVE(box_scn, &probe_args, NULL)); + probe_args.nrealisations = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.nrealisations = N; + probe_args.iprim = 12; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.iprim = 6; + probe_args.uv[0] = probe_args.uv[1] = 1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.uv[0] = probe_args.uv[1] = 0.3; + probe_args.time_range[0] = probe_args.time_range[1] = -1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = 1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[1] = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[1] = INF; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = INF; + OK(SOLVE(box_scn, &probe_args, &estimator)); OK(sdis_estimator_get_type(estimator, &type)); CHK(type == SDIS_ESTIMATOR_FLUX); OK(sdis_scene_get_boundary_position - (box_scn, solve_args.iprim, solve_args.uv, pos)); + (box_scn, probe_args.iprim, probe_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)); - 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)); + probe_args.uv[0] = 0.5; + probe_args.iprim = 4; + BA(SOLVE(square_scn, &probe_args, &estimator)); + probe_args.iprim = 3; + OK(SOLVE(square_scn, &probe_args, &estimator)); OK(sdis_scene_get_boundary_position - (square_scn, solve_args.iprim, solve_args.uv, pos)); + (square_scn, probe_args.iprim, probe_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)); @@ -402,38 +402,52 @@ main(int argc, char** argv) #define SOLVE sdis_solve_boundary_flux prims[0] = 6; prims[1] = 7; - BA(SOLVE(NULL, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, 0, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, NULL, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 0, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 2, NULL, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); - tr[0] = 1; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); - tr[1] = 0; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); + bound_args.nrealisations = N; + bound_args.primitives = prims; + bound_args.nprimitives = 2; + bound_args.time_range[0] = INF; + bound_args.time_range[1] = INF; + bound_args.ambient_radiative_temperature = Trad; + bound_args.reference_temperature = Tref; + + BA(SOLVE(NULL, &bound_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &bound_args, NULL)); + bound_args.primitives = NULL; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.primitives = prims; + bound_args.nprimitives = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.nprimitives = 2; + bound_args.time_range[0] = bound_args.time_range[1] = -1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = 1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[1] = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[1] = INF; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = INF; + prims[0] = 12; + BA(SOLVE(box_scn, &bound_args, &estimator)); + prims[0] = 6; + OK(SOLVE(box_scn, &bound_args, &estimator)); /* Average temperature on the right side of the box */ - OK(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); printf("Average values of the right side of the box = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); /* Average temperature on the right side of the square */ + prims[0] = 4; + bound_args.nprimitives = 1; + BA(SOLVE(square_scn, &bound_args, &estimator)); prims[0] = 3; - OK(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average values of the right side of the square = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); - - /* Check out of bound prims */ - prims[0] = 12; - BA(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - prims[0] = 4; - BA(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); - + /* Average temperature on the left side of the box */ prims[0] = 2; prims[1] = 3; @@ -443,14 +457,16 @@ main(int argc, char** argv) analyticRF = Hrad * (analyticT - Trad); analyticTF = analyticCF + analyticRF; - OK(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); + bound_args.nprimitives = 2; + OK(SOLVE(box_scn, &bound_args, &estimator)); printf("Average values of the left side of the box = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); /* Average temperature on the left/right side of the square */ prims[0] = 1; - OK(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); + bound_args.nprimitives = 1; + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average values of the left side of the square = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator));