commit dcf75a5fddbafcfe70339375f9fc874e0519623a
parent 5cd8799035ee1985539e17f10c612deb994b6576
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 12 Mar 2019 11:18:14 +0100
Add and test sdis_solve_probe_boundary_green_function
Diffstat:
11 files changed, 185 insertions(+), 35 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -793,7 +793,7 @@ sdis_green_function_for_each_path
/* Retrieve the spatio-temporal end point of a path used to estimate the green
* function. Note that this point went back in time from the relative
- * observation time 0. Its time is thus negative ; its absolute value
+ * observation time 0. Its time is thus negative; its absolute value
* represents the time spent by the path into the system. */
SDIS_API res_T
sdis_green_path_get_limit_point
@@ -950,7 +950,7 @@ sdis_solve_medium
* 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 invocations : others interfaces/media are
+ * sdis_green_function_solve invocations: others interfaces/media are
* definitely registered against the green function as interfaces/media with no
* flux/volumic power.
*
@@ -968,6 +968,18 @@ sdis_solve_probe_green_function
struct sdis_green_function** green);
SDIS_API res_T
+sdis_solve_probe_boundary_green_function
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ const enum sdis_side side, /* Side of iprim on which the probe lies */
+ 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_green.c b/src/sdis_green.c
@@ -774,6 +774,27 @@ error:
}
res_T
+green_function_redux_and_clear
+ (struct sdis_green_function* dst,
+ struct sdis_green_function* greens[],
+ const size_t ngreens)
+{
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(dst && greens && ngreens);
+
+ FOR_EACH(i, 0, ngreens) {
+ res = green_function_merge_and_clear(dst, greens[i]);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
green_function_finalize
(struct sdis_green_function* green,
struct ssp_rng_proxy* proxy)
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -42,6 +42,12 @@ green_function_merge_and_clear
(struct sdis_green_function* dst,
struct sdis_green_function* src);
+extern LOCAL_SYM res_T
+green_function_redux_and_clear
+ (struct sdis_green_function* dst,
+ struct sdis_green_function* greens[],
+ const size_t ngreens);
+
/* Finalize the green function state (e.g.: computes the #paths & #failures,
* save the rng state, etc.) */
extern LOCAL_SYM res_T
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 green_path_handle* green_path,
struct sdis_heat_path* heat_path,
double* weight);
@@ -93,6 +94,7 @@ boundary_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct green_path_handle* green_path,
struct sdis_heat_path* heat_path,
double* weight);
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -213,6 +213,7 @@ XD(boundary_realisation)
const double fp_to_meter,
const double Tarad,
const double Tref,
+ struct green_path_handle* green_path,
struct sdis_heat_path* heat_path,
double* weight)
{
@@ -262,6 +263,7 @@ XD(boundary_realisation)
SDIS_HEAT_VERTEX_CONDUCTION);
if(res != RES_OK) goto error;
+ ctx.green_path = green_path;
ctx.heat_path = heat_path;
ctx.Tarad = Tarad;
ctx.Tref3 = Tref*Tref*Tref;
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -235,10 +235,32 @@ sdis_solve_probe_boundary
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, register_paths, out_estimator);
+ side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator);
} else {
return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, time_range,
- side, fp_to_meter, Tarad, Tref, register_paths, out_estimator);
+ side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator);
+ }
+}
+
+res_T
+sdis_solve_probe_boundary_green_function
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ const enum sdis_side side, /* Side of iprim on which the probe lies */
+ 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_probe_boundary_2d(scn, nrealisations, iprim, uv, NULL,
+ side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL);
+ } else {
+ return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, NULL,
+ side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL);
}
}
@@ -300,16 +322,14 @@ sdis_solve_boundary_flux
const double Tref, /* In Kelvin */
struct sdis_estimator** out_estimator)
{
- res_T res = RES_OK;
if(!scn) return RES_BAD_ARG;
if(scene_is_2d(scn)) {
- res = solve_boundary_flux_2d(scn, nrealisations, primitives, nprimitives,
+ return solve_boundary_flux_2d(scn, nrealisations, primitives, nprimitives,
time_range, fp_to_meter, Tarad, Tref, out_estimator);
} else {
- res = solve_boundary_flux_3d(scn, nrealisations, primitives, nprimitives,
+ return solve_boundary_flux_3d(scn, nrealisations, primitives, nprimitives,
time_range, fp_to_meter, Tarad, Tref, out_estimator);
}
- return res;
}
res_T
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -234,7 +234,7 @@ XD(solve_boundary)
/* Invoke the boundary realisation */
res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
- fp_to_meter, Tarad, Tref, pheat_path, &w);
+ fp_to_meter, Tarad, Tref, NULL, pheat_path, &w);
/* Update the MC accumulators */
if(res_local == RES_OK) {
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -367,12 +367,11 @@ XD(solve_medium)
}
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*/
- FOR_EACH(i, 1, scn->dev->nthreads) { /* Merge the per thread green */
- res = green_function_merge_and_clear(green, greens[i]);
- if(res != RES_OK) goto error;
- }
+ 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);
@@ -399,6 +398,10 @@ exit:
if(out_green) *out_green = green;
return (res_T)res;
error:
+ if(green) {
+ SDIS(green_function_ref_put(green));
+ green = NULL;
+ }
if(estimator) {
SDIS(estimator_ref_put(estimator));
estimator = NULL;
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -96,15 +96,14 @@ XD(solve_probe)
res = scene_get_medium(scn, position, NULL, &medium);
if(res != RES_OK) goto error;
+ /* Create the per thread green function */
if(out_green) {
- /* Create the per thread green function */
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 */
@@ -182,12 +181,11 @@ XD(solve_probe)
}
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*/
- FOR_EACH(i, 1, scn->dev->nthreads) { /* Merge the per thread green */
- res = green_function_merge_and_clear(green, greens[i]);
- if(res != RES_OK) goto error;
- }
+ 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);
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -40,9 +40,12 @@ XD(solve_probe_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 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;
@@ -51,13 +54,21 @@ XD(solve_probe_boundary)
ATOMIC res = RES_OK;
if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv
- || !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 || (side != SDIS_FRONT && side != SDIS_BACK)
- || !out_estimator) {
+ || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)) {
+ 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; }
@@ -123,9 +134,21 @@ XD(solve_probe_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;
+ }
/* Here we go! Launch the Monte Carlo estimation */
omp_set_num_threads((int)scn->dev->nthreads);
@@ -134,6 +157,8 @@ XD(solve_probe_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;
double w = NaN;
@@ -142,15 +167,23 @@ XD(solve_probe_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;
}
res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
- fp_to_meter, Tarad, Tref, pheat_path, &w);
+ fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w);
if(res_local == RES_OK) {
accum->sum += w;
accum->sum2 += w*w;
@@ -177,9 +210,24 @@ XD(solve_probe_boundary)
if(res != RES_OK) goto error;
/* Setup the estimated temperature */
- sum_accums(accums, scn->dev->nthreads, &accums[0]);
- 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);
+ }
+
+ 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) {
@@ -188,11 +236,22 @@ exit:
}
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(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:
+ if(green) {
+ SDIS(green_function_ref_put(green));
+ green = NULL;
+ }
if(estimator) {
SDIS(estimator_ref_put(estimator));
estimator = NULL;
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -181,6 +181,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 = SDIS_INTERFACE_SHADER_NULL;
@@ -199,7 +201,7 @@ main(int argc, char** argv)
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
- OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev));
/* Temporary file used to dump heat paths */
CHK(fp = tmpfile());
@@ -291,6 +293,7 @@ main(int argc, char** argv)
ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
#define SOLVE sdis_solve_probe_boundary
+ #define GREEN sdis_solve_probe_boundary_green_function
#define F SDIS_FRONT
uv[0] = 0.3;
uv[1] = 0.3;
@@ -314,7 +317,23 @@ main(int argc, char** argv)
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);
+
+ BA(GREEN(NULL, N, iprim, uv, F, 1, 0, 0, &green));
+ BA(GREEN(box_scn, 0, iprim, uv, F, 1, 0, 0, &green));
+ BA(GREEN(box_scn, N, 12, uv, F, 1, 0, 0, &green));
+ BA(GREEN(box_scn, N, iprim, NULL, F, 1, 0, 0, &green));
+ BA(GREEN(box_scn, N, iprim, uv, -1, 1, 0, 0, &green));
+ BA(GREEN(box_scn, N, iprim, uv, F, 0, 0, 0, &green));
+ BA(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, NULL));
+
+ 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);
+
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Dump paths */
OK(SOLVE(box_scn, N_dump, iprim, uv, time_range, F, 1.0, 0, 0,
@@ -334,7 +353,15 @@ main(int argc, char** argv)
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);
+
+ 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);
+
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 = UNKNOWN_TEMPERATURE;
