commit d24a683d18a634fa1857856c913fbaddd0d0a0bd parent e825330be6d4ac04900d775cd2eeddec51e44636 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Thu, 14 Feb 2019 16:25:43 +0100 Provide a new solver for probe green function Test the probe green function in several tests. Diffstat:
22 files changed, 333 insertions(+), 71 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h @@ -594,7 +594,7 @@ sdis_green_function_solve struct sdis_estimator** estimator); /******************************************************************************* - * Miscellaneous functions + * Solvers ******************************************************************************/ SDIS_API res_T sdis_solve_probe @@ -605,7 +605,6 @@ sdis_solve_probe const double fp_to_meter, /* Scale from floating point units to meters */ const double ambient_radiative_temperature, /* In Kelvin */ const double reference_temperature, /* In Kelvin */ - struct sdis_green_function** green, /* NULL <=> no green registration */ struct sdis_estimator** estimator); SDIS_API res_T @@ -672,6 +671,26 @@ sdis_solve_camera sdis_write_accums_T writer, void* writer_data); +/******************************************************************************* + * 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, not an + * absolute time. As a consequence, on green estimation, the media/interfaces + * parameters cannot use this parameter to vary in time. If these data vary in + * time, the behavior of the estimated green function is undefined. + ******************************************************************************/ +SDIS_API res_T +sdis_solve_probe_green_function + (struct sdis_scene* scn, + const size_t nrealisations, /* #realisations */ + const double position[3], /* Probe position */ + const double fp_to_meter, /* Scale from floating point units to meters */ + const double ambient_radiative_temperature, /* In Kelvin */ + const double reference_temperature, /* In Kelvin */ + struct sdis_green_function** green); + END_DECLS #endif /* SDIS_H */ diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -32,7 +32,7 @@ #include <limits.h> -enum limit_type { LIMIT_FRAGMENT, LIMIT_VERTEX }; +enum limit_type { LIMIT_FRAGMENT, LIMIT_VERTEX, LIMIT_NONE }; struct power_term { double term; /* Power term computed during green estimation */ @@ -99,7 +99,7 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path) darray_power_term_init(allocator, &path->power_terms); path->limit.vertex = SDIS_RWALK_VERTEX_NULL; path->limit_id = UINT_MAX; - path->limit_type = LIMIT_VERTEX; + path->limit_type = LIMIT_NONE; path->ilast_medium = UINT16_MAX; path->ilast_interf = UINT16_MAX; } @@ -288,6 +288,10 @@ green_function_solve_path ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight); path = darray_green_path_cdata_get(&green->paths) + ipath; + if(path->limit_type == LIMIT_NONE) { /* Rejected path */ + res = RES_BAD_OP; + goto error; + } /* Compute medium power terms */ power = 0; @@ -365,6 +369,7 @@ green_function_clear(struct sdis_green_function* green) struct sdis_medium* medium; medium = *htable_medium_iterator_data_get(&it_medium); SDIS(medium_ref_put(medium)); + htable_medium_iterator_next(&it_medium); } htable_medium_clear(&green->media); @@ -375,6 +380,7 @@ green_function_clear(struct sdis_green_function* green) struct sdis_interface* interf; interf = *htable_interf_iterator_data_get(&it_interf); SDIS(interface_ref_put(interf)); + htable_interf_iterator_next(&it_interf); } htable_interf_clear(&green->interfaces); @@ -427,6 +433,7 @@ sdis_green_function_solve struct ssp_rng* rng = NULL; size_t npaths; size_t ipath; + size_t N = 0; /* #realisations */ double accum = 0; double accum2 = 0; res_T res = RES_OK; @@ -442,22 +449,25 @@ sdis_green_function_solve npaths = darray_green_path_size_get(&green->paths); - /* Create the estimator */ - res = estimator_create - (green->dev, SDIS_ESTIMATOR_TEMPERATURE, npaths, npaths, &estimator); - if(res != RES_OK) goto error; - /* Solve the green function */ FOR_EACH(ipath, 0, npaths) { /* TODO add multi-threading */ const double time = sample_time(rng, time_range); double w; res = green_function_solve_path(green, time, ipath, &w); + if(res == RES_BAD_OP) continue; if(res != RES_OK) goto error; + accum += w; accum2 += w*w; + ++N; } + /* Create the estimator */ + res = estimator_create + (green->dev, SDIS_ESTIMATOR_TEMPERATURE, npaths, N, &estimator); + if(res != RES_OK) goto error; + /* Setup the estimated temperature */ estimator_setup_temperature(estimator, accum, accum2); @@ -548,6 +558,7 @@ green_function_merge_and_clear id = medium_get_id(medium); res = htable_medium_set(&dst->media, &id, &medium); if(res != RES_OK) goto error; + htable_medium_iterator_next(&it_medium); } htable_interf_begin(&src->interfaces, &it_interf); @@ -558,6 +569,7 @@ green_function_merge_and_clear id = interface_get_id(interf); res = htable_interf_set(&dst->interfaces, &id, &interf); if(res != RES_OK) goto error; + htable_interf_iterator_next(&it_interf); } green_function_clear(src); @@ -592,7 +604,7 @@ green_path_set_limit_interface_fragment const struct sdis_interface_fragment* frag) { res_T res = RES_OK; - ASSERT(handle && interf && frag); + ASSERT(handle && interf && frag && handle->path->limit_type == LIMIT_NONE); res = ensure_interface_registration(handle->green, interf); if(res != RES_OK) return res; handle->path->limit.fragment = *frag; @@ -608,7 +620,7 @@ green_path_set_limit_vertex const struct sdis_rwalk_vertex* vert) { res_T res = RES_OK; - ASSERT(handle && mdm && vert); + ASSERT(handle && mdm && vert && handle->path->limit_type == LIMIT_NONE); res = ensure_medium_registration(handle->green, mdm); if(res != RES_OK) return res; handle->path->limit.vertex = *vert; diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -65,7 +65,7 @@ fluid_get_temperature (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { ASSERT(mdm && mdm->type == SDIS_FLUID); - ASSERT(vtx->time >= mdm->shader.fluid.t0); + /*ASSERT(vtx->time >= mdm->shader.fluid.t0);*/ return mdm->shader.fluid.temperature(vtx, mdm->data); } @@ -127,7 +127,7 @@ solid_get_temperature (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { ASSERT(mdm && mdm->type == SDIS_SOLID); - ASSERT(vtx->time >= mdm->shader.solid.t0); + /*ASSERT(vtx->time >= mdm->shader.solid.t0);*/ return mdm->shader.solid.temperature(vtx, mdm->data); } diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -167,16 +167,35 @@ sdis_solve_probe const double fp_to_meter,/* Scale factor from floating point unit to meter */ const double Tarad, /* Ambient radiative temperature */ const double Tref, /* Reference temperature */ - struct sdis_green_function** out_green, /* May be NULL<=>Do not store green */ struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { return solve_probe_2d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, out_green, out_estimator); + fp_to_meter, Tarad, Tref, NULL, out_estimator); } else { return solve_probe_3d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, out_green, out_estimator); + fp_to_meter, Tarad, Tref, NULL, out_estimator); + } +} + +res_T +sdis_solve_probe_green_function + (struct sdis_scene* scn, + const size_t nrealisations, + const double position[3], + const double fp_to_meter,/* Scale factor from floating point unit to meter */ + const double Tarad, /* Ambient radiative temperature */ + const double Tref, /* Reference temperature */ + struct sdis_green_function** out_green) +{ + if(!scn) return RES_BAD_ARG; + if(scene_is_2d(scn)) { + return solve_probe_2d(scn, nrealisations, position, NULL, + fp_to_meter, Tarad, Tref, out_green, NULL); + } else { + return solve_probe_3d(scn, nrealisations, position, NULL, + fp_to_meter, Tarad, Tref, out_green, NULL); } } diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -175,13 +175,22 @@ XD(solve_probe) ATOMIC res = RES_OK; if(!scn || !nrealisations || nrealisations > INT64_MAX || !position - || !time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1]) - || fp_to_meter <= 0 || Tref < 0 || !out_estimator) { + || fp_to_meter <= 0 || Tref < 0) { res = RES_BAD_ARG; goto error; } - + if(!out_estimator && !out_green) { + res = RES_BAD_ARG; + goto error; + } + if(out_estimator) { + if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] + || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + res = RES_BAD_ARG; + goto error; + } + } + #if SDIS_XD_DIMENSION == 2 if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; } #else @@ -254,13 +263,15 @@ XD(solve_probe) } if(res != RES_OK) goto error; - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_TEMPERATURE, nrealisations, N, &estimator); - if(res != RES_OK) goto error; + if(out_estimator) { + /* Create the estimator */ + res = estimator_create + (scn->dev, SDIS_ESTIMATOR_TEMPERATURE, nrealisations, N, &estimator); + if(res != RES_OK) goto error; - /* Setup the estimated temperature */ - estimator_setup_temperature(estimator, weight, sqr_weight); + /* Setup the estimated temperature */ + estimator_setup_temperature(estimator, weight, sqr_weight); + } if(out_green) { green = greens[0]; /* Return the green of the 1st thread */ diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -385,6 +385,8 @@ main(int argc, char** argv) FOR_EACH(isimul, 0, nsimuls) { struct sdis_mc T = SDIS_MC_NULL; struct sdis_estimator* estimator; + struct sdis_estimator* estimator2; + struct sdis_green_function* green; double pos[3]; double time_range[2] = { INF, INF }; double ref, u; @@ -396,7 +398,7 @@ main(int argc, char** argv) pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); pos[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); - OK(sdis_solve_probe(scn, N, pos, time_range, 1, -1, Tref, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1, -1, Tref, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -411,7 +413,14 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, 2*T.SE) == 1); + /* Check green function */ + OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, Tref, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); } /* Release memory */ diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -391,6 +391,8 @@ main(int argc, char** argv) FOR_EACH(isimul, 0, nsimuls) { struct sdis_mc T = SDIS_MC_NULL; struct sdis_estimator* estimator; + struct sdis_estimator* estimator2; + struct sdis_green_function* green; double pos[2]; double time_range[2] = { INF, INF }; double ref, u; @@ -401,7 +403,7 @@ main(int argc, char** argv) pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); - OK(sdis_solve_probe(scn, 10000, pos, time_range, 1, -1, Tref, NULL, &estimator)); + OK(sdis_solve_probe(scn, 10000, pos, time_range, 1, -1, Tref, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -416,7 +418,14 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + /* Check green function */ + OK(sdis_solve_probe_green_function(scn, 10000, pos, 1, -1, Tref, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); } /* Release memory */ diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -69,27 +69,30 @@ ******************************************************************************/ static double fluid_get_temperature - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) { - (void)data; CHK(vtx != NULL); - return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE; + if(*((int*)sdis_data_cget(is_stationary))) { + return UNKNOWN_TEMPERATURE; + } else { + return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE; + } } static double fluid_get_volumic_mass - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) { - (void)data; + (void)is_stationary; CHK(vtx != NULL); return RHO; } static double fluid_get_calorific_capacity - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) { - (void)data; + (void)is_stationary; CHK(vtx != NULL); return CP; } @@ -167,6 +170,7 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_medium* fluid = NULL; struct sdis_medium* solid = NULL; + struct sdis_data* is_stationary = NULL; struct sdis_interface* interf_T0 = NULL; struct sdis_interface* interf_T1 = NULL; struct sdis_interface* interf_T2 = NULL; @@ -176,6 +180,8 @@ main(int argc, char** argv) struct sdis_scene* box_scn = NULL; struct sdis_scene* square_scn = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; @@ -194,10 +200,12 @@ main(int argc, char** argv) OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev)); /* Create the fluid medium */ + OK(sdis_data_create(dev, sizeof(int), ALIGNOF(int), NULL, &is_stationary)); + *((int*)sdis_data_get(is_stationary)) = 0; fluid_shader.temperature = fluid_get_temperature; fluid_shader.calorific_capacity = fluid_get_calorific_capacity; fluid_shader.volumic_mass = fluid_get_volumic_mass; - OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + OK(sdis_fluid_create(dev, &fluid_shader, is_stationary, &fluid)); /* Create the solid_medium */ OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); @@ -265,17 +273,29 @@ main(int argc, char** argv) time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); + /* Setup stationary state */ + *((int*)sdis_data_get(is_stationary)) = IS_INF(time); + /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_ref_put(estimator)); printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE); if(nfails) printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(eq_eps(T.E, ref, T.SE * 3)); + + if(IS_INF(time)) { /* Check green function */ + OK(sdis_solve_probe_green_function(box_scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); + } + + OK(sdis_estimator_ref_put(estimator)); } /* Test in 2D for various time values. */ @@ -288,22 +308,35 @@ main(int argc, char** argv) time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); + /* Setup stationnary state */ + *((int*)sdis_data_get(is_stationary)) = IS_INF(time); + /* Solve in 2D */ - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_ref_put(estimator)); printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE); if(nfails) printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(eq_eps(T.E, ref, T.SE * 3)); + + if(IS_INF(time)) { /* Check green function */ + OK(sdis_solve_probe_green_function(square_scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); + } + + OK(sdis_estimator_ref_put(estimator)); } OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(sdis_data_ref_put(is_stationary)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -75,11 +75,15 @@ ******************************************************************************/ static double fluid_get_temperature - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) { - (void)data; CHK(vtx != NULL); - return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE; + CHK(is_stationary != NULL); + if(*((int*)sdis_data_cget(is_stationary))) { + return UNKNOWN_TEMPERATURE; + } else { + return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE; + } } static double @@ -176,6 +180,7 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_medium* fluid = NULL; struct sdis_medium* solid = NULL; + struct sdis_data* is_stationary = NULL; struct sdis_interface* interf_T0 = NULL; struct sdis_interface* interf_T1 = NULL; struct sdis_interface* interf_T2 = NULL; @@ -185,6 +190,8 @@ main(int argc, char** argv) struct sdis_scene* box_scn = NULL; struct sdis_scene* square_scn = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; @@ -202,11 +209,14 @@ main(int argc, char** argv) OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev)); + OK(sdis_data_create(dev, sizeof(int), ALIGNOF(int), NULL, &is_stationary)); + *((int*)sdis_data_get(is_stationary)) = 0; + /* Create the fluid medium */ fluid_shader.temperature = fluid_get_temperature; fluid_shader.calorific_capacity = fluid_get_calorific_capacity; fluid_shader.volumic_mass = fluid_get_volumic_mass; - OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + OK(sdis_fluid_create(dev, &fluid_shader, is_stationary, &fluid)); /* Create the solid_medium */ OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); @@ -280,17 +290,28 @@ main(int argc, char** argv) time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); + *((int*)sdis_data_get(is_stationary)) = IS_INF(time); + /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_ref_put(estimator)); printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE); if(nfails) printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N); CHK(eq_eps(T.E, ref, T.SE * 3)); + + if(IS_INF(time)) { /* Check green function */ + OK(sdis_solve_probe_green_function(box_scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); + } + + OK(sdis_estimator_ref_put(estimator)); } /* Test in 2D for various time values. */ @@ -303,21 +324,33 @@ main(int argc, char** argv) time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + *((int*)sdis_data_get(is_stationary)) = IS_INF(time); + + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_ref_put(estimator)); printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE); if(nfails) printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N); CHK(eq_eps(T.E, ref, T.SE * 3)); + + if(IS_INF(time)) { /* Check green function */ + OK(sdis_solve_probe_green_function(square_scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); + } + + OK(sdis_estimator_ref_put(estimator)); } OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(sdis_data_ref_put(is_stationary)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -234,7 +234,7 @@ main(int argc, char** argv) ref = T0 + (1 - pos[0]) * PHI/LAMBDA; /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -247,7 +247,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE*3)); /* Solve in 2D */ - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -182,6 +182,8 @@ main(int argc, char** argv) struct sdis_scene* scn = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; @@ -261,13 +263,13 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - BA(sdis_solve_probe(NULL, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); - BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, NULL, &estimator)); - BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, NULL, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, NULL, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, NULL, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, NULL)); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(NULL, N, pos, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); BA(sdis_estimator_get_type(estimator, NULL)); BA(sdis_estimator_get_type(NULL, &type)); @@ -316,7 +318,34 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + + fluid_param->temperature = 300; + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + + BA(sdis_solve_probe_green_function(NULL, N, pos, 1.0, 0, 0, &green)); + BA(sdis_solve_probe_green_function(scn, 0, pos, 1.0, 0, 0, &green)); + BA(sdis_solve_probe_green_function(scn, N, NULL, 1.0, 0, 0, &green)); + BA(sdis_solve_probe_green_function(scn, N, pos, 0.0, 0, 0, &green)); + BA(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, -1, &green)); + BA(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, NULL)); + OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); + + BA(sdis_green_function_solve(NULL, time_range, &estimator2)); + BA(sdis_green_function_solve(green, NULL, &estimator2)); + BA(sdis_green_function_solve(green, time_range, NULL)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + + check_estimator_eq(estimator, estimator2); + + BA(sdis_green_function_ref_get(NULL)); + OK(sdis_green_function_ref_get(green)); + BA(sdis_green_function_ref_put(NULL)); + OK(sdis_green_function_ref_put(green)); + OK(sdis_green_function_ref_put(green)); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -149,12 +149,14 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; struct sdis_medium* solid = NULL; struct sdis_medium* fluid = NULL; struct sdis_interface* Tnone = NULL; struct sdis_interface* T300 = NULL; struct sdis_interface* T350 = NULL; struct sdis_scene* scn = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -240,7 +242,7 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -256,8 +258,15 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, 3*T.SE)); + /* Check green */ + OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + /* Release data */ OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -146,12 +146,14 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; struct sdis_medium* solid = NULL; struct sdis_medium* fluid = NULL; struct sdis_interface* Tnone = NULL; struct sdis_interface* T300 = NULL; struct sdis_interface* T350 = NULL; struct sdis_scene* scn = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -238,7 +240,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -254,9 +256,16 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, T.SE*2)); + /* Check green function */ + OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + /* Release data */ OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); OK(sdis_scene_ref_put(scn)); + OK(sdis_green_function_ref_put(green)); OK(sdis_device_ref_put(dev)); check_memory_allocator(&allocator); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -171,6 +171,7 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; struct sdis_medium* solid = NULL; struct sdis_medium* fluid = NULL; struct sdis_interface* Tnone = NULL; @@ -178,6 +179,7 @@ main(int argc, char** argv) struct sdis_interface* T350 = NULL; struct sdis_interface* solid_solid = NULL; struct sdis_scene* scn = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -296,7 +298,7 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -312,8 +314,15 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, 2*T.SE)); + /* Check green function */ + OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + /* Release data */ OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -168,6 +168,7 @@ main(int argc, char** argv) struct sdis_device* dev = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; struct sdis_medium* solid = NULL; struct sdis_medium* fluid = NULL; struct sdis_interface* Tnone = NULL; @@ -175,6 +176,7 @@ main(int argc, char** argv) struct sdis_interface* T350 = NULL; struct sdis_interface* solid_solid = NULL; struct sdis_scene* scn = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -289,7 +291,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -305,8 +307,15 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, 3*T.SE)); + /* Check green function */ + OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + /* Release data */ OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -142,6 +142,8 @@ main(int argc, char** argv) struct sdis_scene* scn = NULL; struct sdis_data* data = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator2 = NULL; + struct sdis_green_function* green = NULL; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -199,7 +201,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); @@ -214,12 +216,18 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, T.SE)); + OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_green_function_solve(green, time_range, &estimator2)); + check_estimator_eq(estimator, estimator2); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -240,6 +240,50 @@ dump_segments } static INLINE void +check_estimator_eq + (const struct sdis_estimator* e1, const struct sdis_estimator* e2) +{ + struct sdis_mc mc1, mc2; + size_t n1, n2; + enum sdis_estimator_type type1, type2; + ASSERT(e1 && e2); + + OK(sdis_estimator_get_type(e1, &type1)); + OK(sdis_estimator_get_type(e2, &type2)); + CHK(type1 == type2); + + OK(sdis_estimator_get_realisation_count(e1, &n1)); + OK(sdis_estimator_get_realisation_count(e2, &n2)); + CHK(n1 == n2); + + OK(sdis_estimator_get_failure_count(e1, &n1)); + OK(sdis_estimator_get_failure_count(e2, &n2)); + CHK(n1 == n2); + + OK(sdis_estimator_get_temperature(e1, &mc1)); + OK(sdis_estimator_get_temperature(e2, &mc2)); + CHK(mc1.E + mc1.SE >= mc2.E - mc2.SE); + CHK(mc1.E - mc1.SE <= mc2.E + mc2.SE); + + if(type1 == SDIS_ESTIMATOR_FLUX) { + OK(sdis_estimator_get_convective_flux(e1, &mc1)); + OK(sdis_estimator_get_convective_flux(e2, &mc2)); + CHK(mc1.E + mc1.SE >= mc2.E - mc2.SE); + CHK(mc1.E - mc1.SE <= mc2.E + mc2.SE); + + OK(sdis_estimator_get_radiative_flux(e1, &mc1)); + OK(sdis_estimator_get_radiative_flux(e2, &mc2)); + CHK(mc1.E + mc1.SE >= mc2.E - mc2.SE); + CHK(mc1.E - mc1.SE <= mc2.E + mc2.SE); + + OK(sdis_estimator_get_total_flux(e1, &mc1)); + OK(sdis_estimator_get_total_flux(e2, &mc2)); + CHK(mc1.E + mc1.SE >= mc2.E - mc2.SE); + CHK(mc1.E - mc1.SE <= mc2.E + mc2.SE); + } +} + +static INLINE void check_memory_allocator(struct mem_allocator* allocator) { if(MEM_ALLOCATED_SIZE(allocator)) { diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -268,7 +268,7 @@ main(int argc, char** argv) ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -282,7 +282,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE)); /* Solve in 2D */ - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -242,7 +242,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) pos[1] = refs[i].pos[1]; pos[2] = refs[i].pos[2]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -261,7 +261,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) pos[0] = refs[i].pos[0]; pos[1] = refs[i].pos[1]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -441,7 +441,7 @@ main(int argc, char** argv) FATAL("Unreachable code.\n"); } - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c @@ -348,7 +348,7 @@ main(int argc, char** argv) Tref = T2 + (T1-T2)/L * (pos[1] + vertices[3]); #endif - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, NULL, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails));