commit 1d6c1d47a62c5b6ccdf2a275339df7948652a123
parent 81dc3fd6a23f9e5541590932cc286d191bc5ad74
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 15 Feb 2019 14:57:15 +0100
Test the probe green function with flux and volumic power
Diffstat:
6 files changed, 230 insertions(+), 90 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -686,10 +686,21 @@ sdis_solve_camera
*
* 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.
+ * 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
+ * sdis_green_function_solve invokations : other interfaces/media are
+ * definitely registered against the green function as interfaces/media with no
+ * flux/volumic power.
+ *
+ * If the aforementionned 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
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -642,7 +642,7 @@ green_path_add_power_term
size_t iterm;
unsigned id;
res_T res = RES_OK;
- ASSERT(handle && mdm && vtx && val >= 0);
+ ASSERT(handle && mdm && vtx);
/* Unused position and time: the current implementation of the green function
* assumes that the power is constant in space and time per medium. */
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -34,6 +34,8 @@ XD(conductive_path)
struct XD(temperature)* T)
{
double position_start[DIM];
+ double green_power_factor = 0;
+ double power_ref = SDIS_VOLUMIC_POWER_NONE;
struct sdis_medium* mdm;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
@@ -51,6 +53,12 @@ XD(conductive_path)
/* Save the submitted position */
dX(set)(position_start, rwalk->vtx.P);
+ if(ctx->green_path) {
+ /* Retrieve the power of the medium. Use it to check that it is effectively
+ * constant along the random walk */
+ power_ref = solid_get_volumic_power(mdm, &rwalk->vtx);
+ }
+
do { /* Solid random walk */
struct get_medium_info info = GET_MEDIUM_INFO_NULL;
struct sXd(hit) hit0, hit1;
@@ -75,7 +83,7 @@ XD(conductive_path)
res = green_path_set_limit_vertex(ctx->green_path, rwalk->mdm, &rwalk->vtx);
if(res != RES_OK) goto error;
}
- goto exit;
+ break;
}
/* Fetch solid properties */
@@ -85,6 +93,15 @@ XD(conductive_path)
cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
power = solid_get_volumic_power(mdm, &rwalk->vtx);
+ if(ctx->green_path && power_ref != power) {
+ log_err(scn->dev,
+ "%s: invalid non constant volumic power term. Expecting a constant "
+ "volumic power in time and space on green function estimation.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
#if DIM == 2
/* Sample a direction around 2PI */
ssp_ran_circle_uniform_float(rng, dir0, NULL);
@@ -159,11 +176,9 @@ XD(conductive_path)
}
}
- /* Register the power term against the green function */
+ /* Register the power term for the green function */
if(ctx->green_path && power != SDIS_VOLUMIC_POWER_NONE) {
- res = green_path_add_power_term
- (ctx->green_path, mdm, &rwalk->vtx, power_factor);
- if(res != RES_OK) goto error;
+ green_power_factor += power_factor;
}
/* Sample the time */
@@ -179,7 +194,7 @@ XD(conductive_path)
if(tmp >= 0) {
T->value += tmp;
T->done = 1;
- goto exit;
+ break;
}
/* The initial condition should have been reached */
log_err(scn->dev,
@@ -245,6 +260,13 @@ XD(conductive_path)
/* Keep going while the solid random walk does not hit an interface */
} while(SXD_HIT_NONE(&rwalk->hit));
+ /* Register the power term for the green function */
+ if(ctx->green_path && power_ref != SDIS_VOLUMIC_POWER_NONE) {
+ res = green_path_add_power_term
+ (ctx->green_path, mdm, &rwalk->vtx, green_power_factor);
+ if(res != RES_OK) goto error;
+ }
+
T->func = XD(boundary_path);
rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -16,6 +16,7 @@
#include "sdis.h"
#include "test_sdis_utils.h"
+#include <rsys/clock_time.h>
#include <rsys/double3.h>
/*
@@ -123,13 +124,98 @@ interface_get_flux
}
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+solve(struct sdis_scene* scn, const double pos[])
+{
+ char dump[128];
+ struct time t0, t1, t2;
+ struct sdis_estimator* estimator;
+ struct sdis_estimator* estimator2;
+ struct sdis_green_function* green;
+ struct sdis_mc T;
+ size_t nreals;
+ size_t nfails;
+ double ref;
+ const double time_range[2] = {INF, INF};
+ enum sdis_scene_dimension dim;
+ ASSERT(scn && pos);
+
+ ref = T0 + (1 - pos[0]) * PHI/LAMBDA;
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &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_scene_get_dimension(scn, &dim));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), 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\n", dump);
+
+ 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, N, pos, 1.0, 0, 0, &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) {
+ case SDIS_SCENE_2D:
+ printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), ref, 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\n", dump);
+
+ check_estimator_eq(estimator, estimator2);
+
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ OK(sdis_green_function_ref_put(green));
+}
+
+/*******************************************************************************
* Test
******************************************************************************/
int
main(int argc, char** argv)
{
struct mem_allocator allocator;
- struct sdis_mc T = SDIS_MC_NULL;
struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
@@ -139,7 +225,6 @@ main(int argc, char** argv)
struct sdis_interface* interf_phi = NULL;
struct sdis_scene* box_scn = NULL;
struct sdis_scene* square_scn = NULL;
- struct sdis_estimator* estimator = 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 = SDIS_INTERFACE_SHADER_NULL;
@@ -147,10 +232,6 @@ main(int argc, char** argv)
struct sdis_interface* square_interfaces[4/*#segments*/];
struct interf* interf_props = NULL;
double pos[3];
- double time_range[2] = { INF, INF };
- double ref;
- size_t nreals;
- size_t nfails;
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
@@ -230,34 +311,12 @@ main(int argc, char** argv)
OK(sdis_interface_ref_put(interf_T0));
OK(sdis_interface_ref_put(interf_phi));
+ /* Solve */
d3_splat(pos, 0.25);
- ref = T0 + (1 - pos[0]) * PHI/LAMBDA;
-
- /* Solve in 3D */
- 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));
- OK(sdis_estimator_ref_put(estimator));
- printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.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, T.SE*3));
-
- /* Solve in 2D */
- 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));
- OK(sdis_estimator_ref_put(estimator));
- printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.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, T.SE*3));
+ printf(">> Box scene\n");
+ solve(box_scn, pos);
+ printf(">> Square Scene\n");
+ solve(square_scn, pos);
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
@@ -267,5 +326,4 @@ main(int argc, char** argv)
mem_shutdown_proxy_allocator(&allocator);
CHK(mem_allocated_size() == 0);
return 0;
-
}
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -244,7 +244,6 @@ 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);
@@ -252,14 +251,6 @@ check_estimator_eq
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);
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -16,6 +16,7 @@
#include "sdis.h"
#include "test_sdis_utils.h"
+#include <rsys/clock_time.h>
#include <rsys/double3.h>
#include <rsys/math.h>
@@ -145,13 +146,100 @@ interface_get_convection_coef
}
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+solve(struct sdis_scene* scn, const double pos[])
+{
+ char dump[128];
+ struct time t0, t1, t2;
+ struct sdis_estimator* estimator;
+ struct sdis_estimator* estimator2;
+ struct sdis_green_function* green;
+ struct sdis_mc T;
+ size_t nreals;
+ size_t nfails;
+ double x;
+ double ref;
+ const double time_range[2] = {INF, INF};
+ enum sdis_scene_dimension dim;
+ ASSERT(scn && pos);
+
+ x = pos[0] - 0.5;
+ ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0;
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &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_scene_get_dimension(scn, &dim));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), 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\n", dump);
+
+ 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, N, pos, 1.0, 0, 0, &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) {
+ case SDIS_SCENE_2D:
+ printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), ref, 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\n", dump);
+
+ check_estimator_eq(estimator, estimator2);
+
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ OK(sdis_green_function_ref_put(green));
+}
+
+/*******************************************************************************
* Test
******************************************************************************/
int
main(int argc, char** argv)
{
struct mem_allocator allocator;
- struct sdis_mc T = SDIS_MC_NULL;
struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
@@ -161,7 +249,6 @@ main(int argc, char** argv)
struct sdis_interface* interf_T0 = NULL;
struct sdis_scene* box_scn = NULL;
struct sdis_scene* square_scn = NULL;
- struct sdis_estimator* estimator = 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 = SDIS_INTERFACE_SHADER_NULL;
@@ -170,11 +257,6 @@ main(int argc, char** argv)
struct interf* interf_props = NULL;
struct solid* solid_props = NULL;
double pos[3];
- double time_range[2] = { INF, INF };
- double x;
- double ref;
- size_t nreals;
- size_t nfails;
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
@@ -263,36 +345,12 @@ main(int argc, char** argv)
OK(sdis_interface_ref_put(interf_adiabatic));
OK(sdis_interface_ref_put(interf_T0));
+ /* Solve */
d3_splat(pos, 0.25);
- x = pos[0] - 0.5;
- 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, &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("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.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));
-
- /* Solve in 2D */
- 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));
- OK(sdis_estimator_ref_put(estimator));
- printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.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));
+ printf(">> Box scene\n");
+ solve(box_scn, pos);
+ printf(">> Square scene\n");
+ solve(square_scn, pos);
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
