commit 7db538678420d047796cd092b306de31f860ac4c
parent dcf75a5fddbafcfe70339375f9fc874e0519623a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 12 Mar 2019 14:56:11 +0100
Add and test sdis_solve_boundary_green_function
Diffstat:
4 files changed, 161 insertions(+), 21 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -980,6 +980,18 @@ sdis_solve_probe_boundary_green_function
struct sdis_green_function** green);
SDIS_API res_T
+sdis_solve_boundary_green_function
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const enum sdis_side sides[], /* Per primitive side to consider */
+ const size_t nprimitives, /* #primitives */
+ 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);
+
+SDIS_API res_T
sdis_solve_medium_green_function
(struct sdis_scene* scn,
const size_t nrealisations, /* #realisations */
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -281,10 +281,34 @@ sdis_solve_boundary
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, register_paths, out_estimator);
+ time_range, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator);
} else {
return solve_boundary_3d(scn, nrealisations, primitives, sides, nprimitives,
- time_range, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
+ time_range, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator);
+ }
+}
+
+res_T
+sdis_solve_boundary_green_function
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const enum sdis_side sides[], /* Per primitive side to consider */
+ const size_t nprimitives, /* #primitives */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_green_function** green)
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return solve_boundary_2d(scn, nrealisations, primitives, sides,
+ nprimitives, NULL, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green,
+ NULL);
+ } else {
+ return solve_boundary_3d(scn, nrealisations, primitives, sides,
+ nprimitives, NULL, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green,
+ NULL);
}
}
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -84,6 +84,7 @@ XD(solve_boundary)
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_green_function** out_green,
struct sdis_estimator** out_estimator)
{
struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL);
@@ -92,6 +93,8 @@ XD(solve_boundary)
struct sXd(shape)* shape = NULL;
struct sXd(scene_view)* view = NULL;
struct sdis_estimator* estimator = NULL;
+ struct sdis_green_function* green = NULL;
+ struct sdis_green_function** greens = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
struct accum* accums = NULL;
@@ -101,13 +104,21 @@ XD(solve_boundary)
ATOMIC res = RES_OK;
if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives
- || !time_range || time_range[0] < 0 || time_range[1] < time_range[0]
- || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])
- || !sides || !nprimitives || fp_to_meter < 0 || Tref < 0
- || !out_estimator) {
+ || !sides || !nprimitives || 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; }
@@ -178,9 +189,21 @@ XD(solve_boundary)
accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
if(!accums) { res = RES_MEM_ERR; goto error; }
+ /* Create the per thread green function */
+ if(out_green) {
+ greens = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*greens));
+ if(!greens) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = green_function_create(scn->dev, &greens[i]);
+ if(res != RES_OK) goto error;
+ }
+ }
+
/* Create the estimator */
- res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
- if(res != RES_OK) goto error;
+ if(out_estimator) {
+ res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
+ if(res != RES_OK) goto error;
+ }
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
@@ -188,6 +211,8 @@ XD(solve_boundary)
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
struct accum* accum = &accums[ithread];
+ struct green_path_handle* pgreen_path = NULL;
+ struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
struct sdis_heat_path* pheat_path = NULL;
struct sdis_heat_path heat_path;
struct sXd(primitive) prim;
@@ -201,11 +226,19 @@ XD(solve_boundary)
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
- time = sample_time(rng, time_range);
-
- if(register_paths) {
- heat_path_init(scn->dev->allocator, &heat_path);
- pheat_path = &heat_path;
+ if(!out_green) {
+ time = sample_time(rng, 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;
+ res_local = green_function_create_path(greens[ithread], &green_path);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+ pgreen_path = &green_path;
}
/* Sample a position onto the boundary */
@@ -234,7 +267,7 @@ XD(solve_boundary)
/* Invoke the boundary realisation */
res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
- fp_to_meter, Tarad, Tref, NULL, pheat_path, &w);
+ fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w);
/* Update the MC accumulators */
if(res_local == RES_OK) {
@@ -261,22 +294,46 @@ XD(solve_boundary)
}
}
- sum_accums(accums, scn->dev->nthreads, &accums[0]);
-
/* Setup the estimated temperature */
- estimator_setup_realisations_count(estimator, nrealisations, accums[0].count);
- estimator_setup_temperature(estimator, accums[0].sum, accums[0].sum2);
+ if(out_estimator) {
+ struct accum acc;
+ sum_accums(accums, scn->dev->nthreads, &acc);
+ estimator_setup_realisations_count(estimator, nrealisations, acc.count);
+ estimator_setup_temperature(estimator, acc.sum, acc.sum2);
+ }
+
+ /* Setup the green function */
+ if(out_green) {
+ /* Redux the per thread green function into the green of the 1st thread */
+ green = greens[0]; /* Return the green of the 1st thread */
+ greens[0] = NULL; /* Make invalid the 1st green for 'on exit' clean up*/
+ res = green_function_redux_and_clear(green, greens+1, scn->dev->nthreads-1);
+ if(res != RES_OK) goto error;
+
+ /* Finalize the estimated green */
+ res = green_function_finalize(green, rng_proxy);
+ if(res != RES_OK) goto error;
+ }
exit:
if(rngs) {
- FOR_EACH(i, 0, scn->dev->nthreads) {if(rngs[i]) SSP(rng_ref_put(rngs[i]));}
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
MEM_RM(scn->dev->allocator, rngs);
}
+ if(greens) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(greens[i]) SDIS(green_function_ref_put(greens[i]));
+ }
+ MEM_RM(scn->dev->allocator, greens);
+ }
if(accums) MEM_RM(scn->dev->allocator, accums);
if(scene) SXD(scene_ref_put(scene));
if(shape) SXD(shape_ref_put(shape));
if(view) SXD(scene_view_ref_put(view));
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_green) *out_green = green;
if(out_estimator) *out_estimator = estimator;
return (res_T)res;
error:
@@ -284,6 +341,10 @@ error:
SDIS(estimator_ref_put(estimator));
estimator = NULL;
}
+ if(green) {
+ SDIS(green_function_ref_put(green));
+ green = NULL;
+ }
goto exit;
}
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -329,7 +329,7 @@ main(int argc, char** argv)
OK(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, &green));
check_green_function(green);
OK(sdis_green_function_solve(green, time_range, &estimator2));
- check_estimator_eq(estimator, estimator2);
+ check_estimator(estimator2, N, ref);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
@@ -357,7 +357,7 @@ main(int argc, char** argv)
OK(GREEN(square_scn, N, iprim, uv, F, 1, 0, 0, &green));
check_green_function(green);
OK(sdis_green_function_solve(green, time_range, &estimator2));
- check_estimator_eq(estimator, estimator2);
+ check_estimator(estimator2, N, ref);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
@@ -370,6 +370,7 @@ main(int argc, char** argv)
#undef F
#undef SOLVE
+ #undef GREEN
sides[0] = SDIS_FRONT;
sides[1] = SDIS_FRONT;
@@ -377,6 +378,7 @@ main(int argc, char** argv)
sides[3] = SDIS_FRONT;
#define SOLVE sdis_solve_boundary
+ #define GREEN sdis_solve_boundary_green_function
prims[0] = 6;
prims[1] = 7;
BA(SOLVE(NULL, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator));
@@ -397,7 +399,23 @@ main(int argc, char** argv)
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);
+
+ BA(GREEN(NULL, N, prims, sides, 2, 1.0, 0, 0, &green));
+ BA(GREEN(box_scn, 0, prims, sides, 2, 1.0, 0, 0, &green));
+ BA(GREEN(box_scn, N, NULL, sides, 2, 1.0, 0, 0, &green));
+ BA(GREEN(box_scn, N, prims, NULL, 2, 1.0, 0, 0, &green));
+ BA(GREEN(box_scn, N, prims, sides, 0, 1.0, 0, 0, &green));
+ BA(GREEN(box_scn, N, prims, sides, 2, 0.0, 0, 0, &green));
+ BA(GREEN(box_scn, N, prims, sides, 2, 1.0, 0, 0, NULL));
+
+ OK(GREEN(box_scn, N, prims, sides, 2, 1.0, 0, 0, &green));
+ check_green_function(green);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ check_estimator(estimator2, N, ref);
+
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Dump path */
OK(SOLVE(box_scn, N_dump, prims, sides, 2, time_range, 1.0, 0, 0,
@@ -411,7 +429,15 @@ main(int argc, char** argv)
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(GREEN(square_scn, N, prims, sides, 1, 1.0, 0, 0, &green));
+ check_green_function(green);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ check_estimator(estimator2, N, ref);
+
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Dump path */
OK(SOLVE(square_scn, N_dump, prims, sides, 1, time_range, 1.0, 0, 0,
@@ -436,7 +462,15 @@ main(int argc, char** argv)
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(GREEN(box_scn, N, prims, sides, 4, 1.0, 0, 0, &green));
+ check_green_function(green);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ check_estimator(estimator2, N, ref);
+
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Average temperature on the left+right sides of the square */
prims[0] = 1;
@@ -444,8 +478,17 @@ main(int argc, char** argv)
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(GREEN(square_scn, N, prims, sides, 2, 1.0, 0, 0, &green));
+ check_green_function(green);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ check_estimator(estimator2, N, ref);
+
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
#undef SOLVE
+ #undef GREEN
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
