commit 19dc1374381595c7505ffc4c7c0ea39bb96e9bbb parent f980c7743b18510d20cf926bbc539e84816d9d2d Author: Christophe Coustet <christophe.coustet@meso-star.com> Date: Mon, 7 Sep 2020 16:06:44 +0200 Change Green API to allow unsteady computations Diffstat:
26 files changed, 530 insertions(+), 402 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h @@ -988,7 +988,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); /* Retrieve the number of valid paths used to estimate the green function. It @@ -1133,23 +1132,11 @@ sdis_compute_mean_power /******************************************************************************* * Green solvers. * - * The caller should ensure that green solvers are invoked on scenes whose data - * do not depend on time. Indeed, on green estimation, the time parameter along - * the random walks registers the relative time spent in the system rather than - * an absolute time. As a consequence, the media/interfaces parameters cannot - * vary in time with respect to an absolute time value. - * - * 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 values 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. - * - * If the aforementioned assumptions are not ensured by the caller, the - * behavior of the estimated green function is undefined. ******************************************************************************/ SDIS_API res_T sdis_solve_probe_green_function @@ -1178,4 +1165,3 @@ sdis_solve_medium_green_function END_DECLS #endif /* SDIS_H */ - diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -276,7 +276,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) { @@ -294,7 +293,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 */ @@ -326,26 +324,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->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: @@ -442,7 +430,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; @@ -454,8 +441,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; } @@ -477,10 +463,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_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -211,7 +211,7 @@ XD(convective_path) double mu, tau, t0; mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V; tau = ssp_ran_exp(rng, mu); - t0 = ctx->green_path ? -INF : fluid_get_t0(rwalk->mdm); + t0 = fluid_get_t0(rwalk->mdm); rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); /* Register the new vertex against the heat path */ 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, @@ -49,7 +50,7 @@ XD(time_rewind) cp = solid_get_calorific_capacity(mdm, &rwalk->vtx); /* Fetch the limit time */ - t0 = ctx->green_path ? -INF : solid_get_t0(mdm); + t0 = solid_get_t0(mdm); /* Sample the time to reroll */ mu = (2*DIM*lambda)/(rho*cp*delta_in_meter*delta_in_meter); @@ -82,6 +83,11 @@ XD(time_rewind) vtx->weight = T->value; } + if(ctx->green_path) { + res = green_path_set_limit_vertex(ctx->green_path, mdm, &rwalk->vtx); + 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. Simply takes 0 as relative time */ - time = 0; + 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. Simply takes 0 as relative time */ - time = 0; + + 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. Simply takes 0 as relative time */ - time = 0; + 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. Simply takes 0 as relative time */ - time = 0; + 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,19 +438,27 @@ 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); @@ -436,7 +466,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 */ p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10; OK(sdis_solve_probe(scn, &solve_args, &estimator)); @@ -448,17 +478,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); @@ -472,7 +511,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,7 +426,7 @@ 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); @@ -440,7 +440,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,7 +297,7 @@ 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); OK(sdis_estimator_ref_put(estimator2)); @@ -338,7 +338,7 @@ 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); OK(sdis_estimator_ref_put(estimator2)); diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -313,7 +313,7 @@ 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); OK(sdis_estimator_ref_put(estimator2)); @@ -353,7 +353,7 @@ 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); OK(sdis_estimator_ref_put(estimator2)); 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,138 +146,194 @@ 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); - - 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); + + 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"); + } } /******************************************************************************* @@ -298,7 +358,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)); @@ -379,15 +439,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,7 +346,7 @@ 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); OK(sdis_green_function_ref_put(green)); @@ -376,7 +376,7 @@ 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); OK(sdis_estimator_ref_put(estimator)); @@ -468,7 +468,7 @@ 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); OK(sdis_green_function_ref_put(green)); @@ -504,7 +504,7 @@ 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); OK(sdis_green_function_ref_put(green)); @@ -536,7 +536,7 @@ 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); OK(sdis_green_function_ref_put(green)); @@ -553,7 +553,7 @@ 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); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -456,7 +456,7 @@ 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); diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -399,7 +399,7 @@ 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); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -436,10 +436,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); @@ -448,7 +448,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)); @@ -466,7 +466,7 @@ 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); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -265,7 +265,7 @@ 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); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -264,7 +264,7 @@ 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); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -321,7 +321,7 @@ 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); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -313,7 +313,7 @@ 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); 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)); 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,141 +173,202 @@ 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); - - 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); + + 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"); + } } /******************************************************************************* @@ -324,7 +393,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)); @@ -349,6 +418,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)); @@ -405,15 +476,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);