commit 97b6e62c27c7a9c7d7590612a5c4bb40b2be58a4
parent bb0c8f4c24b48cd1630feffe15b8e5940c69930b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 2 Oct 2020 10:11:35 +0200
Merge branch 'release_0.10.1'
Diffstat:
23 files changed, 163 insertions(+), 78 deletions(-)
diff --git a/README.md b/README.md
@@ -25,6 +25,18 @@ variable the install directories of its dependencies.
## Release notes
+### Version 0.10.1
+
+- In green function estimation, the time sent to the user callbacks is no more
+ the elapsed time from the beginning of the realisation: as in a regular
+ computation, it is now the observation time.
+- Fix the flux computation for boundaries with an imposed flux: it was
+ previously ignored. The new `sdis_estimator_get_imposed_flux` function
+ returns this estimated flux component.
+- Return an error if the flux is computed at a boundary whose temperature is
+ known: this configuration is not currently supported.
+- Fix build with the CL compiler.
+
### Version 0.10
- Add support of green function [de]serialization. The
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -57,7 +57,7 @@ rcmake_append_runtime_dirs(_runtime_dirs
###############################################################################
set(VERSION_MAJOR 0)
set(VERSION_MINOR 10)
-set(VERSION_PATCH 0)
+set(VERSION_PATCH 1)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SDIS_FILES_SRC
diff --git a/src/sdis.h b/src/sdis.h
@@ -940,6 +940,11 @@ sdis_estimator_get_radiative_flux
struct sdis_mc* flux);
SDIS_API res_T
+sdis_estimator_get_imposed_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);
@@ -1144,11 +1149,8 @@ sdis_compute_power
/*******************************************************************************
* Green solvers.
*
- * The caller should ensure that green solvers are invoked on scenes whose data
- * do not depend on time. Indeed, on green estimation, the time parameter along
- * the random walks registers the relative time spent in the system rather than
- * an absolute time. As a consequence, the media/interfaces parameters cannot
- * vary in time with respect to an absolute time value.
+ * Currently only steady computations are supported. As a consequence, the
+ * observation time is always fixed to infinity.
*
* In addition, the green solvers assumes that the interface fluxes are
* constants in time and space. In the same way the volumic power of the solid
@@ -1159,8 +1161,8 @@ sdis_compute_power
* definitely registered against the green function as interfaces/media with no
* flux/volumic power.
*
- * If the aforementioned assumptions are not ensured by the caller, the
- * behavior of the estimated green function is undefined.
+ * If these assumptions are not ensured by the caller, the behavior of the
+ * estimated green function is undefined.
******************************************************************************/
SDIS_API res_T
sdis_solve_probe_green_function
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -95,10 +95,11 @@ struct XD(rwalk) {
struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
struct sdis_medium* mdm; /* Medium in which the random walk lies */
struct sXd(hit) hit; /* Hit of the random walk */
+ double elapsed_time;
enum sdis_side hit_side;
};
static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, SDIS_SIDE_NULL__
+ SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, 0, SDIS_SIDE_NULL__
};
struct XD(temperature) {
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -136,6 +136,16 @@ sdis_estimator_get_radiative_flux
}
res_T
+sdis_estimator_get_imposed_flux
+ (const struct sdis_estimator* estimator, struct sdis_mc* flux)
+{
+ if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
+ return RES_BAD_ARG;
+ SETUP_MC(flux, &estimator->fluxes[FLUX_IMPOSED]);
+ return RES_OK;
+}
+
+res_T
sdis_estimator_get_total_flux
(const struct sdis_estimator* estimator, struct sdis_mc* flux)
{
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -31,6 +31,7 @@ enum sdis_estimator_type;
enum flux_name {
FLUX_CONVECTIVE,
FLUX_RADIATIVE,
+ FLUX_IMPOSED,
FLUX_TOTAL,
FLUX_NAMES_COUNT__
};
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -1462,7 +1462,8 @@ res_T
green_path_set_limit_interface_fragment
(struct green_path_handle* handle,
struct sdis_interface* interf,
- const struct sdis_interface_fragment* frag)
+ const struct sdis_interface_fragment* frag,
+ const double elapsed_time)
{
res_T res = RES_OK;
ASSERT(handle && interf && frag);
@@ -1470,6 +1471,7 @@ green_path_set_limit_interface_fragment
res = ensure_interface_registration(handle->green, interf);
if(res != RES_OK) return res;
handle->path->limit.fragment = *frag;
+ handle->path->limit.fragment.time = -elapsed_time;
handle->path->limit_id = interface_get_id(interf);
handle->path->limit_type = SDIS_FRAGMENT;
return RES_OK;
@@ -1479,7 +1481,8 @@ res_T
green_path_set_limit_vertex
(struct green_path_handle* handle,
struct sdis_medium* mdm,
- const struct sdis_rwalk_vertex* vert)
+ const struct sdis_rwalk_vertex* vert,
+ const double elapsed_time)
{
res_T res = RES_OK;
ASSERT(handle && mdm && vert);
@@ -1487,6 +1490,7 @@ green_path_set_limit_vertex
res = ensure_medium_registration(handle->green, mdm);
if(res != RES_OK) return res;
handle->path->limit.vertex = *vert;
+ handle->path->limit.vertex.time = -elapsed_time;
handle->path->limit_id = medium_get_id(mdm);
handle->path->limit_type = SDIS_VERTEX;
return RES_OK;
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -71,13 +71,15 @@ extern LOCAL_SYM res_T
green_path_set_limit_interface_fragment
(struct green_path_handle* path,
struct sdis_interface* interf,
- const struct sdis_interface_fragment* fragment);
+ const struct sdis_interface_fragment* fragment,
+ const double elapsed_time);
extern LOCAL_SYM res_T
green_path_set_limit_vertex
(struct green_path_handle* path,
struct sdis_medium* mdm,
- const struct sdis_rwalk_vertex* vertex);
+ const struct sdis_rwalk_vertex* vertex,
+ const double elapsed_time);
extern LOCAL_SYM res_T
green_path_add_power_term
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -977,7 +977,7 @@ XD(boundary_path)
if(ctx->green_path) {
res = green_path_set_limit_interface_fragment
- (ctx->green_path, interf, &frag);
+ (ctx->green_path, interf, &frag, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
if(ctx->heat_path) {
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -251,7 +251,7 @@ XD(conductive_path)
if(ctx->green_path) {
res = green_path_set_limit_vertex
- (ctx->green_path, rwalk->mdm, &rwalk->vtx);
+ (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -103,7 +103,8 @@ XD(convective_path)
T->done = 1;
if(ctx->green_path) {
- res = green_path_set_limit_vertex(ctx->green_path, rwalk->mdm, &rwalk->vtx);
+ res = green_path_set_limit_vertex
+ (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
@@ -201,18 +202,23 @@ XD(convective_path)
for(;;) {
struct sdis_interface_fragment frag;
struct sXd(primitive) prim;
+ double mu, tau, t0;
/* Fetch other physical properties. */
cp = fluid_get_calorific_capacity(rwalk->mdm, &rwalk->vtx);
rho = fluid_get_volumic_mass(rwalk->mdm, &rwalk->vtx);
+ t0 = fluid_get_t0(rwalk->mdm); /* Limit time */
/* Sample the time using the upper bound. */
+ mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V;
+ tau = ssp_ran_exp(rng, mu);
+
+ /* Increment the elapsed time */
+ ASSERT(rwalk->vtx.time > t0);
+ rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0);
+
if(rwalk->vtx.time != INF) {
- double mu, tau, t0;
- mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V;
- 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);
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); /* Time rewind */
/* Register the new vertex against the heat path */
res = XD(register_heat_vertex_in_fluid)(scn, ctx, rwalk, T->value);
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -84,7 +84,8 @@ XD(trace_radiative_path)
struct sdis_rwalk_vertex vtx;
d3_splat(vtx.P, INF);
vtx.time = rwalk->vtx.time;
- res = green_path_set_limit_vertex(ctx->green_path, rwalk->mdm, &vtx);
+ res = green_path_set_limit_vertex
+ (ctx->green_path, rwalk->mdm, &vtx, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
if(ctx->heat_path) {
diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h
@@ -41,20 +41,22 @@ XD(time_rewind)
ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
ASSERT(T->done == 0);
- if(IS_INF(rwalk->vtx.time)) goto exit;
-
- /* Fetch phyisical properties */
+ /* Fetch physical properties */
lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
-
- /* Fetch the limit time */
- t0 = ctx->green_path ? -INF : solid_get_t0(mdm);
+ t0 = solid_get_t0(mdm); /* Limit time */
/* Sample the time to reroll */
mu = (2*DIM*lambda)/(rho*cp*delta_in_meter*delta_in_meter);
tau = ssp_ran_exp(rng, mu);
+ /* Increment the elapsed time */
+ ASSERT(rwalk->vtx.time > t0);
+ rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0);
+
+ if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */
+
/* Time rewind */
rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0);
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -31,7 +31,9 @@
#include <float.h>
#include <limits.h>
+#ifdef COMPILER_GCC
#include <sys/mman.h>
+#endif
/*******************************************************************************
* Helper function
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -264,8 +264,8 @@ XD(solve_boundary)
}
} else {
/* Do not take care of the submitted time when registering the green
- * function. Simply takes 0 as relative time */
- time = 0;
+ * function. Only steady systems are supported yet */
+ time = INF;
res_local = green_function_create_path(greens[ithread], &green_path);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
@@ -446,6 +446,7 @@ XD(solve_boundary_flux)
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
+ struct accum* acc_fi = NULL; /* Per thread imposed flux accumulator */
size_t nrealisations = 0;
int64_t irealisation;
size_t i;
@@ -552,6 +553,7 @@ XD(solve_boundary_flux)
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
+ ALLOC_ACCUMS(acc_fi);
#undef ALLOC_ACCUMS
/* Create the estimator */
@@ -570,6 +572,7 @@ XD(solve_boundary_flux)
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
+ struct accum* acc_fimp = &acc_fi[ithread];
struct sXd(primitive) prim;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
const struct sdis_interface* interf;
@@ -578,7 +581,7 @@ XD(solve_boundary_flux)
double T_brf[3] = { 0, 0, 0 };
const double Tref = args->reference_temperature;
const double Tarad = args->ambient_radiative_temperature;
- double epsilon, hc, hr;
+ double epsilon, hc, hr, imposed_flux, imposed_temp;
size_t iprim;
double uv[DIM - 1];
float st[DIM - 1];
@@ -624,8 +627,11 @@ XD(solve_boundary_flux)
bmd = interface_get_medium(interf, SDIS_BACK);
if(!fmd || !bmd
|| ( !(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
- && !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID)))
- {
+ && !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID))) {
+ log_err(scn->dev, "%s: Attempt to compute a flux at a %s-%s interface.\n",
+ FUNC_NAME,
+ (fmd->type == SDIS_FLUID ? "fluid" : "solid"),
+ (bmd->type == SDIS_FLUID ? "fluid" : "solid"));
ATOMIC_SET(&res, RES_BAD_ARG);
continue;
}
@@ -640,14 +646,22 @@ XD(solve_boundary_flux)
/* Fetch interface parameters */
epsilon = interface_side_get_emissivity(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
-
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ frag.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag);
+ imposed_temp = interface_side_get_temperature(interf, &frag);
+ if(imposed_temp >= 0) {
+ /* Flux computation on T boundaries is not supported yet */
+ log_err(scn->dev, "%s: Attempt to compute a flux at a Dirichlet boundary "
+ "(not available yet).\n", FUNC_NAME);
+ ATOMIC_SET(&res, RES_BAD_ARG);
+ continue;
+ }
/* Fluid, Radiative and Solid temperatures */
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
-
res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
@@ -664,7 +678,8 @@ XD(solve_boundary_flux)
const double Tfluid = T_brf[2];
const double w_conv = hc * (Tboundary - Tfluid);
const double w_rad = hr * (Tboundary - Tradiative);
- const double w_total = w_conv + w_rad;
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
+ const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
@@ -685,6 +700,10 @@ XD(solve_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
+ /* Imposed flux */
+ acc_fimp->sum += w_imp;
+ acc_fimp->sum2 += w_imp*w_imp;
+ ++acc_fimp->count;
}
/* Update progress */
@@ -707,10 +726,12 @@ XD(solve_boundary_flux)
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
+ sum_accums(acc_fi, scn->dev->nthreads, &acc_fi[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
+ ASSERT(acc_tp[0].count == acc_fi[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
@@ -718,6 +739,7 @@ XD(solve_boundary_flux)
estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
+ estimator_setup_flux(estimator, FLUX_IMPOSED, acc_fi[0].sum, acc_fi[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
res = estimator_save_rng_state(estimator, rng_proxy);
@@ -733,6 +755,7 @@ exit:
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
+ if(acc_fi) MEM_RM(scn->dev->allocator, acc_fi);
if(scene) SXD(scene_ref_put(scene));
if(shape) SXD(shape_ref_put(shape));
if(view) SXD(scene_view_ref_put(view));
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -333,8 +333,8 @@ XD(solve_medium)
}
} else {
/* Do not take care of the submitted time when registering the green
- * function. Simply takes 0 as relative time */
- time = 0;
+ * function. Only steady systems are supported yet */
+ time = INF;
res_local = green_function_create_path(greens[ithread], &green_path);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -161,8 +161,8 @@ XD(solve_probe)
}
} else {
/* Do not take care of the submitted time when registering the green
- * function. Simply takes 0 as relative time */
- time = 0;
+ * function. Only steady systems are supported */
+ time = INF;
res_local = green_function_create_path(greens[ithread], &green_path);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -196,8 +196,8 @@ XD(solve_probe_boundary)
}
} else {
/* Do not take care of the submitted time when registering the green
- * function. Simply takes 0 as relative time */
- time = 0;
+ * function. Only steady systems are supported */
+ time = INF;
res_local = green_function_create_path(greens[ithread], &green_path);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
@@ -346,6 +346,7 @@ XD(solve_probe_boundary_flux)
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
+ struct accum* acc_fi = NULL; /* Per thread imposed flux accumulator */
size_t nrealisations = 0;
int64_t irealisation = 0;
size_t i;
@@ -393,9 +394,9 @@ XD(solve_probe_boundary_flux)
}
} else {
const double w = CLAMP(1 - args->uv[0] - args->uv[1], 0, 1);
- if(args->uv[0] < 0
- || args->uv[1] < 0
- || args->uv[0] > 1
+ if(args->uv[0] < 0
+ || args->uv[1] < 0
+ || args->uv[0] > 1
|| args->uv[1] > 1
|| !eq_eps(w + args->uv[0] + args->uv[1], 1, 1.e-6)) {
log_err(scn->dev,
@@ -413,6 +414,11 @@ XD(solve_probe_boundary_flux)
if(!fmd || !bmd
|| ( !(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID)
&& !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID))) {
+ log_err(scn->dev,
+ "%s: Attempt to compute a flux at a %s-%s interface.\n",
+ FUNC_NAME,
+ (fmd->type == SDIS_FLUID ? "fluid" : "solid"),
+ (bmd->type == SDIS_FLUID ? "fluid" : "solid"));
res = RES_BAD_ARG;
goto error;
}
@@ -449,6 +455,7 @@ XD(solve_probe_boundary_flux)
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
+ ALLOC_ACCUMS(acc_fi);
#undef ALLOC_ACCUMS
/* Prebuild the interface fragment */
@@ -473,7 +480,8 @@ XD(solve_probe_boundary_flux)
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
- double time, epsilon, hc, hr;
+ struct accum* acc_fimp = &acc_fi[ithread];
+ double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
double T_brf[3] = { 0, 0, 0 };
const double Tref = args->reference_temperature;
@@ -491,16 +499,27 @@ XD(solve_probe_boundary_flux)
/* Compute hr and hc */
frag.time = time;
+ frag.side = fluid_side;
epsilon = interface_side_get_emissivity(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ frag.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag);
+ imposed_temp = interface_side_get_temperature(interf, &frag);
+ if(imposed_temp >= 0) {
+ /* Flux computation on T boundaries is not supported yet */
+ log_err(scn->dev,"%s: Attempt to compute a flux at a Dirichlet boundary "
+ "(not available yet).\n", FUNC_NAME);
+ ATOMIC_SET(&res, RES_BAD_ARG);
+ continue;
+ }
/* Fluid, Radiative and Solid temperatures */
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv, time,
- solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
+ res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv,
+ time, solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
@@ -515,7 +534,8 @@ XD(solve_probe_boundary_flux)
const double Tfluid = T_brf[2];
const double w_conv = hc * (Tboundary - Tfluid);
const double w_rad = hr * (Tboundary - Tradiative);
- const double w_total = w_conv + w_rad;
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
+ const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
@@ -536,6 +556,10 @@ XD(solve_probe_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
+ /* Imposed flux */
+ acc_fimp->sum += w_imp;
+ acc_fimp->sum2 += w_imp*w_imp;
+ ++acc_fimp->count;
}
/* Update progress */
@@ -558,10 +582,12 @@ XD(solve_probe_boundary_flux)
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
+ sum_accums(acc_fi, scn->dev->nthreads, &acc_fi[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
+ ASSERT(acc_tp[0].count == acc_fi[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
@@ -569,6 +595,7 @@ XD(solve_probe_boundary_flux)
estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
+ estimator_setup_flux(estimator, FLUX_IMPOSED, acc_fi[0].sum, acc_fi[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
res = estimator_save_rng_state(estimator, rng_proxy);
@@ -584,6 +611,7 @@ exit:
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
+ if(acc_fi) MEM_RM(scn->dev->allocator, acc_fi);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
return (res_T)res;
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -62,7 +62,7 @@
*/
#define UNKNOWN_TEMPERATURE -1
-#define N 10000 /* #realisations */
+#define N 100000 /* #realisations */
#define Tf 300.0
#define Tb 0.0
@@ -447,29 +447,17 @@ main(int argc, char** argv)
printf("Average values of the right side of the square = ");
check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
OK(sdis_estimator_ref_put(estimator));
-
- /* Average temperature on the left side of the box */
+
+ /* Flux computation on Dirichlet boundaries is not available yet.
+ * Once available, the expected total flux is the same we expect on the right
+ * side (as the other sides are adiabatic). */
prims[0] = 2;
prims[1] = 3;
-
- analyticT = Tb;
- analyticCF = H * (analyticT - Tf);
- analyticRF = Hrad * (analyticT - Trad);
- analyticTF = analyticCF + analyticRF;
-
bound_args.nprimitives = 2;
- OK(SOLVE(box_scn, &bound_args, &estimator));
- printf("Average values of the left side of the box = ");
- check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
- OK(sdis_estimator_ref_put(estimator));
-
- /* Average temperature on the left/right side of the square */
+ BA(SOLVE(box_scn, &bound_args, &estimator));
prims[0] = 1;
bound_args.nprimitives = 1;
- OK(SOLVE(square_scn, &bound_args, &estimator));
- printf("Average values of the left side of the square = ");
- check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF);
- OK(sdis_estimator_ref_put(estimator));
+ BA(SOLVE(square_scn, &bound_args, &estimator));
#undef SOLVE
OK(sdis_scene_ref_put(box_scn));
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -472,7 +472,8 @@ main(int argc, char** argv)
check_green_function(green);
check_estimator_eq(estimator, estimator2);
- CHK(stream = tmpfile());
+ stream = tmpfile();
+ CHK(stream);
BA(sdis_green_function_write(NULL, stream));
BA(sdis_green_function_write(green, NULL));
OK(sdis_green_function_write(green, stream));
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -74,7 +74,7 @@ static double
temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
- CHK(vtx != NULL);
+ CHK(vtx != NULL && IS_INF(vtx->time));
return -1;
}
@@ -83,7 +83,8 @@ solid_get_calorific_capacity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
- CHK(vtx != NULL && data == NULL);
+ CHK(vtx != NULL && IS_INF(vtx->time) && data == NULL);
+ CHK(IS_INF(vtx->time));
return 2.0;
}
@@ -92,7 +93,7 @@ solid_get_thermal_conductivity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
- CHK(vtx != NULL);
+ CHK(vtx != NULL && IS_INF(vtx->time));
return 50.0;
}
@@ -101,7 +102,7 @@ solid_get_volumic_mass
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
- CHK(vtx != NULL);
+ CHK(vtx != NULL && IS_INF(vtx->time));
return 25.0;
}
@@ -110,7 +111,7 @@ solid_get_delta
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
- CHK(vtx != NULL);
+ CHK(vtx != NULL && IS_INF(vtx->time));
return 1.0/20.0;
}
@@ -125,7 +126,7 @@ static double
null_interface_value
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
- CHK(frag != NULL);
+ CHK(frag != NULL && IS_INF(frag->time));
(void)data;
return 0;
}
@@ -134,7 +135,7 @@ static double
interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
- CHK(data != NULL && frag != NULL);
+ CHK(data != NULL && frag != NULL && IS_INF(frag->time));
return ((const struct interf*)sdis_data_cget(data))->temperature;
}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -376,7 +376,8 @@ check_green_serialization
struct sdis_green_function* green2 = NULL;
CHK(green && time_range);
- CHK(stream = tmpfile());
+ stream = tmpfile();
+ CHK(stream);
OK(sdis_green_function_write(green, stream));
diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c
@@ -22,7 +22,7 @@
#define Power 10000.0
#define H 50.0
#define LAMBDA 100.0
-#define DELTA 0.4/*(1.0/2.0)*/
+#define DELTA 0.2/*(1.0/2.0)*/
#define N 100000
#define LENGTH 10000.0
