commit 0ef1e6334ca62b13278106d321d2263415852324 parent 441021eb210785b60dba400c9b2ba18a95b8d968 Author: Christophe Coustet <christophe.coustet@meso-star.com> Date: Fri, 16 Oct 2020 15:29:12 +0200 Merge branch 'feature_unsteady_green' into develop Diffstat:
26 files changed, 560 insertions(+), 423 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h @@ -993,7 +993,6 @@ sdis_green_function_ref_put SDIS_API res_T sdis_green_function_solve (struct sdis_green_function* green, - const double time_range[2], /* Observation time */ struct sdis_estimator** estimator); SDIS_API res_T @@ -1154,18 +1153,16 @@ sdis_compute_power /******************************************************************************* * Green solvers. * - * Currently only steady computations are supported. As a consequence, the - * observation time is always fixed to infinity. - * - * In addition, the green solvers assumes that the interface fluxes are - * constants in time and space. In the same way the volumic power of the solid - * media must be constant in time and space too. Furthermore, note that only - * the interfaces/media that had a flux/volumic power during green estimation - * can update their flux/volumic power value for subsequent + * Note that only the interfaces/media with flux/volumic power defined during + * green estimation can update their flux/volumic power value for subsequent * sdis_green_function_solve invocations: others interfaces/media are * definitely registered against the green function as interfaces/media with no * flux/volumic power. * + * Also note that the green solvers assume that the interface fluxes are + * constant in time and space. The same applies to the volumic power of the + * solid media. + * * If these assumptions are not ensured by the caller, the behavior of the * estimated green function is undefined. ******************************************************************************/ @@ -1196,4 +1193,3 @@ sdis_solve_medium_green_function END_DECLS #endif /* SDIS_H */ - diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -411,7 +411,6 @@ green_function_fetch_interf static res_T green_function_solve_path (struct sdis_green_function* green, - const double time, /* Sampled time */ const size_t ipath, double* weight) { @@ -429,7 +428,6 @@ green_function_solve_path size_t i, n; res_T res = RES_OK; ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight); - ASSERT(time > 0); path = darray_green_path_cdata_get(&green->paths) + ipath; if(path->limit_type == SDIS_POINT_NONE) { /* Rejected path */ @@ -461,26 +459,16 @@ green_function_solve_path /* Setup time. */ switch(path->limit_type) { case SDIS_FRAGMENT: - time_curr = time + path->limit.fragment.time; + time_curr = path->limit.fragment.time; interf = green_function_fetch_interf(green, path->limit_id); break; case SDIS_VERTEX: - time_curr = time + path->limit.vertex.time; + time_curr = path->limit.vertex.time; medium = green_function_fetch_medium(green, path->limit_id); break; default: FATAL("Unreachable code.\n"); break; } - if(time_curr <= 0 - || (path->limit_type == SDIS_VERTEX && time_curr <= medium_get_t0(medium))) { - log_err(green->scn->dev, - "%s: invalid observation time \"%g\": the initial condition is reached " - "while instationary system are not supported by the green function.\n", - FUNC_NAME, time); - res = RES_BAD_ARG; - goto error; - } - /* Compute limit condition */ switch(path->limit_type) { case SDIS_FRAGMENT: @@ -851,7 +839,6 @@ sdis_green_function_ref_put(struct sdis_green_function* green) res_T sdis_green_function_solve (struct sdis_green_function* green, - const double time_range[2], struct sdis_estimator** out_estimator) { struct sdis_estimator* estimator = NULL; @@ -863,8 +850,7 @@ sdis_green_function_solve double accum2 = 0; res_T res = RES_OK; - if(!green || !time_range || time_range[0] < 0 - || time_range[1] < time_range[0] || !out_estimator) { + if(!green || !out_estimator) { res = RES_BAD_ARG; goto error; } @@ -886,10 +872,9 @@ sdis_green_function_solve /* Solve the green function */ FOR_EACH(ipath, 0, npaths) { - const double time = sample_time(rng, time_range); double w; - res = green_function_solve_path(green, time, ipath, &w); + res = green_function_solve_path(green, ipath, &w); if(res == RES_BAD_OP) continue; if(res != RES_OK) goto error; diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -136,7 +136,7 @@ register_heat_vertex extern LOCAL_SYM res_T time_rewind_2d - (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */ + (struct sdis_medium* mdm, /* Medium into which the time is rewinded */ struct ssp_rng* rng, const double delta, const double fp_to_meter, @@ -146,7 +146,7 @@ time_rewind_2d extern LOCAL_SYM res_T time_rewind_3d - (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */ + (struct sdis_medium* mdm, /* Medium into which the time is rewinded */ struct ssp_rng* rng, const double delta, const double fp_to_meter, diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h @@ -17,6 +17,7 @@ #include "sdis_log.h" #include "sdis_medium_c.h" #include "sdis_misc.h" +#include "sdis_green.h" #include <star/ssp.h> @@ -24,7 +25,7 @@ res_T XD(time_rewind) - (const struct sdis_medium* mdm, + (struct sdis_medium* mdm, struct ssp_rng* rng, const double delta, const double fp_to_meter, @@ -84,6 +85,12 @@ XD(time_rewind) vtx->weight = T->value; } + if(ctx->green_path) { + res = green_path_set_limit_vertex(ctx->green_path, mdm, &rwalk->vtx, + rwalk->elapsed_time); + if(res != RES_OK) goto error; + } + exit: return res; error: diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -195,6 +195,7 @@ XD(probe_realisation) res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); if(res != RES_OK) goto error; + ASSERT(T.value >= 0); *weight = T.value; exit: diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -111,14 +111,13 @@ XD(solve_boundary) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - if(args->time_range[0] < 0 + 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; - } + || (args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) + { + res = RES_BAD_ARG; + goto error; } #if SDIS_XD_DIMENSION == 2 @@ -256,16 +255,8 @@ XD(solve_boundary) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported yet */ - time = INF; + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -274,6 +265,11 @@ XD(solve_boundary) pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + /* Sample a position onto the boundary */ #if SDIS_XD_DIMENSION == 2 res_local = s2d_scene_view_sample diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h @@ -321,29 +321,22 @@ XD(solve_medium) if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ time_current(&t0); - - if(!out_green) { - /* Sample the time */ - time = sample_time(rng, args->time_range); - - /* Prepare path registration if necessary */ - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported yet */ - time = INF; + + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); goto error_it; } - pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + /* Uniformly Sample an enclosure that surround the submitted medium and * uniformly sample a position into it */ enc = sample_medium_enclosure(&cumul, rng); diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h @@ -62,17 +62,15 @@ XD(solve_probe) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - if(args->time_range[0] < 0 + 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; - } + || (args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) + { + res = RES_BAD_ARG; + goto error; } - #if SDIS_XD_DIMENSION == 2 if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; } #else @@ -153,16 +151,9 @@ XD(solve_probe) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported */ - time = INF; + time = sample_time(rng, args->time_range); + + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -170,6 +161,10 @@ XD(solve_probe) } pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, medium, args->position, time, args->fp_to_meter, diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -61,14 +61,13 @@ XD(solve_probe_boundary) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - if(args->time_range[0] < 0 + 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; - } + || (args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) + { + res = RES_BAD_ARG; + goto error; } #if SDIS_XD_DIMENSION == 2 @@ -188,16 +187,8 @@ XD(solve_probe_boundary) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported */ - time = INF; + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -206,6 +197,11 @@ XD(solve_probe_boundary) pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + 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); diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -126,6 +126,8 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* ctx) ******************************************************************************/ struct solid { double lambda; + double initial_temperature; + double t0; }; static double @@ -168,6 +170,20 @@ solid_get_delta return 1.0/10.0; } +static double +solid_get_temperature +(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + double t0; + CHK(vtx != NULL); + CHK(data != NULL); + t0 = ((const struct solid*)sdis_data_cget(data))->t0; + if(vtx->time > t0) + return UNKNOWN_TEMPERATURE; + else + return ((const struct solid*)sdis_data_cget(data))->initial_temperature; +} + /******************************************************************************* * Interface ******************************************************************************/ @@ -269,8 +285,8 @@ main(int argc, char** argv) struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_scene* scn = NULL; struct ssp_rng* rng = NULL; - const size_t nsimuls = 4; - size_t isimul; + const int nsimuls = 4; + int isimul; const double emissivity = 1;/* Emissivity of the side +/-X of the solid */ const double lambda = 0.1; /* Conductivity of the solid */ const double Tref = 300; /* Reference temperature */ @@ -292,11 +308,13 @@ main(int argc, char** argv) OK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); ((struct solid*)sdis_data_get(data))->lambda = lambda; + ((struct solid*)sdis_data_get(data))->t0 = 0; + ((struct solid*)sdis_data_get(data))->initial_temperature = (T0 + T1) / 2; solid_shader.calorific_capacity = solid_get_calorific_capacity; solid_shader.thermal_conductivity = solid_get_thermal_conductivity; solid_shader.volumic_mass = solid_get_volumic_mass; solid_shader.delta_solid = solid_get_delta; - solid_shader.temperature = temperature_unknown; + solid_shader.temperature = solid_get_temperature; OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); @@ -389,22 +407,26 @@ main(int argc, char** argv) struct sdis_estimator* estimator2; struct sdis_green_function* green; struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; - double ref, u; + double ref = -1; size_t nreals = 0; size_t nfails = 0; const size_t N = 10000; double T1b; + int steady = (isimul % 2) == 0; /* Reset temperature */ p_intface->temperature = T1; solve_args.nrealisations = N; - solve_args.time_range[0] = INF; - solve_args.time_range[1] = INF; + if(steady) + solve_args.time_range[0] = solve_args.time_range[1] = INF; + else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); solve_args.position[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); - u = (solve_args.position[0] + 1) / thickness; solve_args.reference_temperature = Tref; OK(sdis_solve_probe(scn, &solve_args, &estimator)); @@ -416,28 +438,36 @@ main(int argc, char** argv) tmp = lambda / (2 * lambda + thickness * hr) * (T1 - T0); Ts0 = T0 + tmp; Ts1 = T1 - tmp; - ref = u * Ts1 + (1-u) * Ts0; - printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + if(steady) { + double u = (solve_args.position[0] + 1) / thickness; + ref = u * Ts1 + (1 - u) * Ts0; + printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" + " and T1=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + p_intface->temperature, T.E, T.SE); + } printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); printf("\n"); - /* Check green used at a different temperature */ + /* Check same green used at a different temperature */ p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10; OK(sdis_solve_probe(scn, &solve_args, &estimator)); @@ -449,17 +479,26 @@ main(int argc, char** argv) tmp = lambda / (2 * lambda + thickness * hr) * (T1b - T0); Ts0 = T0 + tmp; Ts1 = T1b - tmp; - ref = u * Ts1 + (1 - u) * Ts0; - printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + + if(steady) { + double u = (solve_args.position[0] + 1) / thickness; + ref = u * Ts1 + (1 - u) * Ts0; + printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" + " and T1=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + p_intface->temperature, T.E, T.SE); + } printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -473,7 +512,7 @@ main(int argc, char** argv) OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); - printf("\n"); + printf("\n\n"); } /* Release memory */ diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -426,10 +426,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); @@ -441,7 +441,7 @@ main(int argc, char** argv) OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); - printf("\n"); + printf("\n\n"); } /* Release memory */ diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -297,10 +297,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, box_scn, solve_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } @@ -339,10 +339,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, square_scn, solve_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -313,10 +313,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, box_scn, solve_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } @@ -354,10 +354,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, square_scn, solve_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -18,13 +18,14 @@ #include <rsys/clock_time.h> #include <rsys/double3.h> +#include <star/ssp.h> /* * The scene is composed of a solid cube/square whose temperature is unknown. * The temperature is fixed at T0 on the +X face. The Flux of the -X face is * fixed to PHI. The flux on the other faces is null (i.e. adiabatic). This * test computes the temperature of a probe position pos into the solid and - * check that it is equal to: + * check that at t=inf it is equal to: * * T(pos) = T0 + (A-pos) * PHI/LAMBDA * @@ -94,7 +95,10 @@ solid_get_temperature { (void)data; CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; + if(vtx->time > 0) + return UNKNOWN_TEMPERATURE; + else + return T0; } /******************************************************************************* @@ -129,7 +133,7 @@ interface_get_flux static void solve (struct sdis_scene* scn, - const double pos[], + struct ssp_rng* rng, struct interf* interf) { char dump[128]; @@ -142,139 +146,195 @@ solve struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals; size_t nfails; - double ref; - const double time_range[2] = {INF, INF}; + double ref = -1; + const int nsimuls = 4; + int isimul; enum sdis_scene_dimension dim; - ASSERT(scn && pos && interf); - - /* Restore phi value */ - interf->phi = PHI; - - ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; + ASSERT(scn && rng && interf); OK(sdis_scene_get_dimension(scn, &dim)); - - solve_args.nrealisations = N; - solve_args.position[0] = pos[0]; - solve_args.position[1] = pos[1]; - solve_args.position[2] = dim == SDIS_SCENE_2D ? 0 : pos[2]; - solve_args.time_range[0] = INF; - solve_args.time_range[1] = INF; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - switch(dim) { + FOR_EACH(isimul, 0, nsimuls) { + int steady = (isimul % 2) == 0; + + /* Restore phi value */ + interf->phi = PHI; + + solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[2] = + dim == SDIS_SCENE_2D ? 0 : ssp_rng_uniform_double(rng, 0.1, 0.9); + + solve_args.nrealisations = N; + if(steady) + solve_args.time_range[0] = solve_args.time_range[1] = INF; + else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, T.SE*3)); - - time_current(&t0); - OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - time_current(&t1); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_current(&t2); - - OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); - OK(sdis_estimator_get_failure_count(estimator2, &nfails)); - OK(sdis_estimator_get_temperature(estimator2, &T)); - - switch(dim) { + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); + + time_current(&t0); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + time_current(&t1); + OK(sdis_green_function_solve(green, &estimator2)); + time_current(&t2); + + OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); + OK(sdis_estimator_get_failure_count(estimator2, &nfails)); + OK(sdis_estimator_get_temperature(estimator2, &T)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady Green temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean Green temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady Green temperature at (%g, %g, %g) with Phi=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean Green temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - time_sub(&t0, &t1, &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green estimation time = %s\n", dump); - time_sub(&t1, &t2, &t1); - time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); - - check_green_function(green); - check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, time_range); - - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - printf("\n"); - - /* Check green used at a different phi */ - interf->phi = 3 * PHI; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; - - switch (dim) { - case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); - break; - case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); - break; - default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + time_sub(&t0, &t1, &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green estimation time = %s\n", dump); + time_sub(&t1, &t2, &t1); + time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + check_green_serialization(green, scn); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check same green used at a different flux value */ + interf->phi = 3 * PHI; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + switch(dim) { + case SDIS_SCENE_2D: + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } + break; + case SDIS_SCENE_3D: + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g, %g) with Phi=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, T.SE * 3)); + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); - time_current(&t0); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); + time_current(&t0); + OK(sdis_green_function_solve(green, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); - check_green_function(green); - check_estimator_eq(estimator, estimator2); + check_green_function(green); + check_estimator_eq(estimator, estimator2); - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - OK(sdis_green_function_ref_put(green)); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); - printf("\n"); + printf("\n\n"); + } } /******************************************************************************* @@ -299,7 +359,7 @@ main(int argc, char** argv) struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; - double pos[3]; + struct ssp_rng* rng = NULL; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -380,15 +440,16 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf_phi)); /* Solve */ - d3_splat(pos, 0.25); + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); printf(">> Box scene\n"); - solve(box_scn, pos, interf_props); + solve(box_scn, rng, interf_props); printf(">> Square Scene\n"); - solve(square_scn, pos, interf_props); + solve(square_scn, rng, interf_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(ssp_rng_ref_put(rng)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -346,9 +346,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &probe_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, probe_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -377,9 +377,9 @@ main(int argc, char** argv) OK(GREEN(square_scn, &probe_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, probe_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); @@ -469,9 +469,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, bound_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -506,9 +506,9 @@ main(int argc, char** argv) OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, bound_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -539,9 +539,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, bound_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -557,9 +557,9 @@ main(int argc, char** argv) OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, bound_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -456,10 +456,10 @@ main(int argc, char** argv) solve_args.fp_to_meter = 1; OK(sdis_solve_medium_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -399,10 +399,10 @@ main(int argc, char** argv) BA(sdis_solve_medium_green_function(scn, &solve_args, NULL)); OK(sdis_solve_medium_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -438,10 +438,10 @@ main(int argc, char** argv) solve_args.fp_to_meter = 1; OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - BA(sdis_green_function_solve(NULL, solve_args.time_range, &estimator2)); - BA(sdis_green_function_solve(green, NULL, &estimator2)); - BA(sdis_green_function_solve(green, solve_args.time_range, NULL)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + BA(sdis_green_function_solve(NULL, &estimator2)); + BA(sdis_green_function_solve(green, NULL)); + BA(sdis_green_function_solve(NULL, NULL)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -450,7 +450,7 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); printf("\n"); - /* Check green used at a different temperature */ + /* Check same green used at a different temperature */ fluid_param->temperature = 500; OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -468,10 +468,10 @@ main(int argc, char** argv) CHK(nfails < N / 1000); CHK(eq_eps(T.E, ref, T.SE)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - + stream = tmpfile(); CHK(stream); BA(sdis_green_function_write(NULL, stream)); @@ -491,7 +491,7 @@ main(int argc, char** argv) OK(sdis_green_function_create_from_stream(scn, stream, &green)); CHK(!fclose(stream)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator3)); + OK(sdis_green_function_solve(green, &estimator3)); check_green_function(green); check_estimator_eq_strict(estimator2, estimator3); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -266,10 +266,10 @@ main(int argc, char** argv) /* Check green */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -264,10 +264,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -321,10 +321,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -313,10 +313,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -221,7 +221,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE)); OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c @@ -117,7 +117,6 @@ solve_green_path(struct sdis_green_path* path, void* ctx) switch(pt.type) { case SDIS_FRAGMENT: frag = pt.data.itfrag.fragment; - frag.time = INF; OK(sdis_interface_get_shader(pt.data.itfrag.intface, &interf)); data = sdis_interface_get_data(pt.data.itfrag.intface); temp = frag.side == SDIS_FRONT @@ -126,7 +125,6 @@ solve_green_path(struct sdis_green_path* path, void* ctx) break; case SDIS_VERTEX: vtx = pt.data.mdmvert.vertex; - vtx.time = INF; type = sdis_medium_get_type(pt.data.mdmvert.medium); data = sdis_medium_get_data(pt.data.mdmvert.medium); if(type == SDIS_FLUID) { @@ -229,7 +227,7 @@ check_green_function(struct sdis_green_function* green) time_range[0] = time_range[1] = INF; - OK(sdis_green_function_solve(green, time_range, &estimator)); + OK(sdis_green_function_solve(green, &estimator)); BA(sdis_green_function_get_paths_count(NULL, &n)); BA(sdis_green_function_get_paths_count(green, NULL)); @@ -367,15 +365,14 @@ dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator) void check_green_serialization (struct sdis_green_function* green, - struct sdis_scene* scn, - const double time_range[2]) + struct sdis_scene* scn) { FILE* stream = NULL; struct sdis_estimator *e1 = NULL; struct sdis_estimator *e2 = NULL; struct sdis_green_function* green2 = NULL; - CHK(green && time_range); + CHK(green && scn); stream = tmpfile(); CHK(stream); @@ -386,8 +383,8 @@ check_green_serialization CHK(!fclose(stream)); check_green_function(green2); - OK(sdis_green_function_solve(green, time_range, &e1)); - OK(sdis_green_function_solve(green2, time_range, &e2)); + OK(sdis_green_function_solve(green, &e1)); + OK(sdis_green_function_solve(green2, &e2)); check_estimator_eq_strict(e1, e2); OK(sdis_estimator_ref_put(e1)); diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -354,8 +354,7 @@ dump_heat_paths extern LOCAL_SYM void check_green_serialization (struct sdis_green_function* green, - struct sdis_scene* scn, - const double time_range[2]); + struct sdis_scene* scn); #endif /* TEST_SDIS_UTILS_H */ diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -19,13 +19,14 @@ #include <rsys/clock_time.h> #include <rsys/double3.h> #include <rsys/math.h> +#include <star/ssp.h> /* * The scene is composed of a solid cube/square whose temperature is unknown. * The convection coefficient with the surrounding fluid is null everywhere The * Temperature of the +/- X faces are fixed to T0, and the solid has a volumic * power of P0. This test computes the temperature of a probe position pos into - * the solid and check that it is is equal to: + * the solid and check that at t=inf it is is equal to: * * T(pos) = P0 / (2*LAMBDA) * (A^2/4 - (pos-0.5)^2) + T0 * @@ -61,6 +62,8 @@ struct solid { double cp; double delta; double vpower; + double initial_temperature; + double t0; }; static double @@ -108,9 +111,14 @@ static double solid_get_temperature (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) { - (void)data; + double t0; CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; + CHK(data != NULL); + t0 = ((const struct solid*)sdis_data_cget(data))->t0; + if(vtx->time > t0) + return UNKNOWN_TEMPERATURE; + else + return ((const struct solid*)sdis_data_cget(data))->initial_temperature; } static double @@ -152,7 +160,7 @@ interface_get_convection_coef static void solve (struct sdis_scene* scn, - const double pos[], + struct ssp_rng* rng, struct solid* solid) { char dump[128]; @@ -165,142 +173,203 @@ solve struct sdis_mc time = SDIS_MC_NULL; size_t nreals; size_t nfails; - double x; - double ref; - const double time_range[2] = {INF, INF}; + double ref = -1; enum sdis_scene_dimension dim; - ASSERT(scn && pos && solid); - - /* Restore power value */ - solid->vpower = P0; - - x = pos[0] - 0.5; - ref = solid->vpower / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; + const int nsimuls = 4; + int isimul; + ASSERT(scn && rng && solid); OK(sdis_scene_get_dimension(scn, &dim)); - - solve_args.nrealisations = N; - solve_args.time_range[0] = INF; - solve_args.time_range[1] = INF; - solve_args.position[0] = pos[0]; - solve_args.position[1] = pos[1]; - solve_args.position[2] = dim == SDIS_SCENE_2D ? 0 : pos[2]; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - switch(dim) { + FOR_EACH(isimul, 0, nsimuls) { + int steady = (isimul % 2) == 0; + + /* Restore power value */ + solid->vpower = P0; + + solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[2] = + dim == SDIS_SCENE_2D ? 0 : ssp_rng_uniform_double(rng, 0.1, 0.9); + + solve_args.nrealisations = N; + if(steady) + solve_args.time_range[0] = solve_args.time_range[1] = INF; + else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g)th Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, T.SE*3)); - - /* Check green function */ - time_current(&t0); - OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - time_current(&t1); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_current(&t2); - - OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); - OK(sdis_estimator_get_failure_count(estimator2, &nfails)); - OK(sdis_estimator_get_temperature(estimator2, &T)); - - switch(dim) { + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); + + /* Check green function */ + time_current(&t0); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + time_current(&t1); + OK(sdis_green_function_solve(green, &estimator2)); + time_current(&t2); + + OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); + OK(sdis_estimator_get_failure_count(estimator2, &nfails)); + OK(sdis_estimator_get_temperature(estimator2, &T)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) with Power=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Green Steady temperature at (%g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Green Mean temperature at (%g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Green Steady temperature at (%g, %g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Green Mean temperature at (%g, %g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - time_sub(&t0, &t1, &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green estimation time = %s\n", dump); - time_sub(&t1, &t2, &t1); - time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); - - check_green_function(green); - check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, time_range); - - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - printf("\n"); - - /* Check green used at a different power */ - solid->vpower = 3 * P0; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; - - switch (dim) { - case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); - break; - case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); - break; - default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + time_sub(&t0, &t1, &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green estimation time = %s\n", dump); + time_sub(&t1, &t2, &t1); + time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + check_green_serialization(green, scn); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check same green used at a different power level */ + solid->vpower = 3 * P0; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + switch(dim) { + case SDIS_SCENE_2D: + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } + break; + case SDIS_SCENE_3D: + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); - CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, T.SE * 3)); + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); - time_current(&t0); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); + time_current(&t0); + OK(sdis_green_function_solve(green, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); - check_green_function(green); - check_estimator_eq(estimator, estimator2); + check_green_function(green); + check_estimator_eq(estimator, estimator2); - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - OK(sdis_green_function_ref_put(green)); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); - printf("\n"); + printf("\n\n"); + } } /******************************************************************************* @@ -325,7 +394,7 @@ main(int argc, char** argv) struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; struct solid* solid_props = NULL; - double pos[3]; + struct ssp_rng* rng = NULL; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -350,6 +419,8 @@ main(int argc, char** argv) solid_props->rho = 25; solid_props->delta = DELTA; solid_props->vpower = P0; + solid_props->t0 = 0; + solid_props->initial_temperature = T0; OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); @@ -406,15 +477,16 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf_T0)); /* Solve */ - d3_splat(pos, 0.25); + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); printf(">> Box scene\n"); - solve(box_scn, pos, solid_props); + solve(box_scn, rng, solid_props); printf(">> Square scene\n"); - solve(square_scn, pos, solid_props); + solve(square_scn, rng, solid_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(ssp_rng_ref_put(rng)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator);