commit 5cd8799035ee1985539e17f10c612deb994b6576
parent 1c08c544493d3fd94dae13ef1fd10a84e97049e8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 11 Mar 2019 15:56:01 +0100
Add the register_paths flag to the boundary temperature solvers
Diffstat:
7 files changed, 154 insertions(+), 44 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -855,7 +855,7 @@ 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 */
- const int register_path, /* Combination of enum sdis_heat_path_flag */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** estimator);
SDIS_API res_T
@@ -869,6 +869,7 @@ sdis_solve_probe_boundary
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** estimator);
SDIS_API res_T
@@ -882,6 +883,7 @@ sdis_solve_boundary
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** estimator);
SDIS_API res_T
@@ -931,7 +933,7 @@ sdis_solve_medium
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
- const int register_path, /* Combination of enum sdis_heat_path_flag */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** estimator);
/*******************************************************************************
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -79,6 +79,7 @@ boundary_realisation_2d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct sdis_heat_path* heat_path,
double* weight);
extern LOCAL_SYM res_T
@@ -92,6 +93,7 @@ boundary_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct sdis_heat_path* heat_path,
double* weight);
extern LOCAL_SYM res_T
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -46,10 +46,15 @@ XD(compute_temperature)
size_t istack = 0;
#endif
/* Maximum accepted #failures before stopping the realisation */
+ struct sdis_heat_vertex* heat_vtx = NULL;
const size_t MAX_FAILS = 10;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ if(ctx->heat_path && T->func == XD(boundary_path)) {
+ heat_vtx = heat_path_get_last_vertex(ctx->heat_path);
+ }
+
do {
/* Save the current random walk state */
const struct XD(rwalk) rwalk_bkp = *rwalk;
@@ -71,6 +76,25 @@ XD(compute_temperature)
} while(res == RES_BAD_OP && ++nfails < MAX_FAILS);
if(res != RES_OK) goto error;
+ /* Update the type of the first vertex of the random walks that begin on a
+ * boundary. Indeed, one knows the "right" type of the first vertex only
+ * after the boundary_path execution that defines the sub path to resolve
+ * from the submitted boundary position. Note that if the boundary
+ * temperature is know, the type is let as it. */
+ if(heat_vtx && !T->done) {
+ if(heat_path_get_last_vertex(ctx->heat_path) != heat_vtx) {
+ /* Path was reinjected into a solid */
+ heat_vtx->type = SDIS_HEAT_VERTEX_CONDUCTION;
+ } else if(T->func == XD(convective_path)) {
+ heat_vtx->type = SDIS_HEAT_VERTEX_CONVECTION;
+ } else if(T->func == XD(radiative_path)) {
+ heat_vtx->type = SDIS_HEAT_VERTEX_RADIATIVE;
+ } else {
+ FATAL("Unreachable code.\n");
+ }
+ heat_vtx = NULL; /* Notify that the first vertex is finalized */
+ }
+
} while(!T->done);
exit:
@@ -131,7 +155,7 @@ XD(probe_realisation)
rwalk.vtx.time = time;
/* Register the starting position against the heat path */
- type = medium->type == SDIS_SOLID
+ type = medium->type == SDIS_SOLID
? SDIS_HEAT_VERTEX_CONDUCTION
: SDIS_HEAT_VERTEX_CONVECTION;
res = register_heat_vertex(heat_path, &rwalk.vtx, 0, type);
@@ -189,6 +213,7 @@ XD(boundary_realisation)
const double fp_to_meter,
const double Tarad,
const double Tref,
+ struct sdis_heat_path* heat_path,
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -233,14 +258,23 @@ XD(boundary_realisation)
f2_set(rwalk.hit.uv, st);
#endif
+ res = register_heat_vertex(heat_path, &rwalk.vtx, 0/*weight*/,
+ SDIS_HEAT_VERTEX_CONDUCTION);
+ if(res != RES_OK) goto error;
+
+ ctx.heat_path = heat_path;
ctx.Tarad = Tarad;
ctx.Tref3 = Tref*Tref*Tref;
res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
- if(res != RES_OK) return res;
+ if(res != RES_OK) goto error;
*weight = T.value;
- return RES_OK;
+
+exit:
+ return res;
+error:
+ goto exit;
}
res_T
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -229,15 +229,16 @@ sdis_solve_probe_boundary
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** out_estimator)
{
if(!scn) return RES_BAD_ARG;
if(scene_is_2d(scn)) {
return solve_probe_boundary_2d(scn, nrealisations, iprim, uv, time_range,
- side, fp_to_meter, Tarad, Tref, out_estimator);
+ side, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
} else {
return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, time_range,
- side, fp_to_meter, Tarad, Tref, out_estimator);
+ side, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
}
}
@@ -252,15 +253,16 @@ sdis_solve_boundary
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** out_estimator)
{
if(!scn) return RES_BAD_ARG;
if(scene_is_2d(scn)) {
return solve_boundary_2d(scn, nrealisations, primitives, sides, nprimitives,
- time_range, fp_to_meter, Tarad, Tref, out_estimator);
+ time_range, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
} else {
return solve_boundary_3d(scn, nrealisations, primitives, sides, nprimitives,
- time_range, fp_to_meter, Tarad, Tref, out_estimator);
+ time_range, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
}
}
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -83,6 +83,7 @@ XD(solve_boundary)
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** out_estimator)
{
struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL);
@@ -187,6 +188,8 @@ XD(solve_boundary)
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
struct accum* accum = &accums[ithread];
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
struct sXd(primitive) prim;
enum sdis_side side;
size_t iprim;
@@ -200,6 +203,11 @@ XD(solve_boundary)
time = sample_time(rng, time_range);
+ 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
@@ -225,8 +233,8 @@ XD(solve_boundary)
side = sides[prim.prim_id];
/* Invoke the boundary realisation */
- res_local = XD(boundary_realisation)
- (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w);
+ res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
+ fp_to_meter, Tarad, Tref, pheat_path, &w);
/* Update the MC accumulators */
if(res_local == RES_OK) {
@@ -237,6 +245,20 @@ XD(solve_boundary)
ATOMIC_SET(&res, res_local);
continue;
}
+
+ if(pheat_path) {
+ pheat_path->status = res_local == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res_local = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+ }
+ }
}
sum_accums(accums, scn->dev->nthreads, &accums[0]);
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -39,6 +39,7 @@ XD(solve_probe_boundary)
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** out_estimator)
{
struct sdis_estimator* estimator = NULL;
@@ -133,6 +134,8 @@ XD(solve_probe_boundary)
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
struct accum* accum = &accums[ithread];
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
double w = NaN;
double time;
res_T res_local;
@@ -141,8 +144,13 @@ XD(solve_probe_boundary)
time = sample_time(rng, time_range);
- res_local = XD(boundary_realisation)
- (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w);
+ if(register_paths) {
+ heat_path_init(scn->dev->allocator, &heat_path);
+ pheat_path = &heat_path;
+ }
+
+ res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
+ fp_to_meter, Tarad, Tref, pheat_path, &w);
if(res_local == RES_OK) {
accum->sum += w;
accum->sum2 += w*w;
@@ -151,6 +159,20 @@ XD(solve_probe_boundary)
ATOMIC_SET(&res, res_local);
continue;
}
+
+ if(pheat_path) {
+ pheat_path->status = res_local == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res_local = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+ }
+ }
}
if(res != RES_OK) goto error;
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -44,6 +44,7 @@
#define UNKNOWN_TEMPERATURE -1
#define N 10000 /* #realisations */
+#define N_dump 10 /* #dumped paths */
#define Tf 310.0
#define Tb 300.0
@@ -168,6 +169,7 @@ check_estimator
int
main(int argc, char** argv)
{
+ FILE* fp = NULL;
struct mem_allocator allocator;
struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
@@ -199,6 +201,9 @@ main(int argc, char** argv)
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+ /* Temporary file used to dump heat paths */
+ CHK(fp = tmpfile());
+
/* Create the fluid medium */
OK(sdis_data_create
(dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
@@ -291,35 +296,41 @@ main(int argc, char** argv)
uv[1] = 0.3;
iprim = 6;
- BA(SOLVE(NULL, N, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, 0, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, 12, uv, time_range, F, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, iprim, NULL, time_range, F, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, iprim, uv, NULL, F, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, iprim, uv, time_range, -1, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, NULL));
+ BA(SOLVE(NULL, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, 0, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, 12, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, iprim, NULL, time_range, F, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, iprim, uv, NULL, F, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, iprim, uv, time_range, -1, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, NULL));
tr[0] = tr[1] = -1;
- BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL));
tr[0] = 1;
- BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL));
tr[1] = 0;
- BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL));
- OK(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
+ OK(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
OK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos));
printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos));
check_estimator(estimator, N, ref);
OK(sdis_estimator_ref_put(estimator));
+ /* Dump paths */
+ OK(SOLVE(box_scn, N_dump, iprim, uv, time_range, F, 1.0, 0, 0,
+ SDIS_HEAT_PATH_ALL, &estimator));
+ dump_heat_paths(fp, estimator);
+ OK(sdis_estimator_ref_put(estimator));
+
/* The external fluid cannot have an unknown temperature */
fluid_param->temperature = UNKNOWN_TEMPERATURE;
- BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
fluid_param->temperature = Tf;
uv[0] = 0.5;
iprim = 3;
- BA(SOLVE(square_scn, N, 4, uv, time_range, F, 1.0, 0, 0, &estimator));
- OK(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
+ BA(SOLVE(square_scn, N, 4, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
+ OK(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
OK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos));
printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos));
check_estimator(estimator, N, ref);
@@ -327,8 +338,9 @@ main(int argc, char** argv)
/* The external fluid cannot have an unknown temperature */
fluid_param->temperature = UNKNOWN_TEMPERATURE;
- BA(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, &estimator));
+ BA(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator));
fluid_param->temperature = Tf;
+
#undef F
#undef SOLVE
@@ -340,39 +352,51 @@ main(int argc, char** argv)
#define SOLVE sdis_solve_boundary
prims[0] = 6;
prims[1] = 7;
- BA(SOLVE(NULL, N, prims, sides, 2, time_range, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, 0, prims, sides, 2, time_range, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, NULL, sides, 2, time_range, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, prims, NULL, 2, time_range, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, prims, sides, 0, time_range, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, prims, sides, 2, NULL, 1.0, 0, 0, &estimator));
- BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, NULL));
+ BA(SOLVE(NULL, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, 0, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, NULL, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, prims, NULL, 2, time_range, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, prims, sides, 0, time_range, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, prims, sides, 2, NULL, 1.0, 0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, NULL));
tr[0] = tr[1] = -1;
- BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL));
tr[0] = 1;
- BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL));
tr[1] = 0;
- BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, NULL));
+ BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL));
/* Average temperature on the right side of the box */
- OK(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, &estimator));
+ OK(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
printf("Average temperature of the right side of the box = ");
check_estimator(estimator, N, ref);
OK(sdis_estimator_ref_put(estimator));
+ /* Dump path */
+ OK(SOLVE(box_scn, N_dump, prims, sides, 2, time_range, 1.0, 0, 0,
+ SDIS_HEAT_PATH_ALL, &estimator));
+ dump_heat_paths(fp, estimator);
+ OK(sdis_estimator_ref_put(estimator));
+
/* Average temperature on the right side of the square */
prims[0] = 3;
sides[0] = SDIS_FRONT;
- OK(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, &estimator));
+ OK(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, 0, &estimator));
printf("Average temperature of the right side of the square = ");
check_estimator(estimator, N, ref);
OK(sdis_estimator_ref_put(estimator));
+ /* Dump path */
+ OK(SOLVE(square_scn, N_dump, prims, sides, 1, time_range, 1.0, 0, 0,
+ SDIS_HEAT_PATH_ALL, &estimator));
+ dump_heat_paths(fp, estimator);
+ OK(sdis_estimator_ref_put(estimator));
+
/* Check out of bound prims */
prims[0] = 12;
- BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, &estimator));
+ BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
prims[0] = 4;
- BA(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, &estimator));
+ BA(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, 0, &estimator));
/* Average temperature on the left+right sides of the box */
prims[0] = 2;
@@ -382,7 +406,7 @@ main(int argc, char** argv)
ref = (ref + Tb) / 2;
- OK(SOLVE(box_scn, N, prims, sides, 4, time_range, 1.0, 0, 0, &estimator));
+ OK(SOLVE(box_scn, N, prims, sides, 4, time_range, 1.0, 0, 0, 0, &estimator));
printf("Average temperature of the left+right sides of the box = ");
check_estimator(estimator, N, ref);
OK(sdis_estimator_ref_put(estimator));
@@ -390,7 +414,7 @@ main(int argc, char** argv)
/* Average temperature on the left+right sides of the square */
prims[0] = 1;
prims[1] = 3;
- OK(SOLVE(square_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, &estimator));
+ OK(SOLVE(square_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
printf("Average temperature of the left+right sides of the square = ");
check_estimator(estimator, N, ref);
OK(sdis_estimator_ref_put(estimator));
@@ -400,6 +424,8 @@ main(int argc, char** argv)
OK(sdis_scene_ref_put(square_scn));
OK(sdis_device_ref_put(dev));
+ CHK(fclose(fp) == 0);
+
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
CHK(mem_allocated_size() == 0);
