commit 7c3bcd6a15ebe09c469ede18a18c9e504a4972c1
parent debdd7e6ec7c5af0440e738ce12e27187d76c884
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 2 Oct 2020 09:35:06 +0200
Set the observation time to INF in green estimation
Previously, during green estimation, the time submitted to the shader
callbacks was negative and registered the elapsed time from the
beginning of the realisation. This elapsed time is now internally saved
in a separated member variable of the random walk data structure. The
observation time submitted to the user callbacks is now the right one,
i.e. the infinity since only steady computations are handled during
green estimations.
Diffstat:
16 files changed, 58 insertions(+), 44 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -1149,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
@@ -1164,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_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_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);
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);
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
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_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
