commit 770227f1326854bb4f4d0c07348ac022921a3998 parent f92eddb90dbe15af14c4b66237397f92470496f4 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Thu, 21 Feb 2019 16:48:12 +0100 First implementation of the heat path registration Diffstat:
32 files changed, 368 insertions(+), 111 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h @@ -51,12 +51,11 @@ struct senc2d_descriptor; struct senc_descriptor; /* Forward declaration of the Stardis opaque data types. These data types are - * ref counted. Once created with the appropriated `sdis_<TYPE>_create' - * function, the caller implicitly owns the created data, i.e. its reference - * counter is set to 1. The sdis_<TYPE>_ref_<get|put> functions get or release - * a reference on the data, i.e. they increment or decrement the reference - * counter, respectively. When this counter reaches 0, the object is silently - * destroyed and cannot be used anymore. */ + * ref counted. Once created the caller implicitly owns the created data, i.e. + * its reference counter is set to 1. The sdis_<TYPE>_ref_<get|put> functions + * get or release a reference on the data, i.e. they increment or decrement the + * reference counter, respectively. When this counter reaches 0, the object is + * silently destroyed and cannot be used anymore. */ struct sdis_accum_buffer; struct sdis_camera; struct sdis_data; @@ -97,6 +96,19 @@ enum sdis_point_type { SDIS_POINT_NONE = SDIS_POINT_TYPES_COUNT__ }; +enum sdis_heat_vertex_type { + SDIS_HEAT_VERTEX_CONDUCTION, + SDIS_HEAT_VERTEX_CONVECTION, + SDIS_HEAT_VERTEX_RADIATIVE +}; + +enum sdis_heat_path_flag { + SDIS_HEAT_PATH_OK = BIT(0), + SDIS_HEAT_PATH_FAILED = BIT(1), + SDIS_HEAT_PATH_ALL = SDIS_HEAT_PATH_OK | SDIS_HEAT_PATH_FAILED, + SDIS_HEAT_PATH_NONE = 0 +}; + /* Random walk vertex, i.e. a spatiotemporal position at a given step of the * random walk. */ struct sdis_rwalk_vertex { @@ -245,6 +257,18 @@ typedef res_T const size_t naccums[2], /* #accumulations in X and Y */ const struct sdis_accum* accums); /* List of row ordered accumulations */ +/* Vertex of heat path v*/ +struct sdis_heat_vertex { + double P[3]; + double time; + double weight; + enum sdis_heat_vertex_type type; +}; +#define SDIS_HEAT_VERTEX_NULL__ {{0,0,0}, 0, 0, SDIS_HEAT_VERTEX_CONDUCTION} +static const struct sdis_heat_vertex SDIS_HEAT_VERTEX_NULL = + SDIS_HEAT_VERTEX_NULL__; + +/* Path used to estimate the green function */ struct sdis_green_path { /* Internal data. Should not be accessed */ void* green__; @@ -765,6 +789,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 */ struct sdis_estimator** estimator); SDIS_API res_T diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h @@ -20,13 +20,15 @@ /* Forward declaration */ struct green_path_handle; +struct heat_path; struct rwalk_context { struct green_path_handle* green_path; + struct heat_path* heat_path; double Tarad; /* Ambient radiative temperature */ double Tref3; /* Reference temperature ^ 3 */ }; -#define RWALK_CONTEXT_NULL__ {NULL, 0, 0} +#define RWALK_CONTEXT_NULL__ {NULL, NULL, 0, 0} static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__; #endif /* SDIS_XD_BEGIN_H */ diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -17,6 +17,8 @@ #include "sdis_device_c.h" #include "sdis_estimator_c.h" +#include <rsys/mutex.h> + /******************************************************************************* * Helper functions ******************************************************************************/ @@ -28,6 +30,8 @@ estimator_release(ref_T* ref) ASSERT(ref); estimator = CONTAINER_OF(ref, struct sdis_estimator, ref); dev = estimator->dev; + darray_heat_path_release(&estimator->paths); + if(estimator->mutex) mutex_destroy(estimator->mutex); MEM_RM(dev->allocator, estimator); SDIS(device_ref_put(dev)); } @@ -127,8 +131,6 @@ res_T estimator_create (struct sdis_device* dev, const enum sdis_estimator_type type, - const size_t nrealisations, - const size_t nsuccesses, struct sdis_estimator** out_estimator) { struct sdis_estimator* estimator = NULL; @@ -136,9 +138,6 @@ estimator_create if(!dev || (unsigned)type >= SDIS_ESTIMATOR_TYPES_COUNT__ - || !nrealisations - || !nsuccesses - || nsuccesses > nrealisations || !out_estimator) { res = RES_BAD_ARG; goto error; @@ -151,10 +150,18 @@ estimator_create } ref_init(&estimator->ref); SDIS(device_ref_get(dev)); - estimator->nrealisations = nsuccesses; - estimator->nfailures = nrealisations - nsuccesses; + estimator->nrealisations = 0; + estimator->nfailures = 0; estimator->dev = dev; estimator->type = type; + darray_heat_path_init(dev->allocator, &estimator->paths); + + estimator->mutex = mutex_create(); + if(!estimator->mutex) { + res = RES_MEM_ERR; + goto error; + } + exit: if(out_estimator) *out_estimator = estimator; return res; @@ -166,3 +173,30 @@ error: goto exit; } +res_T +estimator_add_and_release_heat_path + (struct sdis_estimator* estimator, struct heat_path* path) +{ + struct heat_path* dst = NULL; + size_t i; + res_T res = RES_OK; + ASSERT(estimator && path); + + mutex_lock(estimator->mutex); + + i = darray_heat_path_size_get(&estimator->paths); + + res = darray_heat_path_resize(&estimator->paths, i+1); + if(res != RES_OK) goto error; + + dst = darray_heat_path_data_get(&estimator->paths) + i; + res = heat_path_copy_and_release(dst, path); + if(res != RES_OK) goto error; + +exit: + mutex_unlock(estimator->mutex); + return res; +error: + goto exit; +} + diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -16,6 +16,8 @@ #ifndef SDIS_ESTIMATOR_C_H #define SDIS_ESTIMATOR_C_H +#include "sdis_heat_path.h" + #include <rsys/math.h> #include <rsys/ref_count.h> @@ -37,11 +39,16 @@ struct sdis_estimator { size_t nrealisations; size_t nfailures; + struct mutex* mutex; + struct darray_heat_path paths; /* Tracked paths */ + enum sdis_estimator_type type; ref_T ref; struct sdis_device* dev; }; +struct sdis_estimator_handle; + /******************************************************************************* * Estimator local API ******************************************************************************/ @@ -49,10 +56,26 @@ extern LOCAL_SYM res_T estimator_create (struct sdis_device* dev, const enum sdis_estimator_type type, - const size_t nrealisations, - const size_t nsuccesses, struct sdis_estimator** estimator); +/* Thread safe */ +extern LOCAL_SYM res_T +estimator_add_and_release_heat_path + (struct sdis_estimator* estimator, + struct heat_path* path); + +/* Must be invoked before any others "estimator_setup" functions */ +static INLINE void +estimator_setup_realisations_count + (struct sdis_estimator* estimator, + const size_t nrealisations, + const size_t nsuccesses) +{ + ASSERT(estimator && nrealisations && nsuccesses && nsuccesses<=nrealisations); + estimator->nrealisations = nsuccesses; + estimator->nfailures = nrealisations - nsuccesses; +} + static INLINE void estimator_setup_temperature (struct sdis_estimator* estim, diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -459,6 +459,10 @@ sdis_green_function_solve npaths = darray_green_path_size_get(&green->paths); + /* Create the estimator */ + res = estimator_create(green->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator); + if(res != RES_OK) goto error; + /* Solve the green function */ FOR_EACH(ipath, 0, npaths) { /* TODO add multi-threading (?) */ const double time = sample_time(rng, time_range); @@ -473,12 +477,8 @@ sdis_green_function_solve ++N; } - /* Create the estimator */ - res = estimator_create - (green->dev, SDIS_ESTIMATOR_TEMPERATURE, npaths, N, &estimator); - if(res != RES_OK) goto error; - /* Setup the estimated temperature */ + estimator_setup_realisations_count(estimator, npaths, N); estimator_setup_temperature(estimator, accum, accum2); exit: diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -16,8 +16,12 @@ #ifndef SDIS_HEAT_PATH_H #define SDIS_HEAT_PATH_H +#include "sdis.h" + +#include <rsys/dynamic_array.h> #include <rsys/rsys.h> +/* Forward declarations */ struct rwalk_2d; struct rwalk_3d; struct rwalk_context; @@ -26,6 +30,74 @@ struct ssp_rng; struct temperature_2d; struct temperature_3d; +/* Generate the dynamic array of heat vertices */ +#define DARRAY_NAME heat_vertex +#define DARRAY_DATA struct sdis_heat_vertex +#include <rsys/dynamic_array.h> + +/******************************************************************************* + * Heat path data structure + ******************************************************************************/ +struct heat_path { + struct darray_heat_vertex vertices; + enum sdis_heat_path_flag status; +}; + +static INLINE void +heat_path_init(struct mem_allocator* allocator, struct heat_path* path) +{ + ASSERT(path); + path->status = SDIS_HEAT_PATH_NONE; + darray_heat_vertex_init(allocator, &path->vertices); +} + +static INLINE void +heat_path_release(struct heat_path* path) +{ + ASSERT(path); + darray_heat_vertex_release(&path->vertices); +} + +static INLINE res_T +heat_path_copy(struct heat_path* dst, const struct heat_path* src) +{ + ASSERT(dst && src); + dst->status = src->status; + return darray_heat_vertex_copy(&dst->vertices, &src->vertices); +} + +static INLINE res_T +heat_path_copy_and_release(struct heat_path* dst, struct heat_path* src) +{ + ASSERT(dst && src); + dst->status = src->status; + return darray_heat_vertex_copy_and_release(&dst->vertices, &src->vertices); +} + +static INLINE res_T +heat_path_copy_and_clear(struct heat_path* dst, struct heat_path* src) +{ + ASSERT(dst && src); + dst->status = src->status; + return darray_heat_vertex_copy_and_clear(&dst->vertices, &src->vertices); +} + +static INLINE res_T +heat_path_add_vertex(struct heat_path* path, const struct sdis_heat_vertex* vtx) +{ + ASSERT(path && vtx); + return darray_heat_vertex_push_back(&path->vertices, vtx); +} + +/* Generate the dynamic array of heat paths */ +#define DARRAY_NAME heat_path +#define DARRAY_DATA struct heat_path +#define DARRAY_FUNCTOR_INIT heat_path_init +#define DARRAY_FUNCTOR_RELEASE heat_path_release +#define DARRAY_FUNCTOR_COPY heat_path_copy +#define DARRAY_FUNCTOR_COPY_AND_RELEASE heat_path_copy_and_release +#include <rsys/dynamic_array.h> + /******************************************************************************* * Trace or pursue a radiative path ******************************************************************************/ diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -252,6 +252,11 @@ XD(solid_solid_boundary_path) rwalk->hit_side = SDIS_SIDE_NULL__; } + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + if(res != RES_OK) goto error; + exit: return res; error: @@ -411,6 +416,11 @@ XD(solid_fluid_boundary_path) rwalk->hit = SXD_HIT_NULL; rwalk->hit_side = SDIS_SIDE_NULL__; } + + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + if(res != RES_OK) goto error; } exit: @@ -541,8 +551,14 @@ XD(solid_boundary_with_flux_path) rwalk->mdm = mdm; rwalk->hit = SXD_HIT_NULL; rwalk->hit_side = SDIS_SIDE_NULL__; + } + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + if(res != RES_OK) goto error; + exit: return res; error: diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -176,7 +176,8 @@ XD(conductive_path) } } - /* Register the power term for the green function */ + /* Register the power term for the green function. Delay its registration + * until the end of the conductive path, i.e. the path is valid */ if(ctx->green_path && power != SDIS_VOLUMIC_POWER_NONE) { green_power_factor += power_factor; } @@ -194,6 +195,12 @@ XD(conductive_path) if(tmp >= 0) { T->value += tmp; T->done = 1; + + /* Register the initial vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + if(res != RES_OK) goto error; + break; } /* The initial condition should have been reached */ @@ -218,6 +225,11 @@ XD(conductive_path) /* Update the random walk position */ XD(move_pos)(rwalk->vtx.P, dir0, delta); + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION); + if(res != RES_OK) goto error; + /* Fetch the current medium */ if(SXD_HIT_NONE(&rwalk->hit)) { CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK); diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -33,7 +33,6 @@ XD(convective_path) struct XD(temperature)* T) { struct sXd(attrib) attr_P, attr_N; - struct sdis_interface_fragment frag; const struct sdis_interface* interf; const struct enclosure* enc; unsigned enc_ids[2]; @@ -147,30 +146,10 @@ XD(convective_path) goto error; } - /* A trick to force first r test result. - * TODO fix this workaround that seems useless */ - r = 1; - /* Sample time until init condition is reached or a true convection occurs. */ for(;;) { + struct sdis_interface_fragment frag; struct sXd(primitive) prim; - /* Setup the fragment of the interface. */ - XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side); - - /* Fetch hc. */ - hc = interface_get_convection_coef(interf, &frag); - if(hc > enc->hc_upper_bound) { - log_err(scn->dev, - "%s: hc (%g) exceeds its provided upper bound (%g) at %g %g %g.\n", - FUNC_NAME, hc, enc->hc_upper_bound, SPLIT3(rwalk->vtx.P)); - res = RES_BAD_OP; - goto error; - } - - if(r < hc / enc->hc_upper_bound) { - /* True convection. Always true if hc == bound. */ - break; - } /* Fetch other physical properties. */ cp = fluid_get_calorific_capacity(rwalk->mdm, &rwalk->vtx); @@ -183,6 +162,12 @@ XD(convective_path) tau = ssp_ran_exp(rng, mu); t0 = ctx->green_path ? -INF : fluid_get_t0(rwalk->mdm); rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); + + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONVECTION); + if(res != RES_OK) goto error; + if(rwalk->vtx.time == t0) { /* Check the initial condition. */ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx); @@ -238,8 +223,29 @@ XD(convective_path) FATAL("Unexpected fluid interface.\n"); } - /* Renew r for next loop. */ + /* Register the new vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONVECTION); + if(res != RES_OK) goto error; + + /* Setup the fragment of the sampled position into the enclosure. */ + XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side); + + /* Fetch the convection coefficient of the sampled position */ + hc = interface_get_convection_coef(interf, &frag); + if(hc > enc->hc_upper_bound) { + log_err(scn->dev, + "%s: hc (%g) exceeds its provided upper bound (%g) at %g %g %g.\n", + FUNC_NAME, hc, enc->hc_upper_bound, SPLIT3(rwalk->vtx.P)); + res = RES_BAD_OP; + goto error; + } + r = ssp_rng_canonical_float(rng); + if(r < hc / enc->hc_upper_bound) { + /* True convection. Always true if hc == bound. */ + break; + } } rwalk->hit.distance = 0; diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -107,6 +107,11 @@ XD(trace_radiative_path) /* Move the random walk to the hit position */ XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance); + /* Register the random walk vertex against the heat path */ + res = register_heat_vertex + (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_RADIATIVE); + if(res != RES_OK) goto error; + /* Fetch the new interface and setup the hit fragment */ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side); @@ -179,7 +184,7 @@ XD(radiative_path) struct XD(temperature)* T) { /* The radiative random walk is always performed in 3D. In 2D, the geometry - * are assumed to be extruded to the infinty along the Z dimension. */ + * are assumed to be extruded to the infinity along the Z dimension. */ float N[3] = {0, 0, 0}; float dir[3] = {0, 0, 0}; res_T res = RES_OK; diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -87,4 +87,25 @@ sample_time(struct ssp_rng* rng, const double time_range[2]) return ssp_rng_uniform_double(rng, time_range[0], time_range[1]); } +static INLINE res_T +register_heat_vertex + (struct heat_path* path, + const struct sdis_rwalk_vertex* vtx, + const double weight, + const enum sdis_heat_vertex_type type) +{ + struct sdis_heat_vertex heat_vtx = SDIS_HEAT_VERTEX_NULL; + ASSERT(vtx); + + if(!path) return RES_OK; + + heat_vtx.P[0] = vtx->P[0]; + heat_vtx.P[1] = vtx->P[1]; + heat_vtx.P[2] = vtx->P[2]; + heat_vtx.time = vtx->time; + heat_vtx.weight = weight; + heat_vtx.type = type; + return heat_path_add_vertex(path, &heat_vtx); +} + #endif /* SDIS_MISC_H */ diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h @@ -46,6 +46,7 @@ probe_realisation_2d const double ambient_radiative_temperature, const double reference_temperature, struct green_path_handle* green_path, + struct heat_path* heat_path, double* weight); extern LOCAL_SYM res_T @@ -59,6 +60,7 @@ probe_realisation_3d const double ambient_radiative_temperature, const double reference_temperature, struct green_path_handle* green_path, + struct heat_path* heat_path, double* weight); /******************************************************************************* diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -97,11 +97,13 @@ XD(probe_realisation) const double ambient_radiative_temperature, const double reference_temperature, struct green_path_handle* green_path, + struct heat_path* heat_path, double* weight) { struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct XD(rwalk) rwalk = XD(RWALK_NULL); struct XD(temperature) T = XD(TEMPERATURE_NULL); + enum sdis_heat_vertex_type type; double t0; double (*get_initial_temperature) (const struct sdis_medium* mdm, @@ -125,6 +127,14 @@ XD(probe_realisation) dX(set)(rwalk.vtx.P, position); rwalk.vtx.time = time; + + /* Register the starting position against the heat path */ + type = medium->type == SDIS_SOLID + ? SDIS_HEAT_VERTEX_CONDUCTION + : SDIS_HEAT_VERTEX_CONVECTION; + res = register_heat_vertex(heat_path, &rwalk.vtx, 0, type); + if(res != RES_OK) goto error; + /* No initial condition with green */ if(!green_path && t0 >= rwalk.vtx.time) { double tmp; @@ -133,20 +143,22 @@ XD(probe_realisation) tmp = get_initial_temperature(medium, &rwalk.vtx); if(tmp >= 0) { *weight = tmp; - return RES_OK; + goto exit; } /* The initial condition should have been reached */ log_err(scn->dev, "%s: undefined initial condition. " "The time is %f but the temperature remains unknown.\n", FUNC_NAME, t0); - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } rwalk.hit = SXD_HIT_NULL; rwalk.mdm = medium; ctx.green_path = green_path; + ctx.heat_path = heat_path; ctx.Tarad = ambient_radiative_temperature; ctx.Tref3 = reference_temperature @@ -154,10 +166,14 @@ XD(probe_realisation) * reference_temperature; 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 @@ -167,15 +167,16 @@ sdis_solve_probe const double fp_to_meter,/* Scale factor from floating point unit to meter */ 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) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { return solve_probe_2d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, NULL, out_estimator); + fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); } else { return solve_probe_3d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, NULL, out_estimator); + fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); } } @@ -192,10 +193,10 @@ sdis_solve_probe_green_function if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { return solve_probe_2d(scn, nrealisations, position, NULL, - fp_to_meter, Tarad, Tref, out_green, NULL); + fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, out_green, NULL); } else { return solve_probe_3d(scn, nrealisations, position, NULL, - fp_to_meter, Tarad, Tref, out_green, NULL); + fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, out_green, NULL); } } diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -157,6 +157,7 @@ XD(solve_probe) const double fp_to_meter,/* Scale factor from floating point unit to meter */ const double Tarad, /* Ambient radiative temperature */ 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) { @@ -225,6 +226,12 @@ XD(solve_probe) } + /* Create the estimator */ + 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); #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N) @@ -236,11 +243,17 @@ XD(solve_probe) struct ssp_rng* rng = rngs[ithread]; struct green_path_handle* pgreen_path = NULL; struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL; + struct heat_path* pheat_path = NULL; + struct heat_path heat_path; - if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */ + if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ 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 */ @@ -252,7 +265,7 @@ XD(solve_probe) } res_local = XD(probe_realisation)(scn, rng, medium, position, time, - fp_to_meter, Tarad, Tref, pgreen_path, &w); + fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w); if(res_local != RES_OK) { if(res_local != RES_BAD_OP) { ATOMIC_SET(&res, res_local); continue; } } else { @@ -260,16 +273,26 @@ XD(solve_probe) sqr_weight += w*w; ++N; } + + if(pheat_path) { + pheat_path->status = res_local == RES_OK + ? SDIS_HEAT_PATH_OK + : 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; + /* Setup the estimated temperature */ if(out_estimator) { - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_TEMPERATURE, nrealisations, N, &estimator); - if(res != RES_OK) goto error; - - /* Setup the estimated temperature */ + estimator_setup_realisations_count(estimator, nrealisations, N); estimator_setup_temperature(estimator, weight, sqr_weight); } @@ -407,6 +430,10 @@ XD(solve_probe_boundary) if(res != RES_OK) goto error; } + /* Create the 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); #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N) @@ -436,12 +463,8 @@ XD(solve_probe_boundary) } if(res != RES_OK) goto error; - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_TEMPERATURE, nrealisations, N, &estimator); - if(res != RES_OK) goto error; - /* Setup the estimated temperature */ + estimator_setup_realisations_count(estimator, nrealisations, N); estimator_setup_temperature(estimator, weight, sqr_weight); exit: @@ -564,6 +587,10 @@ XD(solve_boundary) if(res != RES_OK) goto error; } + /* Create the 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) reduction(+:weight,sqr_weight,N) for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) { @@ -621,12 +648,8 @@ XD(solve_boundary) } } - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_TEMPERATURE, nrealisations, N, &estimator); - if(res != RES_OK) goto error; - /* Setup the estimated temperature */ + estimator_setup_realisations_count(estimator, nrealisations, N); estimator_setup_temperature(estimator, weight, sqr_weight); exit: @@ -761,6 +784,10 @@ XD(solve_probe_boundary_flux) res = XD(interface_prebuild_fragment) (&frag, scn, (unsigned)iprim, uv, fluid_side); + /* Create the estimator */ + res = estimator_create(scn->dev, SDIS_ESTIMATOR_FLUX, &estimator); + if(res != RES_OK) goto error; + /* Here we go! Launch the Monte Carlo estimation */ omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) reduction(+:weight_t,sqr_weight_t,\ @@ -814,12 +841,8 @@ XD(solve_probe_boundary_flux) } if(res != RES_OK) goto error; - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_FLUX, nrealisations, N, &estimator); - if(res != RES_OK) goto error; - /* Setup the estimated values */ + estimator_setup_realisations_count(estimator, nrealisations, N); estimator_setup_temperature(estimator, weight_t, sqr_weight_t); estimator_setup_flux(estimator, FLUX_CONVECTIVE, weight_fc, sqr_weight_fc); estimator_setup_flux(estimator, FLUX_RADIATIVE, weight_fr, sqr_weight_fr); @@ -947,6 +970,10 @@ XD(solve_boundary_flux) if(res != RES_OK) goto error; } + /* Create the estimator */ + res = estimator_create(scn->dev, SDIS_ESTIMATOR_FLUX, &estimator); + if(res != RES_OK) goto error; + omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) reduction(+:weight_t,sqr_weight_t,\ weight_fc,sqr_weight_fc,weight_fr,sqr_weight_fr,weight_f,sqr_weight_f,N) @@ -1045,12 +1072,8 @@ XD(solve_boundary_flux) } if(res != RES_OK) goto error; - /* Create the estimator */ - res = estimator_create - (scn->dev, SDIS_ESTIMATOR_FLUX, nrealisations, N, &estimator); - if(res != RES_OK) goto error; - /* Setup the estimated values */ + estimator_setup_realisations_count(estimator, nrealisations, N); estimator_setup_temperature(estimator, weight_t, sqr_weight_t); estimator_setup_flux(estimator, FLUX_CONVECTIVE, weight_fc, sqr_weight_fc); estimator_setup_flux(estimator, FLUX_RADIATIVE, weight_fr, sqr_weight_fr); diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -398,7 +398,7 @@ main(int argc, char** argv) pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); pos[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); - OK(sdis_solve_probe(scn, N, pos, time_range, 1, -1, Tref, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1, -1, Tref, 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)); diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -403,7 +403,7 @@ main(int argc, char** argv) pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); - OK(sdis_solve_probe(scn, 10000, pos, time_range, 1, -1, Tref, &estimator)); + OK(sdis_solve_probe(scn, 10000, pos, time_range, 1, -1, Tref, 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)); diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -277,7 +277,7 @@ main(int argc, char** argv) *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 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)); @@ -313,7 +313,7 @@ main(int argc, char** argv) *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 2D */ - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -293,7 +293,7 @@ main(int argc, char** argv) *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -327,7 +327,7 @@ main(int argc, char** argv) *((int*)sdis_data_get(is_stationary)) = IS_INF(time); - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -145,7 +145,7 @@ solve(struct sdis_scene* scn, const double 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)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -263,13 +263,13 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - BA(sdis_solve_probe(NULL, N, pos, time_range, 1.0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL)); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(NULL, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, NULL)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); BA(sdis_estimator_get_type(estimator, NULL)); BA(sdis_estimator_get_type(NULL, &type)); @@ -318,10 +318,10 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); fluid_param->temperature = 300; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); BA(sdis_solve_probe_green_function(NULL, N, pos, 1.0, 0, 0, &green)); BA(sdis_solve_probe_green_function(scn, 0, pos, 1.0, 0, 0, &green)); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -242,7 +242,7 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 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)); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -240,7 +240,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 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)); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -298,7 +298,7 @@ main(int argc, char** argv) pos[1] = 0.5; pos[2] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, &estimator)); + OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 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)); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -291,7 +291,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 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)); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -201,7 +201,7 @@ main(int argc, char** argv) pos[0] = 0.5; pos[1] = 0.5; time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); @@ -228,7 +228,7 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -278,7 +278,7 @@ static INLINE void check_memory_allocator(struct mem_allocator* allocator) { if(MEM_ALLOCATED_SIZE(allocator)) { - char dump[128]; + char dump[1024]; MEM_DUMP(allocator, dump, sizeof(dump)); fprintf(stderr, "%s\n", dump); FATAL("Memory leaks.\n"); diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -169,7 +169,7 @@ solve(struct sdis_scene* scn, const double pos[]) 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)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -242,7 +242,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) pos[1] = refs[i].pos[1]; pos[2] = refs[i].pos[2]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 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)); @@ -256,7 +256,6 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) } } - int main(int argc, char** argv) { diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -261,7 +261,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) pos[0] = refs[i].pos[0]; pos[1] = refs[i].pos[1]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 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)); diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -441,7 +441,7 @@ main(int argc, char** argv) FATAL("Unreachable code.\n"); } - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c @@ -348,7 +348,7 @@ main(int argc, char** argv) Tref = T2 + (T1-T2)/L * (pos[1] + vertices[3]); #endif - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, &estimator)); + OK(sdis_solve_probe(scn, N, pos, time_range, 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));