commit 85f0a06d1399f6c08b44fb51f307416b86d88f5a
parent 8941ba92d7d4480f4f66622235ed1c64dff13e7b
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 2 Oct 2020 11:22:27 +0200
Merge branch 'develop' into feature_unsteady_green
Diffstat:
16 files changed, 83 insertions(+), 70 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
@@ -1149,10 +1149,17 @@ sdis_compute_power
* Green solvers.
*
* Note that only the interfaces/media with flux/volumic power defined during
- * green estimation can update their flux/volumic power values for subsequent
+ * green estimation can update their flux/volumic power value for subsequent
* sdis_green_function_solve invocations: others interfaces/media are
* definitely registered against the green function as interfaces/media with no
* flux/volumic power.
+ *
+ * Also note that the green solvers assume that the interface fluxes are
+ * constant in time and space. The same applies to the volumic power of the
+ * solid media.
+ *
+ * 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
@@ -1447,7 +1447,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);
@@ -1455,6 +1456,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;
@@ -1464,7 +1466,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);
@@ -1472,6 +1475,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 = 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
@@ -42,20 +42,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 = 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);
@@ -84,7 +86,8 @@ XD(time_rewind)
}
if(ctx->green_path) {
- res = green_path_set_limit_vertex(ctx->green_path, mdm, &rwalk->vtx);
+ res = green_path_set_limit_vertex(ctx->green_path, mdm, &rwalk->vtx,
+ rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -448,7 +448,6 @@ XD(solve_boundary_flux)
size_t i;
size_t view_nprims;
int progress = 0;
- int msg1 = 0, msg2 = 0;
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
@@ -559,7 +558,7 @@ XD(solve_boundary_flux)
nrealisations = args->nrealisations;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static) shared(msg1, msg2)
+ #pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
struct time t0, t1;
const int ithread = omp_get_thread_num();
@@ -624,19 +623,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)))
- {
- #pragma omp critical
- {
- if(msg1 == 0) {
- msg1 = 1,
- 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"));
- }
- }
+ && !(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;
}
@@ -657,15 +648,8 @@ XD(solve_boundary_flux)
imposed_temp = interface_side_get_temperature(interf, &frag);
if(imposed_temp >= 0) {
/* Flux computation on T boundaries is not supported yet */
- #pragma omp critical
- {
- if(msg2 == 0) {
- msg2 = 1,
- log_err(scn->dev,
- "%s: Attempt to compute a flux at a Dirichlet boundary (not available yet).\n",
- FUNC_NAME);
- }
- }
+ 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;
}
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -347,7 +347,6 @@ XD(solve_probe_boundary_flux)
int64_t irealisation = 0;
size_t i;
int progress = 0;
- int msg1 = 0;
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
@@ -391,9 +390,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,
@@ -467,7 +466,7 @@ XD(solve_probe_boundary_flux)
/* Here we go! Launch the Monte Carlo estimation */
nrealisations = args->nrealisations;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static) shared(msg1)
+ #pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
struct time t0, t1;
const int ithread = omp_get_thread_num();
@@ -505,15 +504,8 @@ XD(solve_probe_boundary_flux)
imposed_temp = interface_side_get_temperature(interf, &frag);
if(imposed_temp >= 0) {
/* Flux computation on T boundaries is not supported yet */
- #pragma omp critical
- {
- if(msg1 == 0) {
- msg1 = 1,
- log_err(scn->dev,
- "%s: Attempt to compute a flux at a Dirichlet boundary (not available yet).\n",
- FUNC_NAME);
- }
- }
+ 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;
}
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
