commit 1c08c544493d3fd94dae13ef1fd10a84e97049e8
parent 8c9afb07f5524a5f4057182f54710804d766c0ee
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 11 Mar 2019 10:39:16 +0100
Add and test the sdis_solve_medium_green_function function
Diffstat:
6 files changed, 165 insertions(+), 35 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -965,6 +965,16 @@ sdis_solve_probe_green_function
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 */
+ struct sdis_medium* medium, /* Medium to solve */
+ 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_solve.c b/src/sdis_solve.c
@@ -460,15 +460,33 @@ sdis_solve_medium
const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_estimator** estimator)
{
- res_T res = RES_OK;
if(!scn) return RES_BAD_ARG;
if(scene_is_2d(scn)) {
- res = solve_medium_2d(scn, nrealisations, medium, time_range, fp_to_meter, Tarad,
- Tref, register_paths, estimator);
+ return solve_medium_2d(scn, nrealisations, medium, time_range, fp_to_meter,
+ Tarad, Tref, register_paths, NULL, estimator);
} else {
- res = solve_medium_3d(scn, nrealisations, medium, time_range, fp_to_meter, Tarad,
- Tref, register_paths, estimator);
+ return solve_medium_3d(scn, nrealisations, medium, time_range, fp_to_meter,
+ Tarad, Tref, register_paths, NULL, estimator);
+ }
+}
+
+res_T
+sdis_solve_medium_green_function
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ struct sdis_medium* medium, /* Medium to solve */
+ 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_medium_2d(scn, nrealisations, medium, NULL, fp_to_meter, Tarad,
+ Tref, SDIS_HEAT_PATH_NONE, green, NULL);
+ } else {
+ return solve_medium_3d(scn, nrealisations, medium, NULL, fp_to_meter, Tarad,
+ Tref, SDIS_HEAT_PATH_NONE, green, NULL);
}
- return res;
}
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -15,6 +15,7 @@
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
+#include "sdis_green.h"
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
@@ -205,9 +206,12 @@ XD(solve_medium)
const double Tarad, /* Ambient radiative temperature */
const double Tref, /* Reference temperature */
const int register_paths, /* Combination of enum sdis_heat_path_flag */
- struct sdis_estimator** out_estimator)
+ struct sdis_green_function** out_green, /* May be NULL <=> No green func */
+ struct sdis_estimator** out_estimator) /* May be NULL <=> No estimator */
{
struct darray_enclosure_cumul cumul;
+ 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 sdis_estimator* estimator = NULL;
@@ -218,11 +222,21 @@ XD(solve_medium)
ATOMIC res = RES_OK;
if(!scn || !mdm || !nrealisations || nrealisations > INT64_MAX
- || !time_range || 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[0] > time_range[1]
+ || (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; }
@@ -247,23 +261,37 @@ XD(solve_medium)
accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
if(!accums) { res = RES_MEM_ERR; goto error; }
- /* Create the estimator */
- res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
- if(res != RES_OK) goto error;
-
/* Compute the enclosure cumulative */
darray_enclosure_cumul_init(scn->dev->allocator, &cumul);
cumul_is_init = 1;
res = compute_medium_enclosure_cumulative(scn, mdm, &cumul);
if(res != RES_OK) goto error;
+ 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 */
+ 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)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
const int ithread = omp_get_thread_num();
- const struct enclosure* enc = NULL;
- struct accum* accum = accums + ithread;
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;
+ const struct enclosure* enc = NULL;
struct sdis_heat_path* pheat_path = NULL;
struct sdis_heat_path heat_path;
double weight;
@@ -273,13 +301,23 @@ XD(solve_medium)
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
- /* Sample the time */
- time = sample_time(rng, time_range);
+ if(!out_green) {
+ /* Sample the time */
+ time = sample_time(rng, time_range);
+
+ /* Prepare path registration if necessary */
+ 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; }
- /* Prepare path registration if necessary */
- if(register_paths) {
- heat_path_init(scn->dev->allocator, &heat_path);
- pheat_path = &heat_path;
+ pgreen_path = &green_path;
}
/* Uniformly Sample an enclosure that surround the submitted medium and
@@ -294,7 +332,7 @@ XD(solve_medium)
/* Run a probe realisation */
res_local = XD(probe_realisation)((size_t)irealisation, scn, rng, mdm, pos,
- time, fp_to_meter, Tarad, Tref, NULL, pheat_path, &weight);
+ time, fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &weight);
if(res_local != RES_OK) {
if(res_local != RES_BAD_OP) { ATOMIC_SET(&res, res_local); continue; }
} else {
@@ -320,23 +358,45 @@ XD(solve_medium)
}
if(res != RES_OK) goto error;
- /* Merge the per thread accumulators into the accumulator of the thread 0 */
- sum_accums(accums, scn->dev->nthreads, &accums[0]);
- ASSERT(accums[0].count <= nrealisations);
-
/* 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);
+ }
+
+ if(out_green) {
+ 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;
+ }
+
+ /* 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(cumul_is_init) darray_enclosure_cumul_release(&cumul);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
+ if(out_green) *out_green = green;
return (res_T)res;
error:
if(estimator) {
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -39,7 +39,7 @@ XD(solve_probe)
const double Tref, /* Reference temperature */
const int register_paths, /* Combination of enum sdis_heat_path_flag */
struct sdis_green_function** out_green, /* May be NULL <=> No green func */
- struct sdis_estimator** out_estimator)
+ struct sdis_estimator** out_estimator) /* May be NULL <=> No estimator */
{
struct sdis_medium* medium = NULL;
struct sdis_estimator* estimator = NULL;
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -25,6 +25,7 @@
#define Tf0 300.0
#define Tf1 330.0
#define N 1000ul /* #realisations */
+#define Np 10000ul /* #realisations precise */
/*
* The scene is composed of 2 super shapes whose temperature is unknown. The
@@ -213,6 +214,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 fluid* fluid_param = NULL;
struct solid* solid_param = NULL;
struct interf* interface_param = NULL;
@@ -393,15 +396,33 @@ main(int argc, char** argv)
CHK(eq_eps(v, v0+v1, 1.e-6));
BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
- OK(sdis_solve_medium(scn, 10000, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 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));
ref = Tf0 * v0/v + Tf1 * v1/v;
printf("Shape0 + Shape1 temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
- printf("#failures = %lu/10000\n", nfails);
+ printf("#failures = %lu/%lu\n", nfails, Np);
CHK(eq_eps(T.E, ref, T.SE*3));
+
+ /* Solve green */
+ BA(sdis_solve_medium_green_function(NULL, Np, solid0, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, 0, solid0, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, NULL, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 0.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, -1, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, NULL));
+ BA(sdis_solve_medium_green_function(scn, Np, solid1, 1.0, 0, 0, &green));
+ OK(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, &green));
+
+ OK(sdis_green_function_solve(green, trange, &estimator2));
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
+
+ OK(sdis_green_function_ref_put(green));
+
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Release */
OK(s3dut_mesh_ref_put(msh0));
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -24,6 +24,7 @@
#define Tf0 300.0
#define Tf1 330.0
#define N 1000ul /* #realisations */
+#define Np 10000ul /* #realisations precise */
/*
* The scene is composed of a square and a disk whose temperature is unknown.
@@ -198,6 +199,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 fluid* fluid_param = NULL;
struct solid* solid_param = NULL;
struct interf* interface_param = NULL;
@@ -364,16 +367,34 @@ main(int argc, char** argv)
/* Estimate the temperature of the square and disk shapes */
BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
- OK(sdis_solve_medium(scn, 10000, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 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));
ref = Tf0 * a0/a + Tf1 * a1/a;
printf("Square + Disk temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
- printf("#failures = %lu / 10000\n", nfails);
+ printf("#failures = %lu / %lu\n", nfails, Np);
CHK(eq_eps(T.E, ref, 3*T.SE));
- CHK(nreals + nfails == 10000);
+ CHK(nreals + nfails == Np);
+
+ /* Solve green */
+ BA(sdis_solve_medium_green_function(NULL, Np, solid0, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, 0, solid0, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, NULL, 1.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 0.0, 0, 0, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, -1, &green));
+ BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, NULL));
+ BA(sdis_solve_medium_green_function(scn, Np, solid1, 1.0, 0, 0, &green));
+ OK(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, &green));
+
+ OK(sdis_green_function_solve(green, trange, &estimator2));
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
+
+ OK(sdis_green_function_ref_put(green));
+
OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
/* Release */
OK(sdis_device_ref_put(dev));
