commit 5150f84e5648a45de70ef79534b30cfdca070552
parent 9d3a75eed9eff3a91a9ffce1f8651f7fdd79c93c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Nov 2021 10:30:11 +0100
Update the Picard order API
The Picard order was a parameter of the scene when it should have been
an argument of the solve functions. This is what this commit does.
Diffstat:
21 files changed, 584 insertions(+), 312 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -385,9 +385,6 @@ struct sdis_scene_create_args {
/* Min/max temperature used to linearise the radiative temperature */
double t_range[2];
-
- /* Picard order used to estimate the radiative temperature */
- size_t picard_order;
};
#define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \
@@ -399,8 +396,7 @@ struct sdis_scene_create_args {
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
- {0.0, -1.0}, /* Temperature range */ \
- 1 /* Picard order */ \
+ {0.0, -1.0} /* Temperature range */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
SDIS_SCENE_CREATE_ARGS_DEFAULT__;
@@ -412,6 +408,12 @@ struct sdis_solve_probe_args {
size_t nrealisations; /* #realisations */
double position[3]; /* Probe position */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
int register_paths; /* Combination of enum sdis_heat_path_flag */
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
};
@@ -419,6 +421,7 @@ struct sdis_solve_probe_args {
10000, /* #realisations */ \
{0,0,0}, /* Position */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
SDIS_HEAT_PATH_NONE, /* Register paths mask */ \
NULL /* RNG state */ \
}
@@ -431,6 +434,12 @@ struct sdis_solve_probe_boundary_args {
size_t iprim; /* Identifier of the primitive on which the probe lies */
double uv[2]; /* Parametric coordinates of the probe onto the primitve */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
enum sdis_side side; /* Side of iprim on which the probe lies */
int register_paths; /* Combination of enum sdis_heat_path_flag */
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
@@ -440,6 +449,7 @@ struct sdis_solve_probe_boundary_args {
0, /* Primitive identifier */ \
{0,0}, /* UV */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
SDIS_SIDE_NULL__, \
SDIS_HEAT_PATH_NONE, \
NULL /* RNG state */ \
@@ -454,6 +464,12 @@ struct sdis_solve_boundary_args {
const enum sdis_side* sides; /* Per primitive side to consider */
size_t nprimitives; /* #primitives */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
int register_paths; /* Combination of enum sdis_heat_path_flag */
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
};
@@ -463,6 +479,7 @@ struct sdis_solve_boundary_args {
NULL, /* Per primitive side */ \
0, /* #primitives */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
SDIS_HEAT_PATH_NONE, \
NULL /* RNG state */ \
}
@@ -473,6 +490,12 @@ struct sdis_solve_medium_args {
size_t nrealisations; /* #realisations */
struct sdis_medium* medium; /* Medium to solve */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
int register_paths; /* Combination of enum sdis_heat_path_flag */
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
};
@@ -480,6 +503,7 @@ struct sdis_solve_medium_args {
10000, /* #realisations */ \
NULL, /* Medium */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
SDIS_HEAT_PATH_NONE, \
NULL /* RNG state */ \
}
@@ -491,6 +515,12 @@ struct sdis_solve_probe_boundary_flux_args {
size_t iprim; /* Identifier of the primitive on which the probe lies */
double uv[2]; /* Parametric coordinates of the probe onto the primitve */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
};
#define SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \
@@ -498,6 +528,7 @@ struct sdis_solve_probe_boundary_flux_args {
0, /* Primitive identifier */ \
{0,0}, /* UV */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
NULL /* RNG state */ \
}
static const struct sdis_solve_probe_boundary_flux_args
@@ -509,6 +540,12 @@ struct sdis_solve_boundary_flux_args {
const size_t* primitives; /* List of boundary primitives to handle */
size_t nprimitives; /* #primitives */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
};
#define SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \
@@ -516,6 +553,7 @@ struct sdis_solve_boundary_flux_args {
NULL, /* List or primitive ids */ \
0, /* #primitives */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
NULL /* RNG state */ \
}
static const struct sdis_solve_boundary_flux_args
@@ -525,6 +563,12 @@ SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT =
struct sdis_solve_camera_args {
struct sdis_camera* cam; /* Point of view */
double time_range[2]; /* Observation time */
+
+ /* Set the Picard recursion order to estimate the radiative temperature. An
+ * order of one means that the radiative temperature is linearized, while
+ * higher orders allow the estimation of the T4 radiative transfer. */
+ size_t picard_order;
+
size_t image_resolution[2]; /* Image resolution */
size_t spp; /* #samples per pixel */
int register_paths; /* Combination of enum sdis_heat_path_flag */
@@ -532,6 +576,7 @@ struct sdis_solve_camera_args {
#define SDIS_SOLVE_CAMERA_ARGS_DEFAULT__ { \
NULL, /* Camera */ \
{DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* Picard order */ \
{512,512}, /* Image resolution */ \
256, /* #realisations per pixel */ \
SDIS_HEAT_PATH_NONE \
@@ -866,20 +911,6 @@ sdis_scene_set_temperature_range
(struct sdis_scene* scn,
const double t_range[2]);
-/* Get the picard recursion order. */
-SDIS_API res_T
-sdis_scene_get_picard_order
- (const struct sdis_scene* scn,
- size_t* picard_order);
-
-/* Set the Picard recursion order to estimate the radiative temperature. An
- * order of one means that the radiative temperature is linearized, while
- * higher orders allow the estimation of the T4 radiative transfer. */
- SDIS_API res_T
-sdis_scene_set_picard_order
- (struct sdis_scene* scn,
- const size_t picard_order);
-
/* Search the point onto the scene geometry that is the closest of `pos'. The
* `radius' parameter controls the maximum search distance around `pos'. The
* returned closest point is expressed locally to the geometric primitive onto
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -34,6 +34,12 @@ struct rwalk_context {
double That2; /* That^2 */
double That3; /* That^3 */
+ /* Maximum branchings i.e. the maximum number of times
+ * XD(compute_temperature) can be called. It controls the number of
+ * ramifications of the heat path and currently corresponds to the Picard
+ * order used to estimate the radiative temperature. */
+ size_t max_branchings;
+
/* Number of heat path branchings */
size_t nbranchings;
};
@@ -46,6 +52,7 @@ struct rwalk_context {
0, /* That */ \
0, /* That^2 */ \
0, /* That^3 */ \
+ 0, /* Max #branchings */ \
SIZE_MAX, /* #branchings */ \
}
static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -89,10 +89,10 @@ XD(boundary_path)
if(mdm_front->type == mdm_back->type) {
res = XD(solid_solid_boundary_path)(scn, ctx, &frag, rwalk, rng, T);
- } else if(ctx->nbranchings == scn->max_branchings) {
+ } else if(ctx->nbranchings == ctx->max_branchings) {
res = XD(solid_fluid_boundary_picard1_path)(scn, ctx, &frag, rwalk, rng, T);
} else {
- ASSERT(ctx->nbranchings < scn->max_branchings);
+ ASSERT(ctx->nbranchings < ctx->max_branchings);
res = XD(solid_fluid_boundary_picardN_path)(scn, ctx, &frag, rwalk, rng, T);
}
if(res != RES_OK) goto error;
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -13,8 +13,26 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#include "sdis_medium_c.h"
#include "sdis_realisation.h"
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+check_ray_realisatio_args(const struct ray_realisation_args* args)
+{
+ return args
+ && args->rng
+ && args->medium
+ && args->medium->type == SDIS_FLUID
+ && args->time >= 0
+ && args->picard_order > 0;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
/* Generate the generic realisations */
#define SDIS_XD_DIMENSION 2
#include "sdis_realisation_Xd.h"
@@ -24,12 +42,7 @@
res_T
ray_realisation_3d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- struct sdis_medium* medium,
- const double position[3],
- const double direction[3],
- const double time,
- struct sdis_heat_path* heat_path, /* May be NULL */
+ struct ray_realisation_args* args,
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -37,34 +50,35 @@ ray_realisation_3d
struct temperature_3d T = TEMPERATURE_NULL_3d;
float dir[3];
res_T res = RES_OK;
- ASSERT(scn && position && direction && time>=0 && weight);
- ASSERT(medium && medium->type == SDIS_FLUID);
+ ASSERT(scn && weight && check_ray_realisatio_args(args));
- d3_set(rwalk.vtx.P, position);
- rwalk.vtx.time = time;
+ d3_set(rwalk.vtx.P, args->position);
+ rwalk.vtx.time = args->time;
rwalk.hit = S3D_HIT_NULL;
rwalk.hit_side = SDIS_SIDE_NULL__;
- rwalk.mdm = medium;
+ rwalk.mdm = args->medium;
- ctx.heat_path = heat_path;
+ ctx.heat_path = args->heat_path;
ctx.Tmin = scn->tmin;
ctx.Tmin2 = ctx.Tmin * ctx.Tmin;
ctx.Tmin3 = ctx.Tmin * ctx.Tmin2;
ctx.That = scn->tmax;
ctx.That2 = ctx.That * ctx.That;
ctx.That3 = ctx.That * ctx.That2;
+ ctx.max_branchings = args->picard_order - 1;
- f3_set_d3(dir, direction);
+ f3_set_d3(dir, args->direction);
/* Register the starting position against the heat path */
- res = register_heat_vertex(heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE);
+ res = register_heat_vertex
+ (args->heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE);
if(res != RES_OK) goto error;
- res = trace_radiative_path_3d(scn, dir, &ctx, &rwalk, rng, &T);
+ res = trace_radiative_path_3d(scn, dir, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
if(!T.done) {
- res = compute_temperature_3d(scn, &ctx, &rwalk, rng, &T);
+ res = compute_temperature_3d(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -56,91 +56,117 @@ compute_temperature_3d
/*******************************************************************************
* Realisation at a given position and time IN a medium
******************************************************************************/
+struct probe_realisation_args {
+ struct ssp_rng* rng;
+ struct sdis_medium* medium; /* Medium into which the realisation starts */
+ double position[3]; /* Probe position */
+ double time; /* Observation time */
+ size_t picard_order; /* Picard order to estimate radiative temperature */
+ struct green_path_handle* green_path; /* May be NULL */
+ struct sdis_heat_path* heat_path; /* May be NULL */
+ size_t irealisation; /* Id of the realisation (for debug) */
+};
+#define PROBE_REALISATION_ARGS_NULL__ { \
+ NULL, NULL, {0,0,0}, -1, 0, NULL, NULL, 0 \
+}
+static const struct probe_realisation_args PROBE_REALISATION_ARGS_NULL =
+ PROBE_REALISATION_ARGS_NULL__;
+
extern LOCAL_SYM res_T
probe_realisation_2d
- (const size_t irealisation, /* For debug */
- struct sdis_scene* scn,
- struct ssp_rng* rng,
- struct sdis_medium* medium,
- const double position[2],
- const double time,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ (struct sdis_scene* scn,
+ struct probe_realisation_args* args,
double* weight);
extern LOCAL_SYM res_T
probe_realisation_3d
- (const size_t irealisation, /* For debug */
- struct sdis_scene* scn,
- struct ssp_rng* rng,
- struct sdis_medium* medium,
- const double position[3],
- const double time,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ (struct sdis_scene* scn,
+ struct probe_realisation_args* args,
double* weight);
/*******************************************************************************
* Realisation at a given position and time ON a given side of a boundary
******************************************************************************/
+struct boundary_realisation_args {
+ struct ssp_rng* rng;
+ size_t iprim; /* Index of the geometruc primitive */
+ double uv[2]; /* Parametric coordinate into the geometric primitive */
+ double time; /* Observation time */
+ size_t picard_order; /* Picard order to estimate radiative temperature */
+ enum sdis_side side; /* Side of the geometric primitive */
+ struct green_path_handle* green_path; /* May be NULL */
+ struct sdis_heat_path* heat_path; /* May be NULL */
+};
+#define BOUNDARY_REALISATION_ARGS_NULL__ { \
+ NULL, SIZE_MAX, {0,0}, -1, 0, SDIS_SIDE_NULL__, NULL, NULL \
+}
+static const struct boundary_realisation_args BOUNDARY_REALISATION_ARGS_NULL =
+ BOUNDARY_REALISATION_ARGS_NULL__;
+
extern LOCAL_SYM res_T
boundary_realisation_2d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double u[1],
- const double time,
- const enum sdis_side side,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ struct boundary_realisation_args* args,
double* weight);
extern LOCAL_SYM res_T
boundary_realisation_3d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double uv[2],
- const double time,
- const enum sdis_side side,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ struct boundary_realisation_args* args,
double* weight);
+/*******************************************************************************
+ * Realisation at a given position and time ON a given side of a boundary
+ ******************************************************************************/
+struct boundary_flux_realisation_args {
+ struct ssp_rng* rng;
+ size_t iprim; /* Index of the geometruc primitive */
+ double uv[2]; /* Parametric coordinate into the geometric primitive */
+ double time; /* Observation time */
+ size_t picard_order; /* Picard order to estimate radiative temperature */
+ enum sdis_side solid_side; /* Side of the geometric primitive */
+ int flux_mask; /* Combination of enum flux_flag */
+};
+#define BOUNDARY_FLUX_REALISATION_ARGS_NULL__ { \
+ NULL, SIZE_MAX, {0,0}, -1, 0, SDIS_SIDE_NULL__, 0 \
+}
+static const struct boundary_flux_realisation_args
+BOUNDARY_FLUX_REALISATION_ARGS_NULL = BOUNDARY_FLUX_REALISATION_ARGS_NULL__;
+
extern LOCAL_SYM res_T
boundary_flux_realisation_2d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double uv[1],
- const double time,
- const enum sdis_side solid_side,
- const int flux_mask, /* Combination of enum flux_flag */
+ struct boundary_flux_realisation_args* args,
struct bound_flux_result* result);
extern LOCAL_SYM res_T
boundary_flux_realisation_3d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double uv[2],
- const double time,
- const enum sdis_side solid_side,
- const int flux_mask, /* Combination of enum flux_flag */
+ struct boundary_flux_realisation_args* args,
struct bound_flux_result* result);
/*******************************************************************************
* Realisation along a given ray at a given time. Available only in 3D.
******************************************************************************/
+struct ray_realisation_args {
+ struct ssp_rng* rng;
+ struct sdis_medium* medium; /* Medium into which the realisation starts */
+ double position[3]; /* Ray position */
+ double direction[3]; /* Ray direction */
+ double time; /* Observation time */
+ size_t picard_order; /* Picard order to estimate radiative temperature */
+ struct sdis_heat_path* heat_path; /* May be NULL */
+};
+#define RAY_REALISATION_ARGS_NULL__ { \
+ NULL, NULL, {0,0,0}, {0,0,0}, -1, 0, NULL \
+}
+static const struct ray_realisation_args RAY_REALISATION_ARGS_NULL =
+ RAY_REALISATION_ARGS_NULL__;
+
extern LOCAL_SYM res_T
ray_realisation_3d
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- struct sdis_medium* medium,
- const double position[3],
- const double direction[3],
- const double time,
- struct sdis_heat_path* heat_path, /* May be NULL */
+ struct ray_realisation_args* args,
double* weight);
#endif /* SDIS_REALISATION_H */
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -27,6 +27,52 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Non generic helper functions
+ ******************************************************************************/
+#ifndef SDIS_REALISATION_XD_H
+#define SDIS_REALISATION_XD_H
+
+static INLINE int
+check_probe_realisation_args(const struct probe_realisation_args* args)
+{
+ return args
+ && args->rng
+ && args->medium
+ && args->time >= 0
+ && args->picard_order > 0;
+}
+
+static INLINE int
+check_boundary_realisation_args(const struct boundary_realisation_args* args)
+{
+ return args
+ && args->rng
+ && args->uv[0] >= 0
+ && args->uv[0] <= 1
+ && args->uv[1] >= 0
+ && args->uv[1] <= 1
+ && args->time >= 0
+ && args->picard_order > 0
+ && (args->side == SDIS_FRONT || args->side == SDIS_BACK);
+}
+
+static INLINE int
+check_boundary_flux_realisation_args
+ (const struct boundary_flux_realisation_args* args)
+{
+ return args
+ && args->rng
+ && args->uv[0] >= 0
+ && args->uv[0] <= 1
+ && args->uv[1] >= 0
+ && args->uv[1] <= 1
+ && args->time >= 0
+ && args->picard_order > 0
+ && (args->solid_side == SDIS_FRONT || args->solid_side == SDIS_BACK);
+}
+#endif /* SDIS_REALISATION_XD_H */
+
+/*******************************************************************************
* Local functions
******************************************************************************/
res_T
@@ -52,7 +98,7 @@ XD(compute_temperature)
ASSERT(scn && ctx && rwalk && rng && T);
ctx->nbranchings += 1;
- CHK(ctx->nbranchings <= scn->max_branchings);
+ CHK(ctx->nbranchings <= ctx->max_branchings);
if(ctx->heat_path && T->func == XD(boundary_path)) {
heat_vtx = heat_path_get_last_vertex(ctx->heat_path);
@@ -112,14 +158,8 @@ error:
res_T
XD(probe_realisation)
- (const size_t irealisation, /* For debug */
- struct sdis_scene* scn,
- struct ssp_rng* rng,
- struct sdis_medium* medium,
- const double position[DIM],
- const double time,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ (struct sdis_scene* scn,
+ struct probe_realisation_args* args,
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -131,38 +171,37 @@ XD(probe_realisation)
(const struct sdis_medium* mdm,
const struct sdis_rwalk_vertex* vtx);
res_T res = RES_OK;
- ASSERT(medium && position && weight && time >= 0);
- (void)irealisation;
+ ASSERT(scn && weight && check_probe_realisation_args(args));
- switch(medium->type) {
+ switch(args->medium->type) {
case SDIS_FLUID:
T.func = XD(convective_path);
get_initial_temperature = fluid_get_temperature;
- t0 = fluid_get_t0(medium);
+ t0 = fluid_get_t0(args->medium);
break;
case SDIS_SOLID:
T.func = XD(conductive_path);
get_initial_temperature = solid_get_temperature;
- t0 = solid_get_t0(medium);
+ t0 = solid_get_t0(args->medium);
break;
default: FATAL("Unreachable code\n"); break;
}
- dX(set)(rwalk.vtx.P, position);
- rwalk.vtx.time = time;
+ dX(set)(rwalk.vtx.P, args->position);
+ rwalk.vtx.time = args->time;
/* Register the starting position against the heat path */
- type = medium->type == SDIS_SOLID
+ type = args->medium->type == SDIS_SOLID
? SDIS_HEAT_VERTEX_CONDUCTION
: SDIS_HEAT_VERTEX_CONVECTION;
- res = register_heat_vertex(heat_path, &rwalk.vtx, 0, type);
+ res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0, type);
if(res != RES_OK) goto error;
if(t0 >= rwalk.vtx.time) {
double tmp;
/* Check the initial condition. */
rwalk.vtx.time = t0;
- tmp = get_initial_temperature(medium, &rwalk.vtx);
+ tmp = get_initial_temperature(args->medium, &rwalk.vtx);
if(tmp >= 0) {
*weight = tmp;
goto exit;
@@ -177,18 +216,19 @@ XD(probe_realisation)
}
rwalk.hit = SXD_HIT_NULL;
- rwalk.mdm = medium;
+ rwalk.mdm = args->medium;
- ctx.green_path = green_path;
- ctx.heat_path = heat_path;
+ ctx.green_path = args->green_path;
+ ctx.heat_path = args->heat_path;
ctx.Tmin = scn->tmin;
ctx.Tmin2 = ctx.Tmin * ctx.Tmin;
ctx.Tmin3 = ctx.Tmin * ctx.Tmin2;
ctx.That = scn->tmax;
ctx.That2 = ctx.That * ctx.That;
ctx.That3 = ctx.That * ctx.That2;
+ ctx.max_branchings = args->picard_order - 1;
- res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
+ res = XD(compute_temperature)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
ASSERT(T.value >= 0);
@@ -203,13 +243,7 @@ error:
res_T
XD(boundary_realisation)
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double uv[DIM-1],
- const double time,
- const enum sdis_side side,
- struct green_path_handle* green_path, /* May be NULL */
- struct sdis_heat_path* heat_path, /* May be NULL */
+ struct boundary_realisation_args* args,
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -222,23 +256,23 @@ XD(boundary_realisation)
float st[2];
#endif
res_T res = RES_OK;
- ASSERT(uv && weight && time >= 0);
+ ASSERT(scn && weight && check_boundary_realisation_args(args));
T.func = XD(boundary_path);
- rwalk.hit_side = side;
+ rwalk.hit_side = args->side;
rwalk.hit.distance = 0;
- rwalk.vtx.time = time;
+ rwalk.vtx.time = args->time;
rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */
#if SDIS_XD_DIMENSION == 2
- st = (float)uv[0];
+ st = (float)args->uv[0];
#else
- f2_set_d2(st, uv);
+ f2_set_d2(st, args->uv);
#endif
/* Fetch the primitive */
SXD(scene_view_get_primitive
- (scn->sXd(view), (unsigned int)iprim, &rwalk.hit.prim));
+ (scn->sXd(view), (unsigned int)args->iprim, &rwalk.hit.prim));
/* Retrieve the world space position of the probe onto the primitive */
SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_POSITION, st, &attr));
@@ -254,20 +288,21 @@ XD(boundary_realisation)
f2_set(rwalk.hit.uv, st);
#endif
- res = register_heat_vertex(heat_path, &rwalk.vtx, 0/*weight*/,
+ res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0/*weight*/,
SDIS_HEAT_VERTEX_CONDUCTION);
if(res != RES_OK) goto error;
- ctx.green_path = green_path;
- ctx.heat_path = heat_path;
+ ctx.green_path = args->green_path;
+ ctx.heat_path = args->heat_path;
ctx.Tmin = scn->tmin;
ctx.Tmin2 = ctx.Tmin * ctx.Tmin;
ctx.Tmin3 = ctx.Tmin * ctx.Tmin2;
ctx.That = scn->tmax;
ctx.That2 = ctx.That * ctx.That;
ctx.That3 = ctx.That * ctx.That2;
+ ctx.max_branchings = args->picard_order - 1;
- res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
+ res = XD(compute_temperature)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
*weight = T.value;
@@ -281,12 +316,7 @@ error:
res_T
XD(boundary_flux_realisation)
(struct sdis_scene* scn,
- struct ssp_rng* rng,
- const size_t iprim,
- const double uv[DIM-1],
- const double time,
- const enum sdis_side solid_side,
- const int flux_mask,
+ struct boundary_flux_realisation_args* args,
struct bound_flux_result* result)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -304,29 +334,36 @@ XD(boundary_flux_realisation)
#endif
double P[SDIS_XD_DIMENSION];
float N[SDIS_XD_DIMENSION];
- const double Tmin = scn->tmin;
- const double Tmin2 = Tmin * Tmin;
- const double Tmin3 = Tmin * Tmin2;
- const double That = scn->tmax;
- const double That2 = That * That;
- const double That3 = That * That2;
- const enum sdis_side fluid_side =
- (solid_side == SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT;
+ double Tmin, Tmin2, Tmin3;
+ double That, That2, That3;
+ enum sdis_side fluid_side;
res_T res = RES_OK;
- const char compute_radiative = (flux_mask & FLUX_FLAG_RADIATIVE) != 0;
- const char compute_convective = (flux_mask & FLUX_FLAG_CONVECTIVE) != 0;
- ASSERT(uv && result && time >= 0 );
+ char compute_radiative;
+ char compute_convective;
+ ASSERT(scn && result && check_boundary_flux_realisation_args(args));
#if SDIS_XD_DIMENSION == 2
#define SET_PARAM(Dest, Src) (Dest).u = (Src);
- st = (float)uv[0];
+ st = (float)args->uv[0];
#else
#define SET_PARAM(Dest, Src) f2_set((Dest).uv, (Src));
- f2_set_d2(st, uv);
+ f2_set_d2(st, args->uv);
#endif
+ Tmin = scn->tmin;
+ Tmin2 = Tmin * Tmin;
+ Tmin3 = Tmin * Tmin2;
+ That = scn->tmax;
+ That2 = That * That;
+ That3 = That * That2;
+
+ fluid_side = (args->solid_side/*solid*/==SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT;
+
+ compute_radiative = (args->flux_mask & FLUX_FLAG_RADIATIVE) != 0;
+ compute_convective = (args->flux_mask & FLUX_FLAG_CONVECTIVE) != 0;
+
/* Fetch the primitive */
- SXD(scene_view_get_primitive(scn->sXd(view), (unsigned int)iprim, &prim));
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned int)args->iprim, &prim));
/* Retrieve the world space position of the probe onto the primitive */
SXD(primitive_get_attrib(&prim, SXD_POSITION, st, &attr));
@@ -340,7 +377,7 @@ XD(boundary_flux_realisation)
rwalk = XD(RWALK_NULL); \
rwalk.hit_side = (Side); \
rwalk.hit.distance = 0; \
- rwalk.vtx.time = time; \
+ rwalk.vtx.time = args->time; \
rwalk.mdm = (Mdm); \
rwalk.hit.prim = prim; \
SET_PARAM(rwalk.hit, st); \
@@ -349,27 +386,28 @@ XD(boundary_flux_realisation)
ctx.That = That; \
ctx.That2 = That2; \
ctx.That3 = That3; \
+ ctx.max_branchings = args->picard_order - 1; \
dX(set)(rwalk.vtx.P, P); \
fX(set)(rwalk.hit.normal, N); \
T = XD(TEMPERATURE_NULL); \
} (void)0
/* Compute boundary temperature */
- RESET_WALK(solid_side, NULL);
+ RESET_WALK(args->solid_side, NULL);
T.func = XD(boundary_path);
- res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
+ res = XD(compute_temperature)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
result->Tboundary = T.value;
/* Fetch the fluid medium */
- interf = scene_get_interface(scn, (unsigned)iprim);
+ interf = scene_get_interface(scn, (unsigned)args->iprim);
fluid_mdm = interface_get_medium(interf, fluid_side);
/* Compute radiative temperature */
if(compute_radiative) {
RESET_WALK(fluid_side, fluid_mdm);
T.func = XD(radiative_path);
- res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
+ res = XD(compute_temperature)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
ASSERT(T.value >= 0);
result->Tradiative = T.value;
@@ -379,7 +417,7 @@ XD(boundary_flux_realisation)
if(compute_convective) {
RESET_WALK(fluid_side, fluid_mdm);
T.func = XD(convective_path);
- res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
+ res = XD(compute_temperature)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
result->Tfluid = T.value;
}
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -243,26 +243,6 @@ sdis_scene_set_temperature_range
}
res_T
-sdis_scene_get_picard_order
- (const struct sdis_scene* scn,
- size_t* picard_order)
-{
- if(!scn || !picard_order) return RES_BAD_ARG;
- *picard_order = scene_get_picard_order(scn);
- return RES_OK;
-}
-
-res_T
-sdis_scene_set_picard_order
- (struct sdis_scene* scn,
- const size_t picard_order)
-{
- if(!scn || picard_order < 1) return RES_BAD_ARG;
- scn->max_branchings = picard_order-1;
- return RES_OK;
-}
-
-res_T
sdis_scene_find_closest_point
(const struct sdis_scene* scn,
const double pos[],
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -120,8 +120,7 @@ check_sdis_scene_create_args(const struct sdis_scene_create_args* args)
&& args->nprimitives < UINT_MAX
&& args->nvertices
&& args->nvertices < UINT_MAX
- && args->fp_to_meter > 0
- && args->picard_order > 0;
+ && args->fp_to_meter > 0;
}
#endif /* SDIS_SCENE_XD_H */
@@ -917,7 +916,6 @@ XD(scene_create)
scn->tmin = args->t_range[0];
scn->tmax = args->t_range[1];
scn->outer_enclosure_id = UINT_MAX;
- scn->max_branchings = args->picard_order - 1;
darray_interf_init(dev->allocator, &scn->interfaces);
darray_medium_init(dev->allocator, &scn->media);
darray_prim_prop_init(dev->allocator, &scn->prim_props);
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -213,12 +213,6 @@ struct sdis_scene {
double tmin; /* Minimum temperature of the system (In Kelvin) */
double tmax; /* Maximum temperature of the system (In Kelvin) */
- /* Maximum branchings i.e. the maximum number of times
- * XD(compute_temperature) can be called. It controls the number of
- * ramifications of the heat path and currently corresponds to the Picard
- * order used to estimate the radiative temperature. */
- size_t max_branchings;
-
ref_T ref;
struct sdis_device* dev;
};
@@ -303,13 +297,5 @@ scene_is_2d(const struct sdis_scene* scn)
return scn->s2d_view != NULL;
}
-static FINLINE size_t
-scene_get_picard_order(const struct sdis_scene* scn)
-{
- ASSERT(scn);
- return scn->max_branchings+1;
-}
-
-
#endif /* SDIS_SCENE_C_H */
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -72,6 +72,7 @@ solve_pixel
const size_t nrealisations,
const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
+ const size_t picard_order,
struct sdis_estimator* estimator)
{
struct accum acc_temp = ACCUM_NULL;
@@ -84,6 +85,7 @@ solve_pixel
ASSERT(estimator && time_range);
FOR_EACH(irealisation, 0, nrealisations) {
+ struct ray_realisation_args realis_args = RAY_REALISATION_ARGS_NULL;
struct time t0, t1;
double samp[2]; /* Pixel sample */
double ray_pos[3];
@@ -111,8 +113,14 @@ solve_pixel
camera_ray(cam, samp, ray_pos, ray_dir);
/* Launch the realisation */
- res_simul = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
- time, pheat_path, &w);
+ realis_args.rng = rng;
+ realis_args.medium = mdm;
+ realis_args.time = time;
+ realis_args.picard_order = picard_order;
+ realis_args.heat_path = pheat_path;
+ d3_set(realis_args.position, ray_pos);
+ d3_set(realis_args.direction, ray_dir);
+ res_simul = ray_realisation_3d(scn, &realis_args, &w);
/* Handle fatal error */
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
@@ -172,6 +180,7 @@ solve_tile
const size_t spp, /* #samples per pixel */
const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
+ const size_t picard_order,
struct sdis_estimator_buffer* buf)
{
size_t mcode; /* Morton code of the tile pixel */
@@ -201,7 +210,7 @@ solve_tile
estimator = estimator_buffer_grab(buf, ipix[0], ipix[1]);
res = solve_pixel(scn, rng, mdm, cam, time_range,
- ipix, spp, register_paths, pix_sz, estimator);
+ ipix, spp, register_paths, pix_sz, picard_order, estimator);
if(res != RES_OK) goto error;
}
@@ -439,7 +448,8 @@ sdis_solve_camera
/* Draw the tile */
res_local = solve_tile(scn, rng, medium, args->cam, args->time_range,
- tile_org, tile_sz, args->spp, args->register_paths, pix_sz, buf);
+ tile_org, tile_sz, args->spp, args->register_paths, pix_sz,
+ args->picard_order, buf);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -1,5 +1,5 @@
/* Copyright (C) 2016-2021 |Meso|Star> (contact@meso-star.com)
- *
+*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -39,6 +39,43 @@ static const struct XD(boundary_context) XD(BOUNDARY_CONTEXT_NULL) = {
/*******************************************************************************
* Help functions
******************************************************************************/
+#ifndef SDIS_SOLVE_BOUNDARY_XD
+#define SDIS_SOLVE_BOUNDARY_XD
+
+static INLINE res_T
+check_solve_boundary_args(const struct sdis_solve_boundary_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Check #realisations */
+ if(!args->nrealisations || args->nrealisations > INT64_MAX) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check the list of primitives */
+ if(!args->primitives || !args->sides || !args->nprimitives) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check time range */
+ if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
+ return RES_BAD_ARG;
+ }
+ if(args->time_range[1] > DBL_MAX
+ && args->time_range[0] != args->time_range[1]) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check picard order */
+ if(args->picard_order < 1) {
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+#endif /* SDIS_SOLVE_BOUNDARY_XD */
+
static INLINE void
XD(boundary_get_indices)(const unsigned iprim, unsigned ids[DIM], void* context)
{
@@ -101,29 +138,24 @@ XD(solve_boundary)
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
- if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX
- || !args->primitives || !args->sides || !args->nprimitives) {
+ if(!scn) {
res = RES_BAD_ARG;
goto error;
}
+
+ res = check_solve_boundary_args(args);
+ if(res != RES_OK) goto error;
+
if(!out_estimator && !out_green) {
res = RES_BAD_ARG;
goto error;
}
- if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(args->time_range[1] > DBL_MAX
- && args->time_range[0] != args->time_range[1]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(out_green && scene_get_picard_order(scn) != 1) {
+
+ if(out_green && args->picard_order != 1) {
log_err(scn->dev, "%s: the evaluation of the green function does not make "
"sense when dealing with the non-linearities of the system; i.e. picard "
"order must be set to 1 while it is currently set to %lu.\n",
- FUNC_NAME, (unsigned long)scene_get_picard_order(scn));
+ FUNC_NAME, (unsigned long)args->picard_order);
res = RES_BAD_ARG;
goto error;
}
@@ -237,6 +269,7 @@ XD(solve_boundary)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) {
+ struct boundary_realisation_args realis_args = BOUNDARY_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -306,8 +339,18 @@ XD(solve_boundary)
side = args->sides[prim.prim_id];
/* Invoke the boundary realisation */
- res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
- pgreen_path, pheat_path, &w);
+ realis_args.rng = rng;
+ realis_args.iprim = iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.side = side;
+ realis_args.green_path = pgreen_path;
+ realis_args.heat_path = pheat_path;
+ realis_args.uv[0] = uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = uv[1];
+#endif
+ res_simul = XD(boundary_realisation)(scn, &realis_args, &w);
/* Fatal error */
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
@@ -566,6 +609,8 @@ XD(solve_boundary_flux)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct boundary_flux_realisation_args realis_args =
+ BOUNDARY_FLUX_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -663,8 +708,19 @@ XD(solve_boundary_flux)
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, flux_mask, &result);
+
+ /* Invoke the boundary flux realisation */
+ realis_args.rng = rng;
+ realis_args.iprim = iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.solid_side = solid_side;
+ realis_args.flux_mask = flux_mask;
+ realis_args.uv[0] = uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = uv[1];
+#endif
+ res_simul = XD(boundary_flux_realisation)(scn, &realis_args, &result);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -135,6 +135,36 @@ sample_medium_enclosure
return enc_cumul_found->enc;
}
+static INLINE res_T
+check_solve_medium_args(const struct sdis_solve_medium_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Check the medium */
+ if(!args->medium) return RES_BAD_ARG;
+
+ /* Check #realisations */
+ if(!args->nrealisations || args->nrealisations > INT64_MAX) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check time range */
+ if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
+ return RES_BAD_ARG;
+ }
+ if(args->time_range[1] > DBL_MAX
+ && args->time_range[0] != args->time_range[1]) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check picard order */
+ if(args->picard_order < 1) {
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
#endif /* !SDIS_SOLVE_MEDIUM_XD_H */
/*******************************************************************************
@@ -221,30 +251,24 @@ XD(solve_medium)
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
- if(!scn || !args || !args->medium || !args->nrealisations
- || args->nrealisations > INT64_MAX) {
+ if(!scn) {
res = RES_BAD_ARG;
goto error;
}
+
+ res = check_solve_medium_args(args);
+ if(res != RES_OK) goto error;
+
if(!out_estimator && !out_green) {
res = RES_BAD_ARG;
goto error;
}
- if(out_estimator) {
- if(args->time_range[0] < 0
- || args->time_range[0] > args->time_range[1]
- || ( args->time_range[1] > DBL_MAX
- && args->time_range[0] != args->time_range[1])) {
- res = RES_BAD_ARG;
- goto error;
- }
- }
- if(out_green && scene_get_picard_order(scn) != 1) {
+ if(out_green && args->picard_order != 1) {
log_err(scn->dev, "%s: the evaluation of the green function does not make "
"sense when dealing with the non-linearities of the system; i.e. picard "
"order must be set to 1 while it is currently set to %lu.\n",
- FUNC_NAME, (unsigned long)scene_get_picard_order(scn));
+ FUNC_NAME, (unsigned long)args->picard_order);
res = RES_BAD_ARG;
goto error;
}
@@ -309,6 +333,7 @@ XD(solve_medium)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct probe_realisation_args realis_args = PROBE_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -330,7 +355,7 @@ XD(solve_medium)
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
time_current(&t0);
-
+
time = sample_time(rng, args->time_range);
if(out_green) {
res_local = green_function_create_path(greens[ithread], &green_path);
@@ -357,9 +382,15 @@ XD(solve_medium)
}
/* Run a probe realisation */
- res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng,
- args->medium, pos, time, pgreen_path, pheat_path, &weight);
-
+ realis_args.rng = rng;
+ realis_args.medium = args->medium;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.green_path = pgreen_path;
+ realis_args.heat_path = pheat_path;
+ realis_args.irealisation = (size_t)irealisation;
+ dX(set)(realis_args.position, pos);
+ res_simul = XD(probe_realisation)(scn, &realis_args, &weight);
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
ATOMIC_SET(&res, res_simul);
goto error_it;
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -28,6 +28,41 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+#ifndef SDIS_SOLVE_PROBE_XD_H
+#define SDIS_SOLVE_PROBE_XD_H
+
+static INLINE res_T
+check_solve_probe_args(const struct sdis_solve_probe_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Check #realisations */
+ if(!args->nrealisations || args->nrealisations > INT64_MAX) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check time range */
+ if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
+ return RES_BAD_ARG;
+ }
+ if(args->time_range[1] > DBL_MAX
+ && args->time_range[0] != args->time_range[1]) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check picard order */
+ if(args->picard_order < 1) {
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+#endif /* SDIS_SOLVE_PROBE_XD_H */
+
+/*******************************************************************************
* Generic solve function
******************************************************************************/
static res_T
@@ -53,28 +88,24 @@ XD(solve_probe)
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
- if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX) {
+ if(!scn) {
res = RES_BAD_ARG;
goto error;
}
+
+ res = check_solve_probe_args(args);
+ if(res != RES_OK) goto error;
+
if(!out_estimator && !out_green) {
res = RES_BAD_ARG;
goto error;
}
- if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(args->time_range[1] > DBL_MAX
- && args->time_range[0] != args->time_range[1]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(out_green && scene_get_picard_order(scn) != 1) {
+
+ if(out_green && args->picard_order != 1) {
log_err(scn->dev, "%s: the evaluation of the green function does not make "
"sense when dealing with the non-linearities of the system; i.e. picard "
"order must be set to 1 while it is currently set to %lu.\n",
- FUNC_NAME, (unsigned long)scene_get_picard_order(scn));
+ FUNC_NAME, (unsigned long)args->picard_order);
res = RES_BAD_ARG;
goto error;
}
@@ -138,6 +169,7 @@ XD(solve_probe)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct probe_realisation_args realis_args = PROBE_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -174,8 +206,16 @@ XD(solve_probe)
pheat_path = &heat_path;
}
- res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, medium,
- args->position, time, pgreen_path, pheat_path, &w);
+ /* Invoke the probe realisation */
+ realis_args.rng = rng;
+ realis_args.medium = medium;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.green_path = pgreen_path;
+ realis_args.heat_path = pheat_path;
+ realis_args.irealisation = (size_t)irealisation;
+ dX(set)(realis_args.position, args->position);
+ res_simul = XD(probe_realisation)(scn, &realis_args, &w);
/* Handle fatal error */
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -28,6 +28,47 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+#ifndef SDIS_SOLVE_PROBE_BOUNDARY_XD_H
+#define SDIS_SOLVE_PROBE_BOUNDARY_XD_H
+
+static INLINE res_T
+check_solve_probe_boundary_args
+ (const struct sdis_solve_probe_boundary_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Check #realisations */
+ if(!args->nrealisations || args->nrealisations > INT64_MAX) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check side */
+ if((unsigned)args->side >= SDIS_SIDE_NULL__) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check time range */
+ if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
+ return RES_BAD_ARG;
+ }
+ if(args->time_range[1] > DBL_MAX
+ && args->time_range[0] != args->time_range[1]) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check picard order */
+ if(args->picard_order < 1) {
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+#endif /* SDIS_SOLVE_PROBE_BOUNDARY_XD_H */
+
+/*******************************************************************************
* Local functions
******************************************************************************/
static res_T
@@ -52,29 +93,24 @@ XD(solve_probe_boundary)
ATOMIC nsolved_realisations = 0;
ATOMIC res = RES_OK;
- if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX
- || ((unsigned)args->side >= SDIS_SIDE_NULL__)) {
+ if(!scn) {
res = RES_BAD_ARG;
goto error;
}
+
+ res = check_solve_probe_boundary_args(args);
+ if(res != RES_OK) goto error;
+
if(!out_estimator && !out_green) {
res = RES_BAD_ARG;
goto error;
}
- if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(args->time_range[1] > DBL_MAX
- && args->time_range[0] != args->time_range[1]) {
- res = RES_BAD_ARG;
- goto error;
- }
- if(out_green && scene_get_picard_order(scn) != 1) {
+
+ if(out_green && args->picard_order != 1) {
log_err(scn->dev, "%s: the evaluation of the green function does not make "
"sense when dealing with the non-linearities of the system; i.e. picard "
"order must be set to 1 while it is currently set to %lu.\n",
- FUNC_NAME, (unsigned long)scene_get_picard_order(scn));
+ FUNC_NAME, (unsigned long)args->picard_order);
res = RES_BAD_ARG;
goto error;
}
@@ -175,6 +211,7 @@ XD(solve_probe_boundary)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct boundary_realisation_args realis_args = BOUNDARY_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -211,8 +248,19 @@ XD(solve_probe_boundary)
pheat_path = &heat_path;
}
- res_simul = XD(boundary_realisation)(scn, rng, args->iprim, args->uv, time,
- args->side, pgreen_path, pheat_path, &w);
+ /* Invoke the boundary realisation */
+ realis_args.rng = rng;
+ realis_args.iprim = args->iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.side = args->side;
+ realis_args.green_path = pgreen_path;
+ realis_args.heat_path = pheat_path;
+ realis_args.uv[0] = args->uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = args->uv[1];
+#endif
+ res_simul = XD(boundary_realisation)(scn, &realis_args, &w);
/* Handle fatal error */
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
@@ -476,6 +524,8 @@ XD(solve_probe_boundary_flux)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct boundary_flux_realisation_args realis_args =
+ BOUNDARY_FLUX_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -522,8 +572,19 @@ XD(solve_probe_boundary_flux)
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, flux_mask, &result);
+
+ /* Invoke the boundary flux realisation */
+ realis_args.rng = rng;
+ realis_args.iprim = args->iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.solid_side = solid_side;
+ realis_args.flux_mask = flux_mask;
+ realis_args.uv[0] = args->uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = args->uv[1];
+#endif
+ res_simul = XD(boundary_flux_realisation)(scn, &realis_args, &result);
/* Stop time registration */
time_sub(&t0, time_current(&t1), &t0);
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -294,34 +294,35 @@ test_invalidity_picardN_green
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
struct sdis_green_function* green = NULL;
- size_t picard_order;
CHK(scn);
- OK(sdis_scene_get_picard_order(scn, &picard_order));
- CHK(picard_order == 1);
-
- OK(sdis_scene_set_picard_order(scn, 2));
+ CHK(probe.picard_order == 1);
+ CHK(probe_bound.picard_order == 1);
+ CHK(bound.picard_order == 1);
+ CHK(mdm.picard_order == 1);
probe.position[0] = 0;
probe.position[1] = 0;
probe.position[2] = 0;
+ probe.picard_order = 2;
BA(sdis_solve_probe_green_function(scn, &probe, &green));
probe_bound.iprim = 2; /* Solid left */
probe_bound.uv[0] = 0.3;
probe_bound.uv[1] = 0.3;
probe_bound.side = SDIS_FRONT;
+ probe_bound.picard_order = 2;
BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green));
bound.primitives = &probe_bound.iprim;
bound.sides = &probe_bound.side;
bound.nprimitives = 1;
+ bound.picard_order = 2;
BA(sdis_solve_boundary_green_function(scn, &bound, &green));
mdm.medium = solid;
+ mdm.picard_order = 2;
BA(sdis_solve_medium_green_function(scn, &mdm, &green));
-
- OK(sdis_scene_set_picard_order(scn, picard_order));
}
/*******************************************************************************
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -301,32 +301,33 @@ test_invalidity_picardN_green
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
struct sdis_green_function* green = NULL;
- size_t picard_order;
CHK(scn);
- OK(sdis_scene_get_picard_order(scn, &picard_order));
- CHK(picard_order == 1);
-
- OK(sdis_scene_set_picard_order(scn, 2));
+ CHK(probe.picard_order == 1);
+ CHK(probe_bound.picard_order == 1);
+ CHK(bound.picard_order == 1);
+ CHK(mdm.picard_order == 1);
probe.position[0] = 0;
probe.position[1] = 0;
+ probe.picard_order = 2;
BA(sdis_solve_probe_green_function(scn, &probe, &green));
probe_bound.iprim = 1; /* Solid left */
probe_bound.uv[0] = 0.5;
probe_bound.side = SDIS_FRONT;
+ probe_bound.picard_order = 2;
BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green));
bound.primitives = &probe_bound.iprim;
bound.sides = &probe_bound.side;
bound.nprimitives = 1;
+ bound.picard_order = 2;
BA(sdis_solve_boundary_green_function(scn, &bound, &green));
mdm.medium = solid;
+ mdm.picard_order = 2;
BA(sdis_solve_medium_green_function(scn, &mdm, &green));
-
- OK(sdis_scene_set_picard_order(scn, picard_order));
}
/*******************************************************************************
diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c
@@ -397,7 +397,10 @@ struct reference_result {
static const struct reference_result REFERENCE_RESULT_NULL = {0,0,0};
static void
-test_picard1(struct sdis_scene* scn, const struct reference_result* ref)
+test_picard
+ (struct sdis_scene* scn,
+ const size_t picard_order,
+ const struct reference_result* ref)
{
struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
struct sdis_solve_boundary_args bound_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
@@ -406,7 +409,7 @@ test_picard1(struct sdis_scene* scn, const struct reference_result* ref)
enum sdis_scene_dimension dim;
size_t prims[2];
enum sdis_side sides[2];
- CHK(scn);
+ CHK(scn && ref && picard_order >= 1);
OK(sdis_scene_get_dimension(scn, &dim));
switch(dim) {
@@ -419,6 +422,7 @@ test_picard1(struct sdis_scene* scn, const struct reference_result* ref)
probe_args.position[0] = 0.05;
probe_args.position[1] = 0;
probe_args.position[2] = 0;
+ probe_args.picard_order = picard_order;
OK(sdis_solve_probe(scn, &probe_args, &estimator));
OK(sdis_estimator_get_temperature(estimator, &mc));
printf("Temperature at `%g %g %g' = %g ~ %g +/- %g\n",
@@ -441,6 +445,7 @@ test_picard1(struct sdis_scene* scn, const struct reference_result* ref)
bound_args.nrealisations = N;
bound_args.primitives = prims;
bound_args.sides = sides;
+ bound_args.picard_order = picard_order;
OK(sdis_solve_boundary(scn, &bound_args, &estimator));
OK(sdis_estimator_get_temperature(estimator, &mc));
printf("T1 = %g ~ %g +/- %g\n", ref->T1, mc.E, mc.SE);
@@ -462,6 +467,7 @@ test_picard1(struct sdis_scene* scn, const struct reference_result* ref)
bound_args.nrealisations = N;
bound_args.primitives = prims;
bound_args.sides = sides;
+ bound_args.picard_order = picard_order;
OK(sdis_solve_boundary(scn, &bound_args, &estimator));
OK(sdis_estimator_get_temperature(estimator, &mc));
printf("T2 = %g ~ %g +/- %g\n", ref->T2, mc.E, mc.SE);
@@ -659,8 +665,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_mX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
+ test_picard(scn_2d, 1/*Picard order*/, &ref);
+ test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
/* Test picard1 using T4 as a reference */
@@ -672,8 +678,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2;
pinterf_props[BOUNDARY_mX]->Tref = 280;
pinterf_props[BOUNDARY_pX]->Tref = 350;
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
+ test_picard(scn_2d, 1/*Picard order*/, &ref);
+ test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
/* Test picard2 */
@@ -685,20 +691,14 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_mX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- OK(sdis_scene_set_picard_order(scn_2d, 2));
- OK(sdis_scene_set_picard_order(scn_3d, 2));
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
- OK(sdis_scene_set_picard_order(scn_2d, 1));
- OK(sdis_scene_set_picard_order(scn_3d, 1));
+ test_picard(scn_2d, 2/*Picard order*/, &ref);
+ test_picard(scn_3d, 2/*Picard order*/, &ref);
printf("\n");
t_range[0] = 200;
t_range[1] = 500;
OK(sdis_scene_set_temperature_range(scn_2d, t_range));
OK(sdis_scene_set_temperature_range(scn_3d, t_range));
- OK(sdis_scene_set_picard_order(scn_2d, 3));
- OK(sdis_scene_set_picard_order(scn_3d, 3));
/* Test picard2 */
printf("Test Picard3 with a delta T of 300K\n");
@@ -711,16 +711,14 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = 450;
pinterf_props[BOUNDARY_mX]->Tref = pinterf_props[BOUNDARY_mX]->temperature;
pinterf_props[BOUNDARY_pX]->Tref = pinterf_props[BOUNDARY_pX]->temperature;
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
+ test_picard(scn_2d, 3/*Picard order*/, &ref);
+ test_picard(scn_3d, 3/*Picard order*/, &ref);
printf("\n");
t_range[0] = 280;
t_range[1] = 350;
OK(sdis_scene_set_temperature_range(scn_2d, t_range));
OK(sdis_scene_set_temperature_range(scn_3d, t_range));
- OK(sdis_scene_set_picard_order(scn_2d, 1));
- OK(sdis_scene_set_picard_order(scn_3d, 1));
pinterf_props[BOUNDARY_mX]->temperature = t_range[0];
pinterf_props[BOUNDARY_pX]->temperature = t_range[1];
@@ -737,8 +735,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_mX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
+ test_picard(scn_2d, 1/*Picard order*/, &ref);
+ test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
/* Test picard1 with a volumic power and T4 a the reference */
@@ -750,8 +748,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2;
pinterf_props[BOUNDARY_mX]->Tref = 280;
pinterf_props[BOUNDARY_pX]->Tref = 350;
- test_picard1(scn_2d, &ref);
- test_picard1(scn_3d, &ref);
+ test_picard(scn_2d, 1/*Picard order*/, &ref);
+ test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
/* Release memory */
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -139,9 +139,6 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
scn_args.fp_to_meter = 0;
BA(sdis_scene_create(dev, &scn_args, &scn));
scn_args.fp_to_meter = 1;
- scn_args.picard_order = 0;
- BA(sdis_scene_create(dev, &scn_args, &scn));
- scn_args.picard_order = 1;
/* Duplicated vertex */
ctx.positions = duplicated_vertices;
ctx.indices = dup_vrtx_indices;
@@ -274,7 +271,6 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
size_t i;
size_t iprim;
size_t dup_vrtx_indices[] = { 0, 1 };
- size_t picard_order;
enum sdis_scene_dimension dim;
ctx.positions = square_vertices;
@@ -309,9 +305,6 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
scn_args.fp_to_meter = 0;
BA(sdis_scene_2d_create(dev, &scn_args, &scn));
scn_args.fp_to_meter = 1;
- scn_args.picard_order = 0;
- BA(sdis_scene_create(dev, &scn_args, &scn));
- scn_args.picard_order = 1;
/* Duplicated vertex */
ctx.positions = duplicated_vertices;
ctx.indices = dup_vrtx_indices;
@@ -396,17 +389,6 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
CHK(t_range[0] == 1);
CHK(t_range[1] == 100);
- BA(sdis_scene_get_picard_order(NULL, &picard_order));
- BA(sdis_scene_get_picard_order(scn, NULL));
- OK(sdis_scene_get_picard_order(scn, &picard_order));
- CHK(picard_order == SDIS_SCENE_CREATE_ARGS_DEFAULT.picard_order);
- CHK(picard_order == 1);
- OK(sdis_scene_set_picard_order(scn, 3));
- OK(sdis_scene_get_picard_order(scn, &picard_order));
- CHK(picard_order == 3);
- BA(sdis_scene_set_picard_order(scn, 0));
- OK(sdis_scene_set_picard_order(scn, 1));
-
BA(sdis_scene_get_boundary_position(NULL, 1, &u0, pos));
BA(sdis_scene_get_boundary_position(scn, 4, &u0, pos));
BA(sdis_scene_get_boundary_position(scn, 1, NULL, pos));
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -333,6 +333,9 @@ main(int argc, char** argv)
probe_args.time_range[1] = 0;
BA(SOLVE(box_scn, &probe_args, &estimator));
probe_args.time_range[0] = probe_args.time_range[1] = INF;
+ probe_args.picard_order = 0;
+ BA(SOLVE(box_scn, &probe_args, &estimator));
+ probe_args.picard_order = 1;
OK(SOLVE(box_scn, &probe_args, &estimator));
OK(sdis_scene_get_boundary_position
@@ -464,6 +467,9 @@ main(int argc, char** argv)
bound_args.time_range[1] = 0;
BA(SOLVE(box_scn, &bound_args, &estimator));
bound_args.time_range[0] = bound_args.time_range[1] = INF;
+ bound_args.picard_order = 0;
+ BA(SOLVE(box_scn, &bound_args, &estimator));
+ bound_args.picard_order = 1;
/* Average temperature on the right side of the box */
OK(SOLVE(box_scn, &bound_args, &estimator));
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -453,6 +453,9 @@ main(int argc, char** argv)
solve_args.medium = solid1;
BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
solve_args.medium = solid0;
+ solve_args.picard_order = 0;
+ BA(sdis_solve_medium_green_function(scn, &solve_args, &green));
+ solve_args.picard_order = 1;
OK(sdis_solve_medium_green_function(scn, &solve_args, &green));
OK(sdis_green_function_solve(green, &estimator2));
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -380,6 +380,9 @@ main(int argc, char** argv)
solve_args.time_range[1] = 0;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
solve_args.time_range[0] = solve_args.time_range[1] = INF;
+ solve_args.picard_order = 0;
+ BA(sdis_solve_probe(scn, &solve_args, &estimator));
+ solve_args.picard_order = 1;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
BA(sdis_estimator_get_type(estimator, NULL));
