commit 7f389490ab19caf4e600e49715880095a5b2d23d parent 2c6465e3e59f6fe15fccb7fc8d9843de73c9a44d Author: Vincent Forest <vincent.forest@meso-star.com> Date: Mon, 11 Feb 2019 15:58:49 +0100 Add green function registration to sdis_solve_probe Diffstat:
34 files changed, 416 insertions(+), 147 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h @@ -571,13 +571,23 @@ sdis_estimator_get_radiative_flux (const struct sdis_estimator* estimator, struct sdis_mc* flux); - SDIS_API res_T sdis_estimator_get_total_flux (const struct sdis_estimator* estimator, struct sdis_mc* flux); /******************************************************************************* + * The green function saves the estimation of the propagator + ******************************************************************************/ +SDIS_API res_T +sdis_green_function_ref_get + (struct sdis_green_function* green); + +SDIS_API res_T +sdis_green_function_ref_put + (struct sdis_green_function* green); + +/******************************************************************************* * Miscellaneous functions ******************************************************************************/ SDIS_API res_T @@ -589,6 +599,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 */ + struct sdis_green_function** green, /* NULL <=> no green registration */ struct sdis_estimator** estimator); SDIS_API res_T diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h @@ -18,10 +18,16 @@ #include <rsys/rsys.h> +/* Forward declaration */ +struct green_path_handle; + struct rwalk_context { + struct green_path_handle* green_path; double Tarad; /* Ambient radiative temperature */ double Tref3; /* Reference temperature ^ 3 */ }; +#define RWALK_CONTEXT_NULL__ {NULL, 0, 0} +static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__; #endif /* SDIS_XD_BEGIN_H */ @@ -80,7 +86,7 @@ struct rwalk_context { /* Current state of the random walk */ struct XD(rwalk) { struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */ - const struct sdis_medium* mdm; /* Medium in which the random walk lies */ + struct sdis_medium* mdm; /* Medium in which the random walk lies */ struct sXd(hit) hit; /* Hit of the random walk */ enum sdis_side hit_side; }; diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -105,6 +105,24 @@ green_path_copy(struct green_path* dst, const struct green_path* src) } static INLINE res_T +green_path_copy_and_clear(struct green_path* dst, struct green_path* src) +{ + res_T res = RES_OK; + ASSERT(dst && src); + + dst->limit_vertex = src->limit_vertex; + dst->end_on_interface = src->end_on_interface; + dst->ilast_medium = src->ilast_medium; + dst->ilast_interf = src->ilast_interf; + res = darray_green_term_copy_and_clear(&dst->flux_terms, &src->flux_terms); + if(res != RES_OK) return res; + res = darray_green_term_copy_and_clear(&dst->power_terms, &src->power_terms); + if(res != RES_OK) return res; + return RES_OK; + +} + +static INLINE res_T green_path_copy_and_release(struct green_path* dst, struct green_path* src) { res_T res = RES_OK; @@ -247,6 +265,25 @@ green_function_release(ref_T* ref) } /******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +sdis_green_function_ref_get(struct sdis_green_function* green) +{ + if(!green) return RES_BAD_ARG; + ref_get(&green->ref); + return RES_OK; +} + +res_T +sdis_green_function_ref_put(struct sdis_green_function* green) +{ + if(!green) return RES_BAD_ARG; + ref_put(&green->ref, green_function_release); + return RES_OK; +} + +/******************************************************************************* * Local functions ******************************************************************************/ res_T @@ -274,24 +311,70 @@ exit: return res; error: if(green) { - green_function_ref_put(green); + SDIS(green_function_ref_put(green)); green = NULL; } goto exit; } -void -green_function_ref_get(struct sdis_green_function* green) +res_T +green_function_merge_and_clear + (struct sdis_green_function* dst, struct sdis_green_function* src) { - ASSERT(green); - ref_get(&green->ref); -} + struct htable_medium_iterator it_medium, end_medium; + struct htable_interf_iterator it_interf, end_interf; + struct green_path* paths_src; + struct green_path* paths_dst; + size_t npaths_src; + size_t npaths_dst; + size_t npaths; + size_t i; + unsigned id; + res_T res = RES_OK; + ASSERT(dst && src); -void -green_function_ref_put(struct sdis_green_function* green) -{ - ASSERT(green); - ref_put(&green->ref, green_function_release); + if(dst == src) goto exit; + + npaths_src = darray_green_path_size_get(&src->paths); + npaths_dst = darray_green_path_size_get(&dst->paths); + npaths = npaths_src + npaths_dst; + + res = darray_green_path_resize(&dst->paths, npaths); + if(res != RES_OK) goto error; + + paths_src = darray_green_path_data_get(&src->paths); + paths_dst = darray_green_path_data_get(&dst->paths) + npaths_dst; + + FOR_EACH(i, 0, darray_green_path_size_get(&src->paths)) { + res = green_path_copy_and_clear(&paths_dst[i], &paths_src[i]); + if(res != RES_OK) goto error; + } + + htable_medium_begin(&src->media, &it_medium); + htable_medium_end(&src->media, &end_medium); + while(!htable_medium_iterator_eq(&it_medium, &end_medium)) { + struct sdis_medium* medium; + medium = *htable_medium_iterator_data_get(&it_medium); + id = medium_get_id(medium); + res = htable_medium_set(&dst->media, &id, &medium); + if(res != RES_OK) goto error; + } + + htable_interf_begin(&src->interfaces, &it_interf); + htable_interf_end(&src->interfaces, &end_interf); + while(!htable_interf_iterator_eq(&it_interf, &end_interf)) { + struct sdis_interface* interf; + interf = *htable_interf_iterator_data_get(&it_interf); + id = interface_get_id(interf); + res = htable_interf_set(&dst->interfaces, &id, &interf); + if(res != RES_OK) goto error; + } + + green_function_clear(src); +exit: + return res; +error: + goto exit; } res_T diff --git a/src/sdis_green.h b/src/sdis_green.h @@ -26,19 +26,20 @@ struct green_path_handle { struct sdis_green_function* green; struct green_path* path; }; +#define GREEN_PATH_HANDLE_NULL__ {NULL, NULL} +static const struct green_path_handle GREEN_PATH_HANDLE_NULL = + GREEN_PATH_HANDLE_NULL__; extern LOCAL_SYM res_T green_function_create (struct sdis_device* dev, struct sdis_green_function** green); -extern LOCAL_SYM void -green_function_ref_get - (struct sdis_green_function* greeN); - -extern LOCAL_SYM void -green_function_ref_put - (struct sdis_green_function* green); +/* Merge `src' into `dst' an clear `src' */ +extern LOCAL_SYM res_T +green_function_merge_and_clear + (struct sdis_green_function* dst, + struct sdis_green_function* src); extern LOCAL_SYM res_T green_function_create_path diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -14,6 +14,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "sdis_device_c.h" +#include "sdis_green.h" #include "sdis_heat_path.h" #include "sdis_interface_c.h" #include "sdis_medium_c.h" @@ -100,7 +101,7 @@ XD(check_rwalk_fragment_consistency) return d2_eq_eps(uv, frag->uv, 1.e-6); } -static void +static res_T XD(solid_solid_boundary_path) (const struct sdis_scene* scn, const double fp_to_meter, @@ -112,10 +113,10 @@ XD(solid_solid_boundary_path) { struct sXd(hit) hit0, hit1, hit2, hit3; struct sXd(hit)* hit; - const struct sdis_interface* interf = NULL; - const struct sdis_medium* solid_front = NULL; - const struct sdis_medium* solid_back = NULL; - const struct sdis_medium* mdm; + struct sdis_interface* interf = NULL; + struct sdis_medium* solid_front = NULL; + struct sdis_medium* solid_back = NULL; + struct sdis_medium* mdm; double lambda_front, lambda_back; double delta_front, delta_back; double delta_boundary_front, delta_boundary_back; @@ -131,6 +132,7 @@ XD(solid_solid_boundary_path) float* dir; float pos[DIM]; int dim = DIM; + res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); (void)frag, (void)ctx; @@ -227,8 +229,13 @@ XD(solid_solid_boundary_path) if(power != SDIS_VOLUMIC_POWER_NONE) { const double delta_in_meter = reinject_dst * fp_to_meter; const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx); - tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda); - T->value += tmp; + tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda); + T->value += power * tmp; + + if(ctx->green_path) { + res = green_path_add_power_term(ctx->green_path, mdm, tmp); + if(res != RES_OK) goto error; + } } /* Reinject */ @@ -244,9 +251,14 @@ XD(solid_solid_boundary_path) rwalk->hit = SXD_HIT_NULL; rwalk->hit_side = SDIS_SIDE_NULL__; } + +exit: + return res; +error: + goto exit; } -static void +static res_T XD(solid_fluid_boundary_path) (const struct sdis_scene* scn, const double fp_to_meter, @@ -256,11 +268,11 @@ XD(solid_fluid_boundary_path) struct ssp_rng* rng, struct XD(temperature)* T) { - const struct sdis_interface* interf = NULL; - const struct sdis_medium* mdm_front = NULL; - const struct sdis_medium* mdm_back = NULL; - const struct sdis_medium* solid = NULL; - const struct sdis_medium* fluid = NULL; + struct sdis_interface* interf = NULL; + struct sdis_medium* mdm_front = NULL; + struct sdis_medium* mdm_back = NULL; + struct sdis_medium* solid = NULL; + struct sdis_medium* fluid = NULL; struct sXd(hit) hit0 = SXD_HIT_NULL; struct sXd(hit) hit1 = SXD_HIT_NULL; struct sdis_interface_fragment frag_fluid; @@ -278,7 +290,7 @@ XD(solid_fluid_boundary_path) float dir0[DIM], dir1[DIM]; float range[2]; int dim = DIM; - + res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); @@ -377,8 +389,13 @@ XD(solid_fluid_boundary_path) const double power = solid_get_volumic_power(solid, &rwalk->vtx); if(power != SDIS_VOLUMIC_POWER_NONE) { const double delta_in_meter = delta_boundary * fp_to_meter; - tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda); - T->value += tmp; + tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda); + T->value += power * tmp; + + if(ctx->green_path) { + res = green_path_add_power_term(ctx->green_path, solid, tmp); + if(res != RES_OK) goto error; + } } /* Reinject */ @@ -395,9 +412,14 @@ XD(solid_fluid_boundary_path) rwalk->hit_side = SDIS_SIDE_NULL__; } } + +exit: + return res; +error: + goto exit; } -static void +static res_T XD(solid_boundary_with_flux_path) (const struct sdis_scene* scn, const double fp_to_meter, @@ -408,8 +430,8 @@ XD(solid_boundary_with_flux_path) struct ssp_rng* rng, struct XD(temperature)* T) { - const struct sdis_interface* interf = NULL; - const struct sdis_medium* mdm = NULL; + struct sdis_interface* interf = NULL; + struct sdis_medium* mdm = NULL; double lambda; double delta; double delta_boundary; @@ -423,6 +445,7 @@ XD(solid_boundary_with_flux_path) float dir1[DIM]; float range[2]; int dim = DIM; + res_T res = RES_OK; ASSERT(frag && phi != SDIS_FLUX_NONE); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); (void)ctx; @@ -487,14 +510,23 @@ XD(solid_boundary_with_flux_path) /* Handle the flux */ delta_in_meter = delta*fp_to_meter; - T->value += phi * delta_in_meter / lambda; + tmp = delta_in_meter / lambda; + T->value += phi * tmp; + if(ctx->green_path) { + res = green_path_add_flux_term(ctx->green_path, interf, tmp); + if(res != RES_OK) goto error; + } /* Handle the volumic power */ power = solid_get_volumic_power(mdm, &rwalk->vtx); if(power != SDIS_VOLUMIC_POWER_NONE) { delta_in_meter = delta_boundary * fp_to_meter; - tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda); - T->value += tmp; + tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda); + T->value += power * tmp; + if(ctx->green_path) { + res = green_path_add_power_term(ctx->green_path, mdm, tmp); + if(res != RES_OK) goto error; + } } /* Reinject into the solid */ @@ -510,6 +542,11 @@ XD(solid_boundary_with_flux_path) rwalk->hit = SXD_HIT_NULL; rwalk->hit_side = SDIS_SIDE_NULL__; } + +exit: + return res; +error: + goto exit; } /******************************************************************************* @@ -525,11 +562,12 @@ XD(boundary_path) struct XD(temperature)* T) { struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; - const struct sdis_interface* interf = NULL; - const struct sdis_medium* mdm_front = NULL; - const struct sdis_medium* mdm_back = NULL; - const struct sdis_medium* mdm = NULL; + struct sdis_interface* interf = NULL; + struct sdis_medium* mdm_front = NULL; + struct sdis_medium* mdm_back = NULL; + struct sdis_medium* mdm = NULL; double tmp; + res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); ASSERT(rwalk->mdm == NULL); ASSERT(!SXD_HIT_NONE(&rwalk->hit)); @@ -546,7 +584,15 @@ XD(boundary_path) if(tmp >= 0) { T->value += tmp; T->done = 1; - return RES_OK; + + if(ctx->green_path) { + double pos[3] = {0,0,0}; + dX(set)(pos, rwalk->vtx.P); + res = green_path_set_interface_limit_vertex + (ctx->green_path, interf, pos, rwalk->vtx.time); + if(res != RES_OK) goto error; + } + goto exit; } /* Check if the boundary flux is known. Note that currently, only solid media @@ -555,9 +601,11 @@ XD(boundary_path) if(sdis_medium_get_type(mdm) == SDIS_SOLID ) { const double phi = interface_side_get_flux(interf, &frag); if(phi != SDIS_FLUX_NONE) { - XD(solid_boundary_with_flux_path) + res = XD(solid_boundary_with_flux_path) (scn, fp_to_meter, ctx, &frag, phi, rwalk, rng, T); - return RES_OK; + if(res != RES_OK) goto error; + + goto exit; } } @@ -565,13 +613,18 @@ XD(boundary_path) mdm_back = interface_get_medium(interf, SDIS_BACK); if(mdm_front->type == mdm_back->type) { - XD(solid_solid_boundary_path) + res = XD(solid_solid_boundary_path) (scn, fp_to_meter, ctx, &frag, rwalk, rng, T); } else { - XD(solid_fluid_boundary_path) + res = XD(solid_fluid_boundary_path) (scn, fp_to_meter, ctx, &frag, rwalk, rng, T); } - return RES_OK; + if(res != RES_OK) goto error; + +exit: + return res; +error: + goto exit; } #include "sdis_Xd_end.h" diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -14,6 +14,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "sdis_device_c.h" +#include "sdis_green.h" #include "sdis_heat_path.h" #include "sdis_medium_c.h" #include "sdis_misc.h" @@ -33,7 +34,8 @@ XD(conductive_path) struct XD(temperature)* T) { double position_start[DIM]; - const struct sdis_medium* mdm; + struct sdis_medium* mdm; + res_T res = RES_OK; ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T); ASSERT(rwalk->mdm->type == SDIS_SOLID); (void)ctx; @@ -43,7 +45,8 @@ XD(conductive_path) if(mdm != rwalk->mdm) { log_err(scn->dev, "%s: invalid solid random walk. " "Unexpected medium at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P)); - return RES_BAD_OP_IRRECOVERABLE; + res = RES_BAD_OP_IRRECOVERABLE; + goto error; } /* Save the submitted position */ dX(set)(position_start, rwalk->vtx.P); @@ -55,6 +58,7 @@ XD(conductive_path) double rho; /* Volumic mass */ double cp; /* Calorific capacity */ double tmp; + double power_factor; double power; float delta, delta_solid; /* Random walk numerical parameter */ float range[2]; @@ -66,7 +70,15 @@ XD(conductive_path) if(tmp >= 0) { T->value += tmp; T->done = 1; - return RES_OK; + + if(ctx->green_path) { + double pos[3] = {0,0,0}; + dX(set)(pos, rwalk->vtx.P); + res = green_path_set_medium_limit_vertex + (ctx->green_path, rwalk->mdm, pos, rwalk->vtx.time); + if(res != RES_OK) goto error; + } + goto exit; } /* Fetch solid properties */ @@ -100,8 +112,8 @@ XD(conductive_path) /* Add the volumic power density to the measured temperature */ if(power != SDIS_VOLUMIC_POWER_NONE) { const double delta_in_meter = delta * fp_to_meter; - tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); - T->value += tmp; + power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); + T->value += power * power_factor; } } else { /* Hit something: move along dir0 of the minimum hit distance */ @@ -128,7 +140,7 @@ XD(conductive_path) h_in_meter = h * fp_to_meter; /* The regular power term at wall */ - tmp = power * h_in_meter * h_in_meter / (2.0 * lambda); + tmp = h_in_meter * h_in_meter / (2.0 * lambda); /* Add the power corrective term */ if(h < delta_solid) { @@ -136,25 +148,32 @@ XD(conductive_path) #if DIM==2 /* tmp1 = sin(2a) / (PI - 2*a) */ const double tmp1 = sin_a * sqrt(1 - sin_a*sin_a)/acos(sin_a); - tmp += -(power*delta_s_in_meter*delta_s_in_meter)/(4.0*lambda) * tmp1; + tmp += -(delta_s_in_meter * delta_s_in_meter)/(4.0*lambda) * tmp1; #else const double tmp1 = (sin_a*sin_a*sin_a - sin_a)/ (1-sin_a); - tmp += (power*delta_s_in_meter*delta_s_in_meter)/(6*lambda) * tmp1; + tmp += (delta_s_in_meter * delta_s_in_meter)/(6.0*lambda) * tmp1; #endif } else if(h == delta_solid) { - tmp += -(delta_s_in_meter*delta_s_in_meter*power)/(2.0*DIM*lambda); + tmp += -(delta_s_in_meter * delta_s_in_meter)/(2.0*DIM*lambda); } - T->value += tmp; + power_factor = tmp; + T->value += power * power_factor; } } + /* Register the power term against the green function */ + if(ctx->green_path && power != SDIS_VOLUMIC_POWER_NONE) { + res = green_path_add_power_term(ctx->green_path, mdm, power_factor); + if(res != RES_OK) goto error; + } + /* Sample the time */ if(rwalk->vtx.time != INF) { double tau, mu, t0; mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter); tau = ssp_ran_exp(rng, mu); - t0 = solid_get_t0(rwalk->mdm); + t0 = ctx->green_path ? -INF : solid_get_t0(rwalk->mdm); rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); if(rwalk->vtx.time == t0) { /* Check the initial condition */ @@ -162,14 +181,15 @@ XD(conductive_path) if(tmp >= 0) { T->value += tmp; T->done = 1; - 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; } } @@ -220,7 +240,8 @@ XD(conductive_path) } #undef VEC_STR #undef VEC_SPLIT - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } /* Keep going while the solid random walk does not hit an interface */ @@ -228,7 +249,11 @@ XD(conductive_path) T->func = XD(boundary_path); rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */ - return RES_OK; + +exit: + return res; +error: + goto exit; } #include "sdis_Xd_end.h" diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -14,6 +14,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "sdis_device_c.h" +#include "sdis_green.h" #include "sdis_heat_path.h" #include "sdis_medium_c.h" #include "sdis_scene_c.h" @@ -47,6 +48,7 @@ XD(convective_path) #else float st[2]; #endif + res_T res = RES_OK; (void)rng, (void)fp_to_meter, (void)ctx; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); ASSERT(rwalk->mdm->type == SDIS_FLUID); @@ -55,7 +57,20 @@ XD(convective_path) if(tmp >= 0) { /* T is known. */ T->value += tmp; T->done = 1; - return RES_OK; + + if(ctx->green_path) { + double pos[3] = {0,0,0}; + dX(set)(pos, rwalk->vtx.P); + res = green_path_set_medium_limit_vertex + (ctx->green_path, rwalk->mdm, pos, rwalk->vtx.time); + if(res != RES_OK) { + log_err(scn->dev, + "%s: could not register the limit vertex of a sampled path " + "against the green function.\n", FUNC_NAME); + goto error; + } + } + goto exit; } if(SXD_HIT_NONE(&rwalk->hit)) { /* The path begins in the fluid */ @@ -75,7 +90,8 @@ XD(convective_path) log_err(scn->dev, "%s: the position %g %g %g lies in the surrounding fluid whose temperature must \n" "be known.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P)); - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } } @@ -103,7 +119,8 @@ XD(convective_path) log_err(scn->dev, "%s: invalid enclosure. The surrounding fluid has an unset temperature.\n", FUNC_NAME); - return RES_BAD_ARG; + res = RES_BAD_ARG; + goto error; } /* The hc upper bound can be 0 if h is uniformly 0. In that case the result @@ -111,12 +128,22 @@ XD(convective_path) if(enc->hc_upper_bound == 0) { /* Cannot be in the fluid without starting there. */ ASSERT(SXD_HIT_NONE(&rwalk->hit)); + + if(ctx->green_path) { + log_err(scn->dev, + "%s: the upper bound of the convection cannot of an enclosure cannot be " + "null when registering the green function; initial condition is not " + "supported.\n", FUNC_NAME); + res = RES_BAD_ARG; + goto error; + } + rwalk->vtx.time = fluid_get_t0(rwalk->mdm); tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx); if(tmp >= 0) { T->value += tmp; T->done = 1; - return RES_OK; + goto exit; } /* At t=0, the initial condition should have been reached. */ @@ -124,10 +151,12 @@ XD(convective_path) "%s: undefined initial condition. " "Time is 0 but the temperature remains unknown.\n", FUNC_NAME); - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } - /* A trick to force first r test result. */ + /* 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. */ @@ -142,7 +171,8 @@ XD(convective_path) 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)); - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } if(r < hc / enc->hc_upper_bound) { @@ -159,7 +189,7 @@ XD(convective_path) double mu, tau, t0; mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V; tau = ssp_ran_exp(rng, mu); - t0 = fluid_get_t0(rwalk->mdm); + t0 = ctx->green_path ? -INF : fluid_get_t0(rwalk->mdm); rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); if(rwalk->vtx.time == t0) { /* Check the initial condition. */ @@ -167,14 +197,15 @@ XD(convective_path) if(tmp >= 0) { T->value += tmp; T->done = 1; - return RES_OK; + goto exit; } /* The initial condition should have been reached. */ log_err(scn->dev, "%s: undefined initial condition. " "Time is %g but the temperature remains unknown.\n", FUNC_NAME, t0); - return RES_BAD_OP; + res = RES_BAD_OP; + goto error; } } @@ -222,7 +253,11 @@ XD(convective_path) rwalk->hit.distance = 0; T->func = XD(boundary_path); rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */ - return RES_OK; + +exit: + return res; +error: + goto exit; } #include "sdis_Xd_end.h" diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -14,6 +14,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "sdis_device_c.h" +#include "sdis_green.h" #include "sdis_heat_path.h" #include "sdis_interface_c.h" #include "sdis_medium_c.h" @@ -52,7 +53,7 @@ XD(trace_radiative_path) for(;;) { const struct sdis_interface* interf = NULL; struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; - const struct sdis_medium* chk_mdm = NULL; + struct sdis_medium* chk_mdm = NULL; double alpha; double epsilon; double r; @@ -74,6 +75,13 @@ XD(trace_radiative_path) if(ctx->Tarad >= 0) { T->value += ctx->Tarad; T->done = 1; + + if(ctx->green_path) { + const double inf_pos[3] = {INF,INF,INF}; + res = green_path_set_medium_limit_vertex + (ctx->green_path, rwalk->mdm, inf_pos, rwalk->vtx.time); + if(res != RES_OK) goto error; + } break; } else { log_err(scn->dev, diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -187,7 +187,7 @@ sdis_interface_ref_put(struct sdis_interface* interf) /******************************************************************************* * Local function ******************************************************************************/ -const struct sdis_medium* +struct sdis_medium* interface_get_medium (const struct sdis_interface* interf, const enum sdis_side side) { diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -36,7 +36,7 @@ struct sdis_interface { struct sdis_device* dev; }; -extern LOCAL_SYM const struct sdis_medium* +extern LOCAL_SYM struct sdis_medium* interface_get_medium (const struct sdis_interface* interf, const enum sdis_side side); diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c @@ -25,7 +25,7 @@ res_T ray_realisation_3d (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* medium, + struct sdis_medium* medium, const double position[], const double direction[], const double time, @@ -34,7 +34,7 @@ ray_realisation_3d const double Tref, double* weight) { - struct rwalk_context ctx; + struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct rwalk_3d rwalk = RWALK_NULL_3d; struct temperature_3d T = TEMPERATURE_NULL_3d; float dir[3]; diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h @@ -21,6 +21,8 @@ #include <rsys/rsys.h> +/* Forward declarations */ +struct green_path_handle; struct sdis_scene; struct ssp_rng; @@ -37,24 +39,26 @@ extern LOCAL_SYM res_T probe_realisation_2d (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* medium, + struct sdis_medium* medium, const double position[2], const double time, const double fp_to_meter,/* Scale factor from floating point unit to meter */ const double ambient_radiative_temperature, const double reference_temperature, + struct green_path_handle* green_path, double* weight); extern LOCAL_SYM res_T probe_realisation_3d (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* medium, + struct sdis_medium* medium, const double position[3], const double time, const double fp_to_meter,/* Scale factor from floating point unit to meter */ const double ambient_radiative_temperature, const double reference_temperature, + struct green_path_handle* green_path, double* weight); /******************************************************************************* @@ -121,7 +125,7 @@ extern LOCAL_SYM res_T ray_realisation_3d (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* medium, + struct sdis_medium* medium, const double position[3], const double direction[3], const double time, diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -90,15 +90,16 @@ res_T XD(probe_realisation) (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* medium, + struct sdis_medium* medium, const double position[], const double time, const double fp_to_meter,/* Scale factor from floating point unit to meter */ const double ambient_radiative_temperature, const double reference_temperature, + struct green_path_handle* green_path, double* weight) { - struct rwalk_context ctx; + struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct XD(rwalk) rwalk = XD(RWALK_NULL); struct XD(temperature) T = XD(TEMPERATURE_NULL); double t0; @@ -124,7 +125,8 @@ XD(probe_realisation) dX(set)(rwalk.vtx.P, position); rwalk.vtx.time = time; - if(t0 >= rwalk.vtx.time) { + /* No initial condition with green */ + if(!green_path && t0 >= rwalk.vtx.time) { double tmp; /* Check the initial condition. */ rwalk.vtx.time = t0; @@ -144,6 +146,7 @@ XD(probe_realisation) rwalk.hit = SXD_HIT_NULL; rwalk.mdm = medium; + ctx.green_path = green_path; ctx.Tarad = ambient_radiative_temperature; ctx.Tref3 = reference_temperature @@ -170,7 +173,7 @@ XD(boundary_realisation) const double Tref, double* weight) { - struct rwalk_context ctx; + struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct XD(rwalk) rwalk = XD(RWALK_NULL); struct XD(temperature) T = XD(TEMPERATURE_NULL); struct sXd(attrib) attr; @@ -236,7 +239,7 @@ XD(boundary_flux_realisation) const int flux_mask, double weight[3]) { - struct rwalk_context ctx; + struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct XD(rwalk) rwalk; struct XD(temperature) T; struct sXd(attrib) attr; @@ -308,9 +311,8 @@ XD(boundary_flux_realisation) /* Compute fluid temperature */ if(compute_convective) { - const struct sdis_interface* interf = - scene_get_interface(scn, (unsigned)iprim); - const struct sdis_medium* mdm = interface_get_medium(interf, fluid_side); + struct sdis_interface* interf = scene_get_interface(scn, (unsigned)iprim); + struct sdis_medium* mdm = interface_get_medium(interf, fluid_side); RESET_WALK(fluid_side, mdm); T.func = XD(convective_path); diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -276,7 +276,7 @@ sdis_scene_boundary_project_position /******************************************************************************* * Local miscellaneous function ******************************************************************************/ -const struct sdis_interface* +struct sdis_interface* scene_get_interface(const struct sdis_scene* scn, const unsigned iprim) { ASSERT(scn && iprim < darray_prim_prop_size_get(&scn->prim_props)); @@ -288,7 +288,7 @@ scene_get_medium (const struct sdis_scene* scn, const double pos[], struct get_medium_info* info, - const struct sdis_medium** out_medium) + struct sdis_medium** out_medium) { return scene_is_2d(scn) ? scene_get_medium_2d(scn, pos, info, out_medium) diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -801,9 +801,9 @@ XD(scene_get_medium) (const struct sdis_scene* scn, const double pos[2], struct get_medium_info* info, /* May be NULL */ - const struct sdis_medium** out_medium) + struct sdis_medium** out_medium) { - const struct sdis_medium* medium = NULL; + struct sdis_medium* medium = NULL; size_t iprim, nprims; size_t nfailures = 0; const size_t max_failures = 10; diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -203,7 +203,7 @@ scene_get_primitives_count(const struct sdis_scene* scn) return darray_prim_prop_size_get(&scn->prim_props); } -extern LOCAL_SYM const struct sdis_interface* +extern LOCAL_SYM struct sdis_interface* scene_get_interface (const struct sdis_scene* scene, const unsigned iprim); @@ -213,7 +213,7 @@ scene_get_medium (const struct sdis_scene* scene, const double position[], struct get_medium_info* info, /* May be NULL */ - const struct sdis_medium** medium); + struct sdis_medium** medium); static INLINE void scene_get_enclosure_ids diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -46,7 +46,7 @@ static res_T solve_pixel (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* mdm, + struct sdis_medium* mdm, const struct sdis_camera* cam, const double time, /* Observation time */ const double fp_to_meter, /* Scale from floating point units to meters */ @@ -108,7 +108,7 @@ static res_T solve_tile (struct sdis_scene* scn, struct ssp_rng* rng, - const struct sdis_medium* mdm, + struct sdis_medium* mdm, const struct sdis_camera* cam, const double time, const double fp_to_meter, @@ -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 */ + struct sdis_green_function** out_green, /* May be NULL<=>Do not store green */ 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, out_estimator); + fp_to_meter, Tarad, Tref, out_green, out_estimator); } else { return solve_probe_3d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, out_estimator); + fp_to_meter, Tarad, Tref, out_green, out_estimator); } } @@ -288,7 +289,7 @@ sdis_solve_camera #define TILE_SIZE 32 /* definition in X & Y of a tile */ STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2); - const struct sdis_medium* medium = NULL; + struct sdis_medium* medium = NULL; struct darray_accum* tiles = NULL; struct ssp_rng_proxy* rng_proxy = NULL; struct ssp_rng** rngs = NULL; diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -15,6 +15,7 @@ #include "sdis_device_c.h" #include "sdis_estimator_c.h" +#include "sdis_green.h" #include "sdis_medium_c.h" #include "sdis_misc.h" #include "sdis_realisation.h" @@ -170,10 +171,13 @@ 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 */ + struct sdis_green_function** out_green, /* May be NULL <=> No green func */ struct sdis_estimator** out_estimator) { - const struct sdis_medium* medium = NULL; + struct sdis_medium* medium = 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; double weight = 0; @@ -204,12 +208,8 @@ XD(solve_probe) if(res != RES_OK) goto error; /* Create the per thread RNG */ - rngs = MEM_CALLOC - (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*)); - if(!rngs) { - res = RES_MEM_ERR; - goto error; - } + rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); + if(!rngs) { res = RES_MEM_ERR; goto error; } FOR_EACH(i, 0, scn->dev->nthreads) { res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i); if(res != RES_OK) goto error; @@ -219,6 +219,17 @@ XD(solve_probe) res = scene_get_medium(scn, position, NULL, &medium); 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; + } + + } + /* 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) @@ -228,18 +239,27 @@ XD(solve_probe) double time; const int ithread = omp_get_thread_num(); struct ssp_rng* rng = rngs[ithread]; + struct green_path_handle* pgreen_path = NULL; + struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL; if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */ - time = sample_time(rng, time_range); + if(!out_green) { + time = sample_time(rng, time_range); + } 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(probe_realisation) - (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w); + res_local = XD(probe_realisation)(scn, rng, medium, position, time, + fp_to_meter, Tarad, Tref, pgreen_path, &w); if(res_local != RES_OK) { - if(res_local != RES_BAD_OP) { - ATOMIC_SET(&res, res_local); - continue; - } + if(res_local != RES_BAD_OP) { ATOMIC_SET(&res, res_local); continue; } } else { weight += w; sqr_weight += w*w; @@ -256,6 +276,15 @@ XD(solve_probe) /* Setup the estimated temperature */ estimator_setup_temperature(estimator, weight, sqr_weight); + 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; + } + } + exit: if(rngs) { FOR_EACH(i, 0, scn->dev->nthreads) { @@ -263,10 +292,21 @@ 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(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_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -396,7 +396,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, NULL, &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 @@ -401,7 +401,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, NULL, &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 @@ -266,7 +266,7 @@ main(int argc, char** argv) ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * 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, NULL, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -289,7 +289,7 @@ main(int argc, char** argv) ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * 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, NULL, &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 @@ -281,7 +281,7 @@ main(int argc, char** argv) ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * 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, NULL, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -303,7 +303,7 @@ main(int argc, char** argv) time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * 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, NULL, &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 @@ -234,7 +234,7 @@ main(int argc, char** argv) ref = T0 + (1 - pos[0]) * PHI/LAMBDA; /* 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, NULL, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -247,7 +247,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE*3)); /* 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, NULL, &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.c b/src/test_sdis_solve_probe.c @@ -261,13 +261,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, NULL, &estimator)); + BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, NULL, &estimator)); + BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, NULL)); + OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, NULL, &estimator)); BA(sdis_estimator_get_type(estimator, NULL)); BA(sdis_estimator_get_type(NULL, &type)); @@ -316,7 +316,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, NULL, &estimator)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -240,7 +240,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, NULL, &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 @@ -238,7 +238,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, NULL, &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 @@ -296,7 +296,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, NULL, &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 @@ -289,7 +289,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, NULL, &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 @@ -199,7 +199,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, NULL, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); @@ -219,7 +219,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, NULL, &estimator)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -268,7 +268,7 @@ main(int argc, char** argv) ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; /* 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, NULL, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -282,7 +282,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE)); /* 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, NULL, &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_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, NULL, &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_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, NULL, &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, NULL, &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, NULL, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails));