commit 9d3a75eed9eff3a91a9ffce1f8651f7fdd79c93c
parent 312238434d4392b2dbba78bf75ce48f6cff805d6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 8 Nov 2021 14:38:35 +0100
Merge branch 'feature_picardN' into develop
Diffstat:
33 files changed, 1563 insertions(+), 850 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -192,7 +192,7 @@ if(NOT NO_TEST)
new_test(test_sdis_flux)
new_test(test_sdis_interface)
new_test(test_sdis_medium)
- new_test(test_sdis_picard1)
+ new_test(test_sdis_picard)
new_test(test_sdis_scene)
new_test(test_sdis_solid_random_walk_robustness)
new_test(test_sdis_solve_camera)
diff --git a/src/sdis.h b/src/sdis.h
@@ -362,7 +362,7 @@ typedef void
struct sdis_ambient_radiative_temperature {
double temperature; /* In Kelvin */
- double reference; /* Used to linearise the radiative transfert */
+ double reference; /* Used to linearise the radiative transfer */
};
#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ {-1, -1}
static const struct sdis_ambient_radiative_temperature
@@ -382,7 +382,12 @@ struct sdis_scene_create_args {
size_t nvertices; /* #vertices */
double fp_to_meter; /* Scale factor used to convert 1.0 in 1 meter */
struct sdis_ambient_radiative_temperature trad; /* Ambient radiative temp */
- double tmax; /* Max temperature used to linearize the radiative temperature */
+
+ /* 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__ { \
@@ -394,7 +399,8 @@ struct sdis_scene_create_args {
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
- -1.0, /* Maximum temperature */ \
+ {0.0, -1.0}, /* Temperature range */ \
+ 1 /* Picard order */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
SDIS_SCENE_CREATE_ARGS_DEFAULT__;
@@ -847,18 +853,32 @@ sdis_scene_set_ambient_radiative_temperature
(struct sdis_scene* scn,
const struct sdis_ambient_radiative_temperature* trad);
-/* Get scene's maximum temperature */
+/* Get scene's minimum/maximum temperature */
SDIS_API res_T
-sdis_scene_get_maximum_temperature
+sdis_scene_get_temperature_range
(const struct sdis_scene* scn,
- double* tmax);
+ double t_range[2]);
-/* Set scene's maximum temperature. Must be correctly defined if there is any
- * radiative transfert in the scene. */
+/* Set scene's minimum/maximum temperature. Must be correctly defined if there
+ * is any radiative transfer in the scene */
SDIS_API res_T
-sdis_scene_set_maximum_temperature
+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 double tmax);
+ 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
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -26,15 +26,27 @@ struct rwalk_context {
struct green_path_handle* green_path;
struct sdis_heat_path* heat_path;
- /* That is the upper bound temperature */
+ double Tmin; /* Lower bound temperature */
+ double Tmin2; /* Tmin^2 */
+ double Tmin3; /* Tmin^3 */
+
+ double That; /* Upper bound temperature */
double That2; /* That^2 */
double That3; /* That^3 */
+
+ /* Number of heat path branchings */
+ size_t nbranchings;
};
#define RWALK_CONTEXT_NULL__ { \
NULL, /* Green path */ \
NULL, /* Heat path */ \
- 0, /* (Temperature upper bound)^2 */ \
- 0 /* (Temperature upper bound)^3 */ \
+ 0, /* Tmin */ \
+ 0, /* Tmin^2 */ \
+ 0, /* Tmin^3 */ \
+ 0, /* That */ \
+ 0, /* That^2 */ \
+ 0, /* That^3 */ \
+ SIZE_MAX, /* #branchings */ \
}
static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
@@ -112,7 +124,7 @@ static const struct XD(rwalk) XD(RWALK_NULL) = {
struct XD(temperature) {
res_T (*func)/* Next function to invoke in order to compute the temperature */
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* temp);
diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h
@@ -150,7 +150,7 @@ extern LOCAL_SYM res_T
trace_radiative_path_2d
(struct sdis_scene* scn,
const float ray_dir[3],
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* temperature);
@@ -159,7 +159,7 @@ extern LOCAL_SYM res_T
trace_radiative_path_3d
(struct sdis_scene* scn,
const float ray_dir[3],
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* temperature);
@@ -167,7 +167,7 @@ trace_radiative_path_3d
extern LOCAL_SYM res_T
radiative_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* temperature);
@@ -175,7 +175,7 @@ radiative_path_2d
extern LOCAL_SYM res_T
radiative_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* temperature);
@@ -186,7 +186,7 @@ radiative_path_3d
extern LOCAL_SYM res_T
convective_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* temperature);
@@ -194,7 +194,7 @@ convective_path_2d
extern LOCAL_SYM res_T
convective_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* temperature);
@@ -205,7 +205,7 @@ convective_path_3d
extern LOCAL_SYM res_T
conductive_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* temperature);
@@ -213,7 +213,7 @@ conductive_path_2d
extern LOCAL_SYM res_T
conductive_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* temperature);
@@ -224,7 +224,7 @@ conductive_path_3d
extern LOCAL_SYM res_T
boundary_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* temperature);
@@ -232,7 +232,7 @@ boundary_path_2d
extern LOCAL_SYM res_T
boundary_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* temperature);
diff --git a/src/sdis_heat_path_boundary.c b/src/sdis_heat_path_boundary.c
@@ -32,6 +32,10 @@
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_boundary_Xd_solid_fluid_picard1.h"
#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd_solid_fluid_picardN.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd_solid_fluid_picardN.h"
+#define SDIS_XD_DIMENSION 2
#include "sdis_heat_path_boundary_Xd_solid_solid.h"
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_boundary_Xd_solid_solid.h"
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -30,7 +30,7 @@
res_T
XD(boundary_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
@@ -89,8 +89,11 @@ XD(boundary_path)
if(mdm_front->type == mdm_back->type) {
res = XD(solid_solid_boundary_path)(scn, ctx, &frag, rwalk, rng, T);
- } else {
+ } else if(ctx->nbranchings == scn->max_branchings) {
res = XD(solid_fluid_boundary_picard1_path)(scn, ctx, &frag, rwalk, rng, T);
+ } else {
+ ASSERT(ctx->nbranchings < scn->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_heat_path_boundary_Xd_fixed_flux.h b/src/sdis_heat_path_boundary_Xd_fixed_flux.h
@@ -31,7 +31,7 @@
res_T
XD(solid_boundary_with_flux_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
const double phi,
struct XD(rwalk)* rwalk,
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -113,7 +113,7 @@ error:
res_T
XD(solid_fluid_boundary_picard1_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
@@ -263,8 +263,7 @@ XD(solid_fluid_boundary_picard1_path)
/* From there, we know the path is either a radiative path or a
* null-collision */
- /* Sample a radiative path and get the Tref at its end.
- * TODO handle the registration of the path geometry */
+ /* Sample a radiative path and get the Tref at its end. */
T_s = *T;
rwalk_s = *rwalk;
rwalk_s.mdm = fluid;
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -0,0 +1,375 @@
+/* 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
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * 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_green.h"
+#include "sdis_heat_path_boundary_c.h"
+#include "sdis_interface_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_realisation.h"
+#include "sdis_scene_c.h"
+
+#include <star/ssp.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Non generic helper functions
+ ******************************************************************************/
+#ifndef SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H
+#define SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H
+
+static INLINE res_T
+restart_heat_path
+ (struct sdis_heat_path* path,
+ const struct sdis_heat_vertex* vtx)
+{
+ size_t nverts = 0;
+ size_t nbreaks = 0;
+ res_T res = RES_OK;
+
+ if(!path) goto exit;
+ ASSERT(vtx);
+
+ nbreaks = darray_size_t_size_get(&path->breaks);
+ nverts = darray_heat_vertex_size_get(&path->vertices);
+
+ res = heat_path_add_break(path);
+ if(res != RES_OK) goto error;
+ res = heat_path_add_vertex(path, vtx);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ CHK(darray_size_t_resize(&path->breaks, nbreaks) == RES_OK);
+ CHK(darray_heat_vertex_resize(&path->vertices, nverts) == RES_OK);
+ goto exit;
+}
+
+#endif /* SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARDN_H */
+
+/*******************************************************************************
+ * Generic helper functions
+ ******************************************************************************/
+static INLINE res_T
+XD(sample_path)
+ (struct sdis_scene* scn,
+ const struct XD(rwalk)* rwalk_from,
+ struct rwalk_context* ctx,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ res_T res = RES_OK;
+ ASSERT(rwalk_from && rng && T);
+
+ /* Clean-up the output variable */
+ *T = XD(TEMPERATURE_NULL);
+ T->func = XD(boundary_path);
+
+ /* Init the random walk */
+ rwalk.vtx = rwalk_from->vtx;
+ rwalk.mdm = rwalk_from->mdm;
+ rwalk.hit = rwalk_from->hit;
+ rwalk.hit_side = rwalk_from->hit_side;
+
+ /* Start the registration of a new heat path */
+ if(ctx->heat_path) {
+ struct sdis_heat_vertex heat_vtx = SDIS_HEAT_VERTEX_NULL;
+
+ heat_vtx.P[0] = rwalk.vtx.P[0];
+ heat_vtx.P[1] = rwalk.vtx.P[1];
+ heat_vtx.P[2] = rwalk.vtx.P[2];
+ heat_vtx.time = rwalk.vtx.time;
+ heat_vtx.weight = 0;
+ heat_vtx.type = SDIS_HEAT_VERTEX_RADIATIVE;
+
+ res = restart_heat_path(ctx->heat_path, &heat_vtx);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Sample the path */
+ res = XD(compute_temperature)(scn, ctx, &rwalk, rng, T);
+ if(res != RES_OK) goto error;
+
+ /* Check the returned temperature */
+ ASSERT(T->done);
+ if(T->value < scn->tmin || scn->tmax < T->value) {
+ log_err(scn->dev, "%s: invalid temperature range `[%g, %g]K` regarding the "
+ "retrieved temperature %gK.\n", FUNC_NAME, scn->tmin, scn->tmax, T->value);
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Boundary path between a solid and a fluid
+ ******************************************************************************/
+res_T
+XD(solid_fluid_boundary_picardN_path)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ /* Input/output arguments of the function used to sample a reinjection */
+ struct XD(sample_reinjection_step_args) samp_reinject_step_args =
+ XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
+ struct XD(reinjection_step) reinject_step =
+ XD(REINJECTION_STEP_NULL);
+
+ /* Fragment on the fluid side of the boundary */
+ struct sdis_interface_fragment frag_fluid;
+
+ /* Vertex of the heat path */
+ struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL;
+ struct sdis_heat_vertex hvtx_s = SDIS_HEAT_VERTEX_NULL;
+
+ /* Data attached to the boundary */
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+
+ double h_cond; /* Conductive coefficient */
+ double h_conv; /* Convective coefficient */
+ double h_radi_hat; /* Radiative coefficient with That */
+ double h_hat; /* Sum of h_<conv|cond|rad_hat> */
+ double p_conv; /* Convective proba */
+ double p_cond; /* Conductive proba */
+
+ /* Min/Max Temperatures */
+ double Tmin, Tmin2, Tmin3;
+ double That, That2, That3;
+
+ double epsilon; /* Interface emissivity */
+ double lambda; /* Solid conductivity */
+ double delta_boundary; /* Orthogonal reinjection dst at the boundary */
+ double delta; /* Orthogonal fitted reinjection dst at the boundary */
+
+ double r;
+ enum sdis_side solid_side = SDIS_SIDE_NULL__;
+ enum sdis_side fluid_side = SDIS_SIDE_NULL__;
+ res_T res = RES_OK;
+
+ ASSERT(scn && rwalk && rng && T && ctx);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+
+ /* Fetch the Min/max temperature */
+ Tmin = ctx->Tmin;
+ Tmin2 = ctx->Tmin2;
+ Tmin3 = ctx->Tmin3;
+ That = ctx->That;
+ That2 = ctx->That2;
+ That3 = ctx->That3;
+
+ /* Retrieve the solid and the fluid split by the boundary */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ solid = interface_get_medium(interf, SDIS_FRONT);
+ fluid = interface_get_medium(interf, SDIS_BACK);
+ solid_side = SDIS_FRONT;
+ fluid_side = SDIS_BACK;
+ if(solid->type != SDIS_SOLID) {
+ SWAP(struct sdis_medium*, solid, fluid);
+ SWAP(enum sdis_side, solid_side, fluid_side);
+ ASSERT(fluid->type == SDIS_FLUID);
+ }
+
+ /* Setup a fragment for the fluid side */
+ frag_fluid = *frag;
+ frag_fluid.side = fluid_side;
+
+ /* Fetch the solid properties */
+ lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
+ delta = solid_get_delta(solid, &rwalk->vtx);
+
+ /* Fetch the boundary emissivity */
+ epsilon = interface_side_get_emissivity(interf, &frag_fluid);
+
+ /* Note that the reinjection distance is *FIXED*. It MUST ensure that the
+ * orthogonal distance from the boundary to the reinjection point is at most
+ * equal to delta. */
+ delta_boundary = sqrt(DIM) * delta;
+
+ /* Sample a reinjection step */
+ samp_reinject_step_args.rng = rng;
+ samp_reinject_step_args.solid = solid;
+ samp_reinject_step_args.rwalk = rwalk;
+ samp_reinject_step_args.distance = delta_boundary;
+ samp_reinject_step_args.side = solid_side;
+ res = XD(sample_reinjection_step_solid_fluid)
+ (scn, &samp_reinject_step_args, &reinject_step);
+ if(res != RES_OK) goto error;
+
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = reinject_step.distance / sqrt(DIM);
+
+ /* Compute the convective, conductive and the upper bound radiative coef */
+ h_conv = interface_get_convection_coef(interf, frag);
+ h_cond = lambda / (delta * scn->fp_to_meter);
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon;
+
+ /* Compute a global upper bound coefficient */
+ h_hat = h_conv + h_cond + h_radi_hat;
+
+ /* Compute the probas to switch in convection or conduction */
+ p_conv = h_conv / h_hat;
+ p_cond = h_cond / h_hat;
+
+ /* Fetch the last registered heat path vertex */
+ if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path);
+
+ /* Null collision main loop */
+ for(;;) {
+ /* Temperature and random walk state of the sampled radiative path */
+ struct XD(temperature) T_s;
+ struct XD(rwalk) rwalk_s;
+
+ double h_radi, h_radi_min, h_radi_max; /* Radiative coefficients */
+ double p_radi, p_radi_min, p_radi_max; /* Radiative probas */
+ double T0, T1, T2, T3, T4, T5; /* Computed temperatures */
+
+ r = ssp_rng_canonical(rng);
+
+ /* Switch in convective path */
+ if(r < p_conv) {
+ T->func = XD(convective_path);
+ rwalk->mdm = fluid;
+ rwalk->hit_side = fluid_side;
+ break;
+ }
+
+ /* Switch in conductive path */
+ if(r < p_conv + p_cond) {
+ struct XD(solid_reinjection_args) solid_reinject_args =
+ XD(SOLID_REINJECTION_ARGS_NULL);
+
+ /* Perform the reinjection into the solid */
+ solid_reinject_args.reinjection = &reinject_step;
+ solid_reinject_args.rwalk_ctx = ctx;
+ solid_reinject_args.rwalk = rwalk;
+ solid_reinject_args.rng = rng;
+ solid_reinject_args.T = T;
+ solid_reinject_args.fp_to_meter = scn->fp_to_meter;
+ res = XD(solid_reinjection)(solid, &solid_reinject_args);
+ if(res != RES_OK) goto error;
+ break;
+ }
+
+ /* Sample a radiative path */
+ T_s = *T;
+ rwalk_s = *rwalk;
+ rwalk_s.mdm = fluid;
+ rwalk_s.hit_side = fluid_side;
+ res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s);
+ if(res != RES_OK) goto error;
+
+ /* Fetch the last registered heat path vertex of the radiative path */
+ if(ctx->heat_path) hvtx_s = *heat_path_get_last_vertex(ctx->heat_path);
+
+ h_radi_min = 4.0 * BOLTZMANN_CONSTANT * Tmin3 * epsilon;
+ p_radi_min = h_radi_min / h_hat;
+
+ /* Switch in radiative path */
+ if(r < p_conv + p_cond + p_radi_min) { *rwalk = rwalk_s; *T = T_s; break; }
+
+ /* Define some helper macros */
+ #define SWITCH_IN_RADIATIVE { \
+ *rwalk = rwalk_s; *T = T_s; \
+ res = restart_heat_path(ctx->heat_path, &hvtx_s); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+
+ #define NULL_COLLISION { \
+ res = restart_heat_path(ctx->heat_path, &hvtx); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+
+ #define COMPUTE_TEMPERATURE(Result, RWalk) { \
+ struct XD(temperature) T_p; \
+ res = XD(sample_path)(scn, RWalk, ctx, rng, &T_p); \
+ if(res != RES_OK) goto error; \
+ Result = T_p.value; \
+ } (void)0
+
+ #define CHECK_PMIN_PMAX { \
+ p_radi_min = h_radi_min*epsilon / h_hat; \
+ p_radi_max = h_radi_max*epsilon / h_hat; \
+ if(r < p_conv + p_cond + p_radi_min) { SWITCH_IN_RADIATIVE; break; } \
+ if(r > p_conv + p_cond + p_radi_max) { NULL_COLLISION; continue; } \
+ } (void)0
+
+ /* Sample a 1st heat path at the end of the radiative path */
+ COMPUTE_TEMPERATURE(T0, &rwalk_s);
+ h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + 3*Tmin2*T0);
+ h_radi_max = BOLTZMANN_CONSTANT*(That3 + 3*That2*T0);
+ CHECK_PMIN_PMAX;
+
+ /* Sample a 2nd heat path at the end of the radiative path */
+ COMPUTE_TEMPERATURE(T1, &rwalk_s);
+ h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + Tmin2*T0 + 2*Tmin*T0*T1);
+ h_radi_max = BOLTZMANN_CONSTANT*(That3 + That2*T0 + 2*That*T0*T1);
+ CHECK_PMIN_PMAX;
+
+ /* Sample a 3rd heat path at the end of the radiative path */
+ COMPUTE_TEMPERATURE(T2, &rwalk_s);
+ h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + Tmin2*T0 + Tmin*T0*T1 + T0*T1*T2);
+ h_radi_max = BOLTZMANN_CONSTANT*(That3 + That2*T0 + That*T0*T1 + T0*T1*T2);
+ CHECK_PMIN_PMAX;
+
+ /* Sample a 1st heat path at the current position onto the interface */
+ COMPUTE_TEMPERATURE(T3, rwalk);
+ h_radi_min = BOLTZMANN_CONSTANT*(Tmin2*T3 + Tmin*T0*T3 + T0*T1*T3 + T0*T1*T2);
+ h_radi_max = BOLTZMANN_CONSTANT*(That2*T3 + That*T0*T3 + T0*T1*T3 + T0*T1*T2);
+ CHECK_PMIN_PMAX;
+
+ /* Sample a 2nd heat path at the current position onto the interface */
+ COMPUTE_TEMPERATURE(T4, rwalk);
+ h_radi_min = BOLTZMANN_CONSTANT*(Tmin*T3*T4 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2);
+ h_radi_max = BOLTZMANN_CONSTANT*(That*T3*T4 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2);
+ CHECK_PMIN_PMAX;
+
+ /* Sample a 3rd heat path at the current position onto the interface */
+ COMPUTE_TEMPERATURE(T5, rwalk);
+ h_radi = BOLTZMANN_CONSTANT*(T3*T4*T5 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2);
+ p_radi = h_radi * epsilon / h_hat;
+
+ /* Switch in radiative path */
+ if(r < p_cond + p_conv + p_radi) { SWITCH_IN_RADIATIVE; break; }
+
+ /* Null-collision, looping at the beginning */
+ NULL_COLLISION;
+
+ #undef SWITCH_IN_RADIATIVE
+ #undef NULL_COLLISION
+ #undef COMPUTE_TEMPERATURE
+ #undef CHECK_PMIN_PMAX
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
+
diff --git a/src/sdis_heat_path_boundary_Xd_solid_solid.h b/src/sdis_heat_path_boundary_Xd_solid_solid.h
@@ -31,7 +31,7 @@
res_T
XD(solid_solid_boundary_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -106,7 +106,7 @@ sample_reinjection_step_solid_solid_3d
******************************************************************************/
struct solid_reinjection_args_2d {
const struct reinjection_step_2d* reinjection; /* Reinjection to do */
- const struct rwalk_context* rwalk_ctx;
+ struct rwalk_context* rwalk_ctx;
struct rwalk_2d* rwalk; /* Current state of the random walk */
struct ssp_rng* rng; /* Random number generator */
struct temperature_2d* T;
@@ -115,7 +115,7 @@ struct solid_reinjection_args_2d {
struct solid_reinjection_args_3d {
const struct reinjection_step_3d* reinjection; /* Reinjection to do */
- const struct rwalk_context* rwalk_ctx;
+ struct rwalk_context* rwalk_ctx;
struct rwalk_3d* rwalk; /* Current state of the random walk */
struct ssp_rng* rng; /* Random number generator */
struct temperature_3d* T;
@@ -145,7 +145,7 @@ solid_reinjection_3d
extern LOCAL_SYM res_T
solid_boundary_with_flux_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
const double phi,
struct rwalk_2d* rwalk,
@@ -155,7 +155,7 @@ solid_boundary_with_flux_path_2d
extern LOCAL_SYM res_T
solid_boundary_with_flux_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
const double phi,
struct rwalk_3d* rwalk,
@@ -163,36 +163,36 @@ solid_boundary_with_flux_path_3d
struct temperature_3d* T);
extern LOCAL_SYM res_T
-solid_fluid_boundary_path_2d
+solid_fluid_boundary_picard1_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* T);
extern LOCAL_SYM res_T
-solid_fluid_boundary_path_3d
+solid_fluid_boundary_picard1_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
struct temperature_3d* T);
extern LOCAL_SYM res_T
-solid_fluid_boundary_picard1_path_2d
+solid_fluid_boundary_picardN_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
struct temperature_2d* T);
extern LOCAL_SYM res_T
-solid_fluid_boundary_picard1_path_3d
+solid_fluid_boundary_picardN_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
@@ -201,7 +201,7 @@ solid_fluid_boundary_picard1_path_3d
extern LOCAL_SYM res_T
solid_solid_boundary_path_2d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_2d* rwalk,
struct ssp_rng* rng,
@@ -210,7 +210,7 @@ solid_solid_boundary_path_2d
extern LOCAL_SYM res_T
solid_solid_boundary_path_3d
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
struct rwalk_3d* rwalk,
struct ssp_rng* rng,
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -199,7 +199,7 @@ error:
res_T
XD(conductive_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -29,7 +29,7 @@
static res_T
XD(register_heat_vertex_in_fluid)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
const double weight)
{
@@ -70,7 +70,7 @@ XD(register_heat_vertex_in_fluid)
res_T
XD(convective_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -33,7 +33,7 @@ res_T
XD(trace_radiative_path)
(struct sdis_scene* scn,
const float ray_dir[3],
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
@@ -188,7 +188,7 @@ error:
res_T
XD(radiative_path)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -47,8 +47,12 @@ ray_realisation_3d
rwalk.mdm = medium;
ctx.heat_path = heat_path;
- ctx.That2 = scn->tmax * scn->tmax;
- ctx.That3 = scn->tmax * ctx.That2;
+ 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;
f3_set_d3(dir, direction);
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -35,6 +35,25 @@ enum flux_flag {
};
/*******************************************************************************
+ * Helper function used to compute a temperature
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+compute_temperature_2d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_2d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_2d* T);
+
+extern LOCAL_SYM res_T
+compute_temperature_3d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T);
+
+/*******************************************************************************
* Realisation at a given position and time IN a medium
******************************************************************************/
extern LOCAL_SYM res_T
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -27,12 +27,12 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
- * Helper functions
+ * Local functions
******************************************************************************/
-static res_T
+res_T
XD(compute_temperature)
(struct sdis_scene* scn,
- const struct rwalk_context* ctx,
+ struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
@@ -51,11 +51,14 @@ XD(compute_temperature)
res_T res = RES_OK;
ASSERT(scn && ctx && rwalk && rng && T);
+ ctx->nbranchings += 1;
+ CHK(ctx->nbranchings <= scn->max_branchings);
+
if(ctx->heat_path && T->func == XD(boundary_path)) {
heat_vtx = heat_path_get_last_vertex(ctx->heat_path);
}
- do {
+ while(!T->done) {
/* Save the current random walk state */
const struct XD(rwalk) rwalk_bkp = *rwalk;
const struct XD(temperature) T_bkp = *T;
@@ -94,22 +97,19 @@ XD(compute_temperature)
}
heat_vtx = NULL; /* Notify that the first vertex is finalized */
}
+ }
- } while(!T->done);
exit:
#ifndef NDEBUG
sa_release(stack);
#endif
+ ctx->nbranchings -= 1;
return res == RES_BAD_OP_IRRECOVERABLE ? RES_BAD_OP : res;
error:
goto exit;
}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
res_T
XD(probe_realisation)
(const size_t irealisation, /* For debug */
@@ -181,8 +181,12 @@ XD(probe_realisation)
ctx.green_path = green_path;
ctx.heat_path = heat_path;
- ctx.That2 = scn->tmax * scn->tmax;
- ctx.That3 = scn->tmax * ctx.That2;
+ 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;
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
@@ -256,8 +260,12 @@ XD(boundary_realisation)
ctx.green_path = green_path;
ctx.heat_path = heat_path;
- ctx.That2 = scn->tmax * scn->tmax;
- ctx.That3 = scn->tmax * ctx.That2;
+ 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;
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
@@ -296,6 +304,9 @@ 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;
@@ -333,6 +344,9 @@ XD(boundary_flux_realisation)
rwalk.mdm = (Mdm); \
rwalk.hit.prim = prim; \
SET_PARAM(rwalk.hit, st); \
+ ctx.Tmin = Tmin; \
+ ctx.Tmin3 = Tmin3; \
+ ctx.That = That; \
ctx.That2 = That2; \
ctx.That3 = That3; \
dX(set)(rwalk.vtx.P, P); \
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -221,18 +221,44 @@ sdis_scene_set_ambient_radiative_temperature
}
res_T
-sdis_scene_get_maximum_temperature(const struct sdis_scene* scn, double* tmax)
+sdis_scene_get_temperature_range
+ (const struct sdis_scene* scn,
+ double t_range[2])
+{
+ if(!scn || !t_range) return RES_BAD_ARG;
+ t_range[0] = scn->tmin;
+ t_range[1] = scn->tmax;
+ return RES_OK;
+}
+
+res_T
+sdis_scene_set_temperature_range
+ (struct sdis_scene* scn,
+ const double t_range[2])
{
- if(!scn || !tmax) return RES_BAD_ARG;
- *tmax = scn->tmax;
+ if(!scn || !t_range) return RES_BAD_ARG;
+ scn->tmin = t_range[0];
+ scn->tmax = t_range[1];
return RES_OK;
}
res_T
-sdis_scene_set_maximum_temperature(struct sdis_scene* scn, const double tmax)
+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 || tmax < 0) return RES_BAD_ARG;
- scn->tmax = tmax;
+ if(!scn || picard_order < 1) return RES_BAD_ARG;
+ scn->max_branchings = picard_order-1;
return RES_OK;
}
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -120,7 +120,8 @@ 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->fp_to_meter > 0
+ && args->picard_order > 0;
}
#endif /* SDIS_SCENE_XD_H */
@@ -913,8 +914,10 @@ XD(scene_create)
scn->dev = dev;
scn->fp_to_meter = args->fp_to_meter;
scn->trad = args->trad;
- scn->tmax = args->tmax;
+ 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
@@ -210,8 +210,15 @@ struct sdis_scene {
double fp_to_meter;
struct sdis_ambient_radiative_temperature trad;
+ 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;
};
@@ -296,5 +303,13 @@ 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_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -119,6 +119,14 @@ XD(solve_boundary)
res = RES_BAD_ARG;
goto error;
}
+ if(out_green && scene_get_picard_order(scn) != 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));
+ res = RES_BAD_ARG;
+ goto error;
+ }
#if SDIS_XD_DIMENSION == 2
if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; }
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -240,6 +240,15 @@ XD(solve_medium)
}
}
+ if(out_green && scene_get_picard_order(scn) != 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));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
#if SDIS_XD_DIMENSION == 2
if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; }
#else
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -70,6 +70,14 @@ XD(solve_probe)
res = RES_BAD_ARG;
goto error;
}
+ if(out_green && scene_get_picard_order(scn) != 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));
+ res = RES_BAD_ARG;
+ goto error;
+ }
#if SDIS_XD_DIMENSION == 2
if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; }
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -70,6 +70,14 @@ XD(solve_probe_boundary)
res = RES_BAD_ARG;
goto error;
}
+ if(out_green && scene_get_picard_order(scn) != 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));
+ res = RES_BAD_ARG;
+ goto error;
+ }
#if SDIS_XD_DIMENSION == 2
if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; }
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -278,6 +278,53 @@ create_interface
}
/*******************************************************************************
+ * Test that the evaluation of the green function failed with a picard order
+ * greater than 1, i.e. when one want to handle the non-linearties of the
+ * system.
+ ******************************************************************************/
+static void
+test_invalidity_picardN_green
+ (struct sdis_scene* scn,
+ struct sdis_medium* solid)
+{
+ struct sdis_solve_probe_args probe = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_solve_boundary_args bound = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_solve_medium_args mdm = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT;
+ struct sdis_solve_probe_boundary_args probe_bound =
+ 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));
+
+ probe.position[0] = 0;
+ probe.position[1] = 0;
+ probe.position[2] = 0;
+ 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;
+ BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green));
+
+ bound.primitives = &probe_bound.iprim;
+ bound.sides = &probe_bound.side;
+ bound.nprimitives = 1;
+ BA(sdis_solve_boundary_green_function(scn, &bound, &green));
+
+ mdm.medium = solid;
+ BA(sdis_solve_medium_green_function(scn, &mdm, &green));
+
+ OK(sdis_scene_set_picard_order(scn, picard_order));
+}
+
+/*******************************************************************************
* Test
******************************************************************************/
int
@@ -306,6 +353,7 @@ main(int argc, char** argv)
const double T0 = 300; /* Fixed temperature on the left side of the system */
const double T1 = 310; /* Fixed temperature on the right side of the system */
const double thickness = 2.0; /* Thickness of the solid along X */
+ double t_range[2];
double Ts0, Ts1, hr, tmp;
struct interfac* p_intface;
(void)argc, (void)argv;
@@ -414,7 +462,8 @@ main(int argc, char** argv)
scn_args.get_position = get_position;
scn_args.nprimitives = ntriangles;
scn_args.nvertices = nvertices;
- scn_args.tmax = MMAX(T0, T1);
+ scn_args.t_range[0] = MMIN(T0, T1);
+ scn_args.t_range[1] = MMAX(T0, T1);
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, &scn));
@@ -492,7 +541,9 @@ main(int argc, char** argv)
/* Check same green used at a different temperature */
p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10;
- OK(sdis_scene_set_maximum_temperature(scn, MMAX(T0, T1b)));
+ t_range[0] = MMIN(T0, T1b);
+ t_range[1] = MMAX(T0, T1b);
+ OK(sdis_scene_set_temperature_range(scn, t_range));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
@@ -539,6 +590,8 @@ main(int argc, char** argv)
printf("\n\n");
}
+ test_invalidity_picardN_green(scn, solid);
+
/* Release memory */
OK(sdis_scene_ref_put(scn));
OK(sdis_interface_ref_put(interfaces[0]));
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -285,6 +285,51 @@ create_interface
}
/*******************************************************************************
+ * Test that the evaluation of the green function failed with a picard order
+ * greater than 1, i.e. when one want to handle the non-linearties of the
+ * system.
+ ******************************************************************************/
+static void
+test_invalidity_picardN_green
+ (struct sdis_scene* scn,
+ struct sdis_medium* solid)
+{
+ struct sdis_solve_probe_args probe = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_solve_boundary_args bound = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_solve_medium_args mdm = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT;
+ struct sdis_solve_probe_boundary_args probe_bound =
+ 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));
+
+ probe.position[0] = 0;
+ probe.position[1] = 0;
+ 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;
+ BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green));
+
+ bound.primitives = &probe_bound.iprim;
+ bound.sides = &probe_bound.side;
+ bound.nprimitives = 1;
+ BA(sdis_solve_boundary_green_function(scn, &bound, &green));
+
+ mdm.medium = solid;
+ BA(sdis_solve_medium_green_function(scn, &mdm, &green));
+
+ OK(sdis_scene_set_picard_order(scn, picard_order));
+}
+
+/*******************************************************************************
* Test
******************************************************************************/
int
@@ -408,7 +453,8 @@ main(int argc, char** argv)
scn_args.nprimitives = nsegments;
scn_args.nvertices = nvertices;
scn_args.context = &geom;
- scn_args.tmax = MMAX(T0, T1);
+ scn_args.t_range[0] = MMIN(T0, T1);
+ scn_args.t_range[1] = MMAX(T0, T1);
OK(sdis_scene_2d_create(dev, &scn_args, &scn));
hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
@@ -473,6 +519,8 @@ main(int argc, char** argv)
printf("\n\n");
}
+ test_invalidity_picardN_green(scn, solid);
+
/* Release memory */
OK(sdis_scene_ref_put(scn));
OK(sdis_interface_ref_put(interfaces[0]));
diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c
@@ -0,0 +1,774 @@
+/* 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
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * 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.h"
+#include "test_sdis_utils.h"
+
+#include <string.h>
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000
+
+/* This test consists in solving the stationary temperature profile in a solid
+ * slab surrounded by two different radiative temperatures (left / right). The
+ * conductivity of the solid material is known, as well as its thickness and
+ * the source term (volumic power density).
+ *
+ * The purpose is to test the Picard-1 radiative transfer algorithm, that can
+ * be compared with analytic results. This algorithm can use a possibly
+ * non-uniform reference temperature field. When the reference temperature
+ * field is uniform, results should be identical to the classical Monte-Carlo
+ * algorithm (using a linearized radiative transfer scheme).
+ *
+ * Y
+ * | (0.1,1)
+ * o--- X +------+----------+------+ (1.1,1)
+ * | |##########| |
+ * | |##########| |
+ * 280K | E=1 |##########| E=1 | 350K
+ * | |##########| |
+ * | |##########| |
+ * (-1,-1) +------+----------+------+
+ * (0,-1)
+ *
+ *
+ *
+ * Y (0.1, 1, 1)
+ * | +------+----------+------+ (1.1,1,1)
+ * o--- X /' /##########/' /|
+ * / +------+----------+------+ |
+ * Z | ' |##########|*' | | 350K
+ * | ' |##########|*' | |
+ * 280K | ' E=1|##########|*'E=1 | |
+ * | +....|##########|*+....|.+
+ * |/ |##########|/ |/
+ * (-1,-1,-1) +------+----------+------+
+ * (0,-1,-1)
+ *
+ * lambda = 1.15 W/(m.K)
+ * rho = 1000 kg.m^-3
+ * cp = 800 J/(kg.K)
+ * emissivity = 1
+ *
+ * Basic Tref = 300 K
+ * probe = 0.05 0 0 m
+ * (power = 1000 W.m^-3) */
+
+enum interface_type {
+ ADIABATIC,
+ SOLID_FLUID_mX,
+ SOLID_FLUID_pX,
+ BOUNDARY_mX,
+ BOUNDARY_pX,
+ INTERFACES_COUNT__
+};
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct geometry {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface** interfaces;
+};
+
+static const double vertices_2d[8/*#vertices*/*2/*#coords par vertex*/] = {
+ 0.1, -1.0,
+ 0.0, -1.0,
+ 0.0, 1.0,
+ 0.1, 1.0,
+ 1.1, -1.0,
+ -1.0, -1.0,
+ -1.0, 1.0,
+ 1.1, 1.0
+};
+static const size_t nvertices_2d = sizeof(vertices_2d) / (sizeof(double)*2);
+
+static const size_t indices_2d[10/*#segments*/*2/*#indices per segment*/] = {
+ 0, 1, /* Solid -Y */
+ 1, 2, /* Solid -X */
+ 2, 3, /* Solid +Y */
+ 3, 0, /* Solid +X */
+
+ 1, 5, /* Left fluid -Y */
+ 5, 6, /* Left fluid -X */
+ 6, 2, /* Left fluid +Y */
+
+ 4, 0, /* Right fluid -Y */
+ 3, 7, /* Right fluid +Y */
+ 7, 4 /* Right fluid +X */
+};
+static const size_t nprimitives_2d = sizeof(indices_2d) / (sizeof(size_t)*2);
+
+static const double vertices_3d[16/*#vertices*/*3/*#coords per vertex*/] = {
+ 0.0,-1.0,-1.0,
+ 0.1,-1.0,-1.0,
+ 0.0, 1.0,-1.0,
+ 0.1, 1.0,-1.0,
+ 0.0,-1.0, 1.0,
+ 0.1,-1.0, 1.0,
+ 0.0, 1.0, 1.0,
+ 0.1, 1.0, 1.0,
+ -1.0,-1.0,-1.0,
+ 1.1,-1.0,-1.0,
+ -1.0, 1.0,-1.0,
+ 1.1, 1.0,-1.0,
+ -1.0,-1.0, 1.0,
+ 1.1,-1.0, 1.0,
+ -1.0, 1.0, 1.0,
+ 1.1, 1.0, 1.0,
+};
+static const size_t nvertices_3d = sizeof(vertices_3d) / (sizeof(double)*3);
+
+static const size_t indices_3d[32/*#triangles*/*3/*#indices per triangle*/] = {
+ 0, 2, 1, 1, 2, 3, /* Solid -Z */
+ 0, 4, 2, 2, 4, 6, /* Solid -X */
+ 4, 5, 6, 6, 5, 7, /* Solid +Z */
+ 3, 7, 1, 1, 7, 5, /* Solid +X */
+ 2, 6, 3, 3, 6, 7, /* Solid +Y */
+ 0, 1, 4, 4, 1, 5, /* Solid -Y */
+
+ 8, 10, 0, 0, 10, 2, /* Left fluid -Z */
+ 8, 12, 10, 10, 12, 14, /* Left fluid -X */
+ 12, 4, 14, 14, 4, 6, /* Left fluid +Z */
+ 10, 14, 2, 2, 14, 6, /* Left fluid +Y */
+ 8, 0, 12, 12, 0, 4, /* Left fluid -Y */
+
+ 1, 3, 9, 9, 3, 11, /* Right fluid -Z */
+ 5, 13, 7, 7, 13, 15, /* Right fluid +Z */
+ 11, 15, 9, 9, 15, 13, /* Right fluid +X */
+ 3, 7, 11, 11, 7, 15, /* Right fluid +Y */
+ 1, 9, 5, 5, 9, 13 /* Right fluid -Y */
+};
+static const size_t nprimitives_3d = sizeof(indices_3d) / (sizeof(size_t)*3);
+
+static void
+get_indices_2d(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ ids[0] = geom->indices[iseg*2+0];
+ ids[1] = geom->indices[iseg*2+1];
+}
+
+static void
+get_indices_3d(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ ids[0] = geom->indices[itri*3+0];
+ ids[1] = geom->indices[itri*3+1];
+ ids[2] = geom->indices[itri*3+2];
+}
+
+static void
+get_position_2d(const size_t ivert, double pos[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ pos[0] = geom->positions[ivert*2+0];
+ pos[1] = geom->positions[ivert*2+1];
+}
+
+static void
+get_position_3d(const size_t ivert, double pos[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ pos[0] = geom->positions[ivert*3+0];
+ pos[1] = geom->positions[ivert*3+1];
+ pos[2] = geom->positions[ivert*3+2];
+}
+
+static void
+get_interface(const size_t iprim, struct sdis_interface** bound, void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ *bound = geom->interfaces[iprim];
+}
+
+/*******************************************************************************
+ * media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double volumic_power;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct solid* solid = sdis_data_cget(data);
+ CHK(vtx && solid);
+ return solid->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct solid* solid = sdis_data_cget(data);
+ CHK(vtx && solid);
+ return solid->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct solid* solid = sdis_data_cget(data);
+ CHK(vtx && solid);
+ return solid->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return 0.005;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return UNKNOWN_TEMPERATURE;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct solid* solid = sdis_data_cget(data);
+ CHK(vtx && solid);
+ return solid->volumic_power;
+}
+
+static void
+create_solid
+ (struct sdis_device* dev,
+ const struct solid* solid_props,
+ struct sdis_medium** solid)
+{
+ struct sdis_data* data = NULL;
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ CHK(dev && solid_props && solid);
+
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ memcpy(sdis_data_get(data), solid_props, sizeof(struct solid));
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta_solid = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.volumic_power = solid_get_volumic_power;
+ OK(sdis_solid_create(dev, &shader, data, solid));
+ OK(sdis_data_ref_put(data));
+}
+
+static void
+create_fluid(struct sdis_device* dev, struct sdis_medium** fluid)
+{
+ struct sdis_fluid_shader shader = DUMMY_FLUID_SHADER;
+ OK(sdis_fluid_create(dev, &shader, NULL, fluid));
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double h;
+ double emissivity;
+ double specular_fraction;
+ double Tref;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && interf);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && interf);
+ return interf->h;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && interf);
+ return interf->emissivity;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && interf);
+ return interf->specular_fraction;
+}
+
+static double
+interface_get_Tref
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && interf);
+ return interf->Tref;
+}
+
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct interf* interf,
+ struct sdis_interface** out_interf)
+{
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_data* data = NULL;
+
+ CHK(interf != NULL);
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.convection_coef = interface_get_convection_coef;
+ }
+ if(sdis_medium_get_type(front) == SDIS_FLUID) {
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.specular_fraction = interface_get_specular_fraction;
+ shader.front.reference_temperature = interface_get_Tref;
+ }
+ if(sdis_medium_get_type(back) == SDIS_FLUID) {
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.specular_fraction = interface_get_specular_fraction;
+ shader.back.reference_temperature = interface_get_Tref;
+ }
+ shader.convection_coef_upper_bound = MMAX(0, interf->h);
+
+ OK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
+ memcpy(sdis_data_get(data), interf, sizeof(*interf));
+
+ OK(sdis_interface_create(dev, front, back, &shader, data, out_interf));
+ OK(sdis_data_ref_put(data));
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+struct reference_result {
+ double T; /* Temperature at the center of the solid [K] */
+ double T1; /* Temperature on the left boundary of the solid [K] */
+ double T2; /* Temperature on the right boundary of the solid [K] */
+};
+static const struct reference_result REFERENCE_RESULT_NULL = {0,0,0};
+
+static void
+test_picard1(struct sdis_scene* scn, 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;
+ struct sdis_mc mc = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ enum sdis_scene_dimension dim;
+ size_t prims[2];
+ enum sdis_side sides[2];
+ CHK(scn);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+ switch(dim) {
+ case SDIS_SCENE_2D: printf(">>> 2D\n"); break;
+ case SDIS_SCENE_3D: printf(">>> 3D\n"); break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ probe_args.nrealisations = N;
+ probe_args.position[0] = 0.05;
+ probe_args.position[1] = 0;
+ probe_args.position[2] = 0;
+ 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",
+ SPLIT3(probe_args.position), ref->T, mc.E, mc.SE);
+ CHK(eq_eps(ref->T, mc.E, mc.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ prims[0] = 1; sides[0] = SDIS_BACK;
+ bound_args.nprimitives = 1;
+ break;
+ case SDIS_SCENE_3D:
+ prims[0] = 2; sides[0] = SDIS_BACK;
+ prims[1] = 3; sides[1] = SDIS_BACK;
+ bound_args.nprimitives = 2;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ bound_args.nrealisations = N;
+ bound_args.primitives = prims;
+ bound_args.sides = sides;
+ 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);
+ CHK(eq_eps(ref->T1, mc.E, mc.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ prims[0] = 3; sides[0] = SDIS_BACK;
+ bound_args.nprimitives = 1;
+ break;
+ case SDIS_SCENE_3D:
+ prims[0] = 6; sides[0] = SDIS_BACK;
+ prims[1] = 7; sides[1] = SDIS_BACK;
+ bound_args.nprimitives = 2;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ bound_args.nrealisations = N;
+ bound_args.primitives = prims;
+ bound_args.sides = sides;
+ 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);
+ CHK(eq_eps(ref->T2, mc.E, mc.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+static void
+create_scene_3d
+ (struct sdis_device* dev,
+ struct sdis_interface* interfaces[INTERFACES_COUNT__],
+ struct sdis_scene** scn)
+{
+ struct geometry geom;
+ struct sdis_interface* prim_interfaces[32];
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ CHK(dev && interfaces && scn);
+
+ /* Setup the per primitive interface of the solid medium */
+ prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC];
+ prim_interfaces[2] = prim_interfaces[3] = interfaces[SOLID_FLUID_mX];
+ prim_interfaces[4] = prim_interfaces[5] = interfaces[ADIABATIC];
+ prim_interfaces[6] = prim_interfaces[7] = interfaces[SOLID_FLUID_pX];
+ prim_interfaces[8] = prim_interfaces[9] = interfaces[ADIABATIC];
+ prim_interfaces[10] = prim_interfaces[11] = interfaces[ADIABATIC];
+
+ /* Setup the per primitive interface for the left fluid */
+ prim_interfaces[12] = prim_interfaces[13] = interfaces[BOUNDARY_mX];
+ prim_interfaces[14] = prim_interfaces[15] = interfaces[BOUNDARY_mX];
+ prim_interfaces[16] = prim_interfaces[17] = interfaces[BOUNDARY_mX];
+ prim_interfaces[18] = prim_interfaces[19] = interfaces[BOUNDARY_mX];
+ prim_interfaces[20] = prim_interfaces[21] = interfaces[BOUNDARY_mX];
+
+ /* Setup the per primitive interface for the right fluid */
+ prim_interfaces[22] = prim_interfaces[23] = interfaces[BOUNDARY_pX];
+ prim_interfaces[24] = prim_interfaces[25] = interfaces[BOUNDARY_pX];
+ prim_interfaces[26] = prim_interfaces[27] = interfaces[BOUNDARY_pX];
+ prim_interfaces[28] = prim_interfaces[29] = interfaces[BOUNDARY_pX];
+ prim_interfaces[30] = prim_interfaces[31] = interfaces[BOUNDARY_pX];
+
+ /* Create the scene */
+ geom.positions = vertices_3d;
+ geom.indices = indices_3d;
+ geom.interfaces = prim_interfaces;
+ scn_args.get_indices = get_indices_3d;
+ scn_args.get_interface = get_interface;
+ scn_args.get_position = get_position_3d;
+ scn_args.nprimitives = nprimitives_3d;
+ scn_args.nvertices = nvertices_3d;
+ scn_args.t_range[0] = 280;
+ scn_args.t_range[1] = 350;
+ scn_args.context = &geom;
+ OK(sdis_scene_create(dev, &scn_args, scn));
+}
+
+static void
+create_scene_2d
+ (struct sdis_device* dev,
+ struct sdis_interface* interfaces[INTERFACES_COUNT__],
+ struct sdis_scene** scn)
+{
+ struct geometry geom;
+ struct sdis_interface* prim_interfaces[10/*#segment*/];
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ CHK(dev && interfaces && scn);
+
+ /* Setup the per primitive interface of the solid medium */
+ prim_interfaces[0] = interfaces[ADIABATIC];
+ prim_interfaces[1] = interfaces[SOLID_FLUID_mX];
+ prim_interfaces[2] = interfaces[ADIABATIC];
+ prim_interfaces[3] = interfaces[SOLID_FLUID_pX];
+
+ /* Setup the per primitive interface of the fluid on the left of the medium */
+ prim_interfaces[4] = interfaces[BOUNDARY_mX];
+ prim_interfaces[5] = interfaces[BOUNDARY_mX];
+ prim_interfaces[6] = interfaces[BOUNDARY_mX];
+
+ /* Setup the per primitive interface of the fluid on the right of the medium */
+ prim_interfaces[7] = interfaces[BOUNDARY_pX];
+ prim_interfaces[8] = interfaces[BOUNDARY_pX];
+ prim_interfaces[9] = interfaces[BOUNDARY_pX];
+
+ /* Create the scene */
+ geom.positions = vertices_2d;
+ geom.indices = indices_2d;
+ geom.interfaces = prim_interfaces;
+ scn_args.get_indices = get_indices_2d;
+ scn_args.get_interface = get_interface;
+ scn_args.get_position = get_position_2d;
+ scn_args.nprimitives = nprimitives_2d;
+ scn_args.nvertices = nvertices_2d;
+ scn_args.t_range[0] = 280;
+ scn_args.t_range[1] = 350;
+ scn_args.context = &geom;
+ OK(sdis_scene_2d_create(dev, &scn_args, scn));
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+
+ struct sdis_device* dev = NULL;
+ struct sdis_scene* scn_2d = NULL;
+ struct sdis_scene* scn_3d = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* dummy = NULL;
+ struct sdis_interface* interfaces[INTERFACES_COUNT__];
+
+ struct solid solid_props;
+ struct solid* psolid_props;
+ struct reference_result ref = REFERENCE_RESULT_NULL;
+ struct interf interf_props;
+ struct interf* pinterf_props[INTERFACES_COUNT__];
+
+ double t_range[2];
+
+ size_t i;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Solid medium */
+ solid_props.lambda = 1.15;
+ solid_props.rho = 1000;
+ solid_props.cp = 800;
+ create_solid(dev, &solid_props, &solid);
+
+ /* Dummy solid medium */
+ solid_props.lambda = 0;
+ solid_props.rho = 1000;
+ solid_props.cp = 800;
+ solid_props.volumic_power = SDIS_VOLUMIC_POWER_NONE;
+ create_solid(dev, &solid_props, &dummy);
+
+ /* Fluid medium */
+ create_fluid(dev, &fluid);
+
+ /* Create the adiabatic interface for the solid */
+ interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.h = -1;
+ interf_props.emissivity = -1;
+ interf_props.specular_fraction = -1;
+ interf_props.Tref = UNKNOWN_TEMPERATURE;
+ create_interface(dev, solid, dummy, &interf_props, interfaces+ADIABATIC);
+
+ /* Create the interface between the solid and the fluid */
+ interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.h = 0;
+ interf_props.emissivity = 1;
+ interf_props.specular_fraction = 0;
+ interf_props.Tref = 280;
+ create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_mX);
+ interf_props.Tref = 350;
+ create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_pX);
+
+ /* Create the interface for the fluid on the left */
+ interf_props.temperature = 280;
+ interf_props.h = -1;
+ interf_props.emissivity = 1;
+ interf_props.specular_fraction = -1;
+ interf_props.Tref = 280;
+ create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_mX);
+
+ /* Create the interface for the fluid on the right */
+ interf_props.temperature = 350;
+ interf_props.h = -1;
+ interf_props.emissivity = 1;
+ interf_props.specular_fraction = -1;
+ interf_props.Tref = 350;
+ create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_pX);
+
+ /* Fetch pointers toward the solid and the interfaces */
+ psolid_props = sdis_data_get(sdis_medium_get_data(solid));
+ FOR_EACH(i, 0, INTERFACES_COUNT__) {
+ pinterf_props[i] = sdis_data_get(sdis_interface_get_data(interfaces[i]));
+ }
+
+ create_scene_2d(dev, interfaces, &scn_2d);
+ create_scene_3d(dev, interfaces, &scn_3d);
+
+ /* Test picard1 with a constant Tref <=> regular linearisation */
+ printf("Test Picard1 with a constant Tref of 300 K\n");
+ ref.T = 314.99999999999989;
+ ref.T1 = 307.64122364709766;
+ ref.T2 = 322.35877635290217;
+ pinterf_props[SOLID_FLUID_mX]->Tref = 300;
+ 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);
+ printf("\n");
+
+ /* Test picard1 using T4 as a reference */
+ printf("Test Picard1 using T4 as a reference\n");
+ ref.T = 320.37126474482994;
+ ref.T1 = 312.12650299072266;
+ ref.T2 = 328.61602649893723;
+ pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
+ 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);
+ printf("\n");
+
+ /* Test picard2 */
+ printf("Test Picard2 with a constant Tref of 300K\n");
+ ref.T = 320.37126474482994;
+ ref.T1 = 312.12650299072266;
+ ref.T2 = 328.61602649893723;
+ pinterf_props[SOLID_FLUID_mX]->Tref = 300;
+ 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));
+ 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");
+ ref.T = 416.4023;
+ ref.T1 = 372.7557;
+ ref.T2 = 460.0489;
+ pinterf_props[BOUNDARY_mX]->temperature = t_range[0];
+ pinterf_props[BOUNDARY_pX]->temperature = t_range[1];
+ pinterf_props[SOLID_FLUID_mX]->Tref = 350;
+ 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);
+ 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];
+
+ /* Add volumic power */
+ psolid_props->volumic_power = 1000;
+
+ /* Test picard1 with a volumic power and constant Tref */
+ printf("Test Picard1 with a volumic power of 1000 W/m^3 and a constant Tref "
+ "of 300 K\n");
+ ref.T = 324.25266420769509;
+ ref.T1 = 315.80693133305368;
+ ref.T2 = 330.52448403885825;
+ pinterf_props[SOLID_FLUID_mX]->Tref = 300;
+ 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);
+ printf("\n");
+
+ /* Test picard1 with a volumic power and T4 a the reference */
+ printf("Test Picard1 with a volumic power of 1000 W/m^3 and T4 as a reference\n");
+ ref.T = 327.95981050850446;
+ ref.T1 = 318.75148773193359;
+ ref.T2 = 334.99422024159708;
+ pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
+ 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);
+ printf("\n");
+
+ /* Release memory */
+ OK(sdis_scene_ref_put(scn_2d));
+ OK(sdis_scene_ref_put(scn_3d));
+ OK(sdis_interface_ref_put(interfaces[ADIABATIC]));
+ OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_mX]));
+ OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_pX]));
+ OK(sdis_interface_ref_put(interfaces[BOUNDARY_mX]));
+ OK(sdis_interface_ref_put(interfaces[BOUNDARY_pX]));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_device_ref_put(dev));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_picard1.c b/src/test_sdis_picard1.c
@@ -1,722 +0,0 @@
-/* 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
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * 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.h"
-#include "test_sdis_utils.h"
-
-#include <string.h>
-
-#define UNKNOWN_TEMPERATURE -1
-#define N 10000
-
-/* This test consists in solving the stationary temperature profile in a solid
- * slab surrounded by two different radiative temperatures (left / right). The
- * conductivity of the solid material is known, as well as its thickness and
- * the source term (volumic power density).
- *
- * The purpose is to test the Picard-1 radiative transfer algorithm, that can
- * be compared with analytic results. This algorithm can use a possibly
- * non-uniform reference temperature field. When the reference temperature
- * field is uniform, results should be identical to the classical Monte-Carlo
- * algorithm (using a linearized radiative transfer scheme).
- *
- * Y
- * | (0.1,1)
- * o--- X +------+----------+------+ (1.1,1)
- * | |##########| |
- * | |##########| |
- * 300K | E=1 |##########| E=1 | 350K
- * | |##########| |
- * | |##########| |
- * (-1,-1) +------+----------+------+
- * (0,-1)
- *
- *
- *
- * Y (0.1, 1, 1)
- * | +------+----------+------+ (1.1,1,1)
- * o--- X /' /##########/' /|
- * / +------+----------+------+ |
- * Z | ' |##########|*' | | 350K
- * | ' |##########|*' | |
- * 280K | ' E=1|##########|*'E=1 | |
- * | +....|##########|*+....|.+
- * |/ |##########|/ |/
- * (-1,-1,-1) +------+----------+------+
- * (0,-1,-1)
- *
- * lambda = 1.15 W/(m.K)
- * rho = 1000 kg.m^-3
- * cp = 800 J/(kg.K)
- * emissivity = 1
- *
- * Basic Tref = 300 K
- * probe = 0.05 0 0 m
- * (power = 1000 W.m^-3) */
-
-enum interface_type {
- ADIABATIC,
- SOLID_FLUID_mX,
- SOLID_FLUID_pX,
- BOUNDARY_mX,
- BOUNDARY_pX,
- INTERFACES_COUNT__
-};
-
-/*******************************************************************************
- * Geometry
- ******************************************************************************/
-struct geometry {
- const double* positions;
- const size_t* indices;
- struct sdis_interface** interfaces;
-};
-
-static const double vertices_2d[8/*#vertices*/*2/*#coords par vertex*/] = {
- 0.1, -1.0,
- 0.0, -1.0,
- 0.0, 1.0,
- 0.1, 1.0,
- 1.1, -1.0,
- -1.0, -1.0,
- -1.0, 1.0,
- 1.1, 1.0
-};
-static const size_t nvertices_2d = sizeof(vertices_2d) / (sizeof(double)*2);
-
-static const size_t indices_2d[10/*#segments*/*2/*#indices per segment*/] = {
- 0, 1, /* Solid -Y */
- 1, 2, /* Solid -X */
- 2, 3, /* Solid +Y */
- 3, 0, /* Solid +X */
-
- 1, 5, /* Left fluid -Y */
- 5, 6, /* Left fluid -X */
- 6, 2, /* Left fluid +Y */
-
- 4, 0, /* Right fluid -Y */
- 3, 7, /* Right fluid +Y */
- 7, 4 /* Right fluid +X */
-};
-static const size_t nprimitives_2d = sizeof(indices_2d) / (sizeof(size_t)*2);
-
-static const double vertices_3d[16/*#vertices*/*3/*#coords per vertex*/] = {
- 0.0,-1.0,-1.0,
- 0.1,-1.0,-1.0,
- 0.0, 1.0,-1.0,
- 0.1, 1.0,-1.0,
- 0.0,-1.0, 1.0,
- 0.1,-1.0, 1.0,
- 0.0, 1.0, 1.0,
- 0.1, 1.0, 1.0,
- -1.0,-1.0,-1.0,
- 1.1,-1.0,-1.0,
- -1.0, 1.0,-1.0,
- 1.1, 1.0,-1.0,
- -1.0,-1.0, 1.0,
- 1.1,-1.0, 1.0,
- -1.0, 1.0, 1.0,
- 1.1, 1.0, 1.0,
-};
-static const size_t nvertices_3d = sizeof(vertices_3d) / (sizeof(double)*3);
-
-static const size_t indices_3d[32/*#triangles*/*3/*#indices per triangle*/] = {
- 0, 2, 1, 1, 2, 3, /* Solid -Z */
- 0, 4, 2, 2, 4, 6, /* Solid -X */
- 4, 5, 6, 6, 5, 7, /* Solid +Z */
- 3, 7, 1, 1, 7, 5, /* Solid +X */
- 2, 6, 3, 3, 6, 7, /* Solid +Y */
- 0, 1, 4, 4, 1, 5, /* Solid -Y */
-
- 8, 10, 0, 0, 10, 2, /* Left fluid -Z */
- 8, 12, 10, 10, 12, 14, /* Left fluid -X */
- 12, 4, 14, 14, 4, 6, /* Left fluid +Z */
- 10, 14, 2, 2, 14, 6, /* Left fluid +Y */
- 8, 0, 12, 12, 0, 4, /* Left fluid -Y */
-
- 1, 3, 9, 9, 3, 11, /* Right fluid -Z */
- 5, 13, 7, 7, 13, 15, /* Right fluid +Z */
- 11, 15, 9, 9, 15, 13, /* Right fluid +X */
- 3, 7, 11, 11, 7, 15, /* Right fluid +Y */
- 1, 9, 5, 5, 9, 13 /* Right fluid -Y */
-};
-static const size_t nprimitives_3d = sizeof(indices_3d) / (sizeof(size_t)*3);
-
-static void
-get_indices_2d(const size_t iseg, size_t ids[2], void* ctx)
-{
- struct geometry* geom = ctx;
- CHK(ctx != NULL);
- ids[0] = geom->indices[iseg*2+0];
- ids[1] = geom->indices[iseg*2+1];
-}
-
-static void
-get_indices_3d(const size_t itri, size_t ids[3], void* ctx)
-{
- struct geometry* geom = ctx;
- CHK(ctx != NULL);
- ids[0] = geom->indices[itri*3+0];
- ids[1] = geom->indices[itri*3+1];
- ids[2] = geom->indices[itri*3+2];
-}
-
-static void
-get_position_2d(const size_t ivert, double pos[2], void* ctx)
-{
- struct geometry* geom = ctx;
- CHK(ctx != NULL);
- pos[0] = geom->positions[ivert*2+0];
- pos[1] = geom->positions[ivert*2+1];
-}
-
-static void
-get_position_3d(const size_t ivert, double pos[3], void* ctx)
-{
- struct geometry* geom = ctx;
- CHK(ctx != NULL);
- pos[0] = geom->positions[ivert*3+0];
- pos[1] = geom->positions[ivert*3+1];
- pos[2] = geom->positions[ivert*3+2];
-}
-
-static void
-get_interface(const size_t iprim, struct sdis_interface** bound, void* ctx)
-{
- struct geometry* geom = ctx;
- CHK(ctx != NULL);
- *bound = geom->interfaces[iprim];
-}
-
-/*******************************************************************************
- * media
- ******************************************************************************/
-struct solid {
- double lambda;
- double rho;
- double cp;
- double volumic_power;
-};
-
-static double
-solid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- const struct solid* solid = sdis_data_cget(data);
- CHK(vtx && solid);
- return solid->cp;
-}
-
-static double
-solid_get_thermal_conductivity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- const struct solid* solid = sdis_data_cget(data);
- CHK(vtx && solid);
- return solid->lambda;
-}
-
-static double
-solid_get_volumic_mass
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- const struct solid* solid = sdis_data_cget(data);
- CHK(vtx && solid);
- return solid->rho;
-}
-
-static double
-solid_get_delta
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx && data);
- return 0.005;
-}
-
-static double
-solid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx && data);
- return UNKNOWN_TEMPERATURE;
-}
-
-static double
-solid_get_volumic_power
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- const struct solid* solid = sdis_data_cget(data);
- CHK(vtx && solid);
- return solid->volumic_power;
-}
-
-static void
-create_solid
- (struct sdis_device* dev,
- const struct solid* solid_props,
- struct sdis_medium** solid)
-{
- struct sdis_data* data = NULL;
- struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
- CHK(dev && solid_props && solid);
-
- OK(sdis_data_create
- (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- memcpy(sdis_data_get(data), solid_props, sizeof(struct solid));
- shader.calorific_capacity = solid_get_calorific_capacity;
- shader.thermal_conductivity = solid_get_thermal_conductivity;
- shader.volumic_mass = solid_get_volumic_mass;
- shader.delta_solid = solid_get_delta;
- shader.temperature = solid_get_temperature;
- shader.volumic_power = solid_get_volumic_power;
- OK(sdis_solid_create(dev, &shader, data, solid));
- OK(sdis_data_ref_put(data));
-}
-
-static void
-create_fluid(struct sdis_device* dev, struct sdis_medium** fluid)
-{
- struct sdis_fluid_shader shader = DUMMY_FLUID_SHADER;
- OK(sdis_fluid_create(dev, &shader, NULL, fluid));
-}
-
-/*******************************************************************************
- * Interface
- ******************************************************************************/
-struct interf {
- double temperature;
- double h;
- double emissivity;
- double specular_fraction;
- double Tref;
-};
-
-static double
-interface_get_temperature
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf = sdis_data_cget(data);
- CHK(frag && interf);
- return interf->temperature;
-}
-
-static double
-interface_get_convection_coef
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf = sdis_data_cget(data);
- CHK(frag && interf);
- return interf->h;
-}
-
-static double
-interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf = sdis_data_cget(data);
- CHK(frag && interf);
- return interf->emissivity;
-}
-
-static double
-interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf = sdis_data_cget(data);
- CHK(frag && interf);
- return interf->specular_fraction;
-}
-
-static double
-interface_get_Tref
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf = sdis_data_cget(data);
- CHK(frag && interf);
- return interf->Tref;
-}
-
-static void
-create_interface
- (struct sdis_device* dev,
- struct sdis_medium* front,
- struct sdis_medium* back,
- const struct interf* interf,
- struct sdis_interface** out_interf)
-{
- struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
- struct sdis_data* data = NULL;
-
- CHK(interf != NULL);
-
- shader.front.temperature = interface_get_temperature;
- shader.back.temperature = interface_get_temperature;
- if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
- shader.convection_coef = interface_get_convection_coef;
- }
- if(sdis_medium_get_type(front) == SDIS_FLUID) {
- shader.front.emissivity = interface_get_emissivity;
- shader.front.specular_fraction = interface_get_specular_fraction;
- shader.front.reference_temperature = interface_get_Tref;
- }
- if(sdis_medium_get_type(back) == SDIS_FLUID) {
- shader.back.emissivity = interface_get_emissivity;
- shader.back.specular_fraction = interface_get_specular_fraction;
- shader.back.reference_temperature = interface_get_Tref;
- }
- shader.convection_coef_upper_bound = MMAX(0, interf->h);
-
- OK(sdis_data_create
- (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
- memcpy(sdis_data_get(data), interf, sizeof(*interf));
-
- OK(sdis_interface_create(dev, front, back, &shader, data, out_interf));
- OK(sdis_data_ref_put(data));
-}
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-struct reference_result {
- double T; /* Temperature at the center of the solid [K] */
- double T1; /* Temperature on the left boundary of the solid [K] */
- double T2; /* Temperature on the right boundary of the solid [K] */
-};
-static const struct reference_result REFERENCE_RESULT_NULL = {0,0,0};
-
-static void
-test_picard1(struct sdis_scene* scn, 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;
- struct sdis_mc mc = SDIS_MC_NULL;
- struct sdis_estimator* estimator = NULL;
- enum sdis_scene_dimension dim;
- size_t prims[2];
- enum sdis_side sides[2];
- CHK(scn);
-
- OK(sdis_scene_get_dimension(scn, &dim));
- switch(dim) {
- case SDIS_SCENE_2D: printf(">>> 2D\n"); break;
- case SDIS_SCENE_3D: printf(">>> 3D\n"); break;
- default: FATAL("Unreachable code.\n"); break;
- }
-
- probe_args.nrealisations = N;
- probe_args.position[0] = 0.05;
- probe_args.position[1] = 0;
- probe_args.position[2] = 0;
- 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",
- SPLIT3(probe_args.position), ref->T, mc.E, mc.SE);
- CHK(eq_eps(ref->T, mc.E, mc.SE*3));
- OK(sdis_estimator_ref_put(estimator));
-
- switch(dim) {
- case SDIS_SCENE_2D:
- prims[0] = 1; sides[0] = SDIS_BACK;
- bound_args.nprimitives = 1;
- break;
- case SDIS_SCENE_3D:
- prims[0] = 2; sides[0] = SDIS_BACK;
- prims[1] = 3; sides[1] = SDIS_BACK;
- bound_args.nprimitives = 2;
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- bound_args.nrealisations = N;
- bound_args.primitives = prims;
- bound_args.sides = sides;
- 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);
- CHK(eq_eps(ref->T1, mc.E, mc.SE*3));
- OK(sdis_estimator_ref_put(estimator));
-
- switch(dim) {
- case SDIS_SCENE_2D:
- prims[0] = 3; sides[0] = SDIS_BACK;
- bound_args.nprimitives = 1;
- break;
- case SDIS_SCENE_3D:
- prims[0] = 6; sides[0] = SDIS_BACK;
- prims[1] = 7; sides[1] = SDIS_BACK;
- bound_args.nprimitives = 2;
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- bound_args.nrealisations = N;
- bound_args.primitives = prims;
- bound_args.sides = sides;
- 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);
- CHK(eq_eps(ref->T2, mc.E, mc.SE*3));
- OK(sdis_estimator_ref_put(estimator));
-}
-
-static void
-create_scene_3d
- (struct sdis_device* dev,
- struct sdis_interface* interfaces[INTERFACES_COUNT__],
- struct sdis_scene** scn)
-{
- struct geometry geom;
- struct sdis_interface* prim_interfaces[32];
- struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
-
- CHK(dev && interfaces && scn);
-
- /* Setup the per primitive interface of the solid medium */
- prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC];
- prim_interfaces[2] = prim_interfaces[3] = interfaces[SOLID_FLUID_mX];
- prim_interfaces[4] = prim_interfaces[5] = interfaces[ADIABATIC];
- prim_interfaces[6] = prim_interfaces[7] = interfaces[SOLID_FLUID_pX];
- prim_interfaces[8] = prim_interfaces[9] = interfaces[ADIABATIC];
- prim_interfaces[10] = prim_interfaces[11] = interfaces[ADIABATIC];
-
- /* Setup the per primitive interface for the left fluid */
- prim_interfaces[12] = prim_interfaces[13] = interfaces[BOUNDARY_mX];
- prim_interfaces[14] = prim_interfaces[15] = interfaces[BOUNDARY_mX];
- prim_interfaces[16] = prim_interfaces[17] = interfaces[BOUNDARY_mX];
- prim_interfaces[18] = prim_interfaces[19] = interfaces[BOUNDARY_mX];
- prim_interfaces[20] = prim_interfaces[21] = interfaces[BOUNDARY_mX];
-
- /* Setup the per primitive interface for the right fluid */
- prim_interfaces[22] = prim_interfaces[23] = interfaces[BOUNDARY_pX];
- prim_interfaces[24] = prim_interfaces[25] = interfaces[BOUNDARY_pX];
- prim_interfaces[26] = prim_interfaces[27] = interfaces[BOUNDARY_pX];
- prim_interfaces[28] = prim_interfaces[29] = interfaces[BOUNDARY_pX];
- prim_interfaces[30] = prim_interfaces[31] = interfaces[BOUNDARY_pX];
-
- /* Create the scene */
- geom.positions = vertices_3d;
- geom.indices = indices_3d;
- geom.interfaces = prim_interfaces;
- scn_args.get_indices = get_indices_3d;
- scn_args.get_interface = get_interface;
- scn_args.get_position = get_position_3d;
- scn_args.nprimitives = nprimitives_3d;
- scn_args.nvertices = nvertices_3d;
- scn_args.tmax = 350;
- scn_args.context = &geom;
- OK(sdis_scene_create(dev, &scn_args, scn));
-}
-
-static void
-create_scene_2d
- (struct sdis_device* dev,
- struct sdis_interface* interfaces[INTERFACES_COUNT__],
- struct sdis_scene** scn)
-{
- struct geometry geom;
- struct sdis_interface* prim_interfaces[10/*#segment*/];
- struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
-
- CHK(dev && interfaces && scn);
-
- /* Setup the per primitive interface of the solid medium */
- prim_interfaces[0] = interfaces[ADIABATIC];
- prim_interfaces[1] = interfaces[SOLID_FLUID_mX];
- prim_interfaces[2] = interfaces[ADIABATIC];
- prim_interfaces[3] = interfaces[SOLID_FLUID_pX];
-
- /* Setup the per primitive interface of the fluid on the left of the medium */
- prim_interfaces[4] = interfaces[BOUNDARY_mX];
- prim_interfaces[5] = interfaces[BOUNDARY_mX];
- prim_interfaces[6] = interfaces[BOUNDARY_mX];
-
- /* Setup the per primitive interface of the fluid on the right of the medium */
- prim_interfaces[7] = interfaces[BOUNDARY_pX];
- prim_interfaces[8] = interfaces[BOUNDARY_pX];
- prim_interfaces[9] = interfaces[BOUNDARY_pX];
-
- /* Create the scene */
- geom.positions = vertices_2d;
- geom.indices = indices_2d;
- geom.interfaces = prim_interfaces;
- scn_args.get_indices = get_indices_2d;
- scn_args.get_interface = get_interface;
- scn_args.get_position = get_position_2d;
- scn_args.nprimitives = nprimitives_2d;
- scn_args.nvertices = nvertices_2d;
- scn_args.tmax = 350;
- scn_args.context = &geom;
- OK(sdis_scene_2d_create(dev, &scn_args, scn));
-}
-
-/*******************************************************************************
- * Test
- ******************************************************************************/
-int
-main(int argc, char** argv)
-{
- struct mem_allocator allocator;
-
- struct sdis_device* dev = NULL;
- struct sdis_scene* scn_2d = NULL;
- struct sdis_scene* scn_3d = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_medium* fluid = NULL;
- struct sdis_medium* dummy = NULL;
- struct sdis_interface* interfaces[INTERFACES_COUNT__];
-
- struct solid solid_props;
- struct solid* psolid_props;
- struct reference_result ref = REFERENCE_RESULT_NULL;
- struct interf interf_props;
- struct interf* pinterf_props[INTERFACES_COUNT__];
-
- size_t i;
- (void)argc, (void)argv;
-
- OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
- OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
-
- /* Solid medium */
- solid_props.lambda = 1.15;
- solid_props.rho = 1000;
- solid_props.cp = 800;
- create_solid(dev, &solid_props, &solid);
-
- /* Dummy solid medium */
- solid_props.lambda = 0;
- solid_props.rho = 1000;
- solid_props.cp = 800;
- solid_props.volumic_power = SDIS_VOLUMIC_POWER_NONE;
- create_solid(dev, &solid_props, &dummy);
-
- /* Fluid medium */
- create_fluid(dev, &fluid);
-
- /* Create the adiabatic interface for the solid */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
- interf_props.h = -1;
- interf_props.emissivity = -1;
- interf_props.specular_fraction = -1;
- interf_props.Tref = UNKNOWN_TEMPERATURE;
- create_interface(dev, solid, dummy, &interf_props, interfaces+ADIABATIC);
-
- /* Create the interface between the solid and the fluid */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
- interf_props.h = 0;
- interf_props.emissivity = 1;
- interf_props.specular_fraction = 0;
- interf_props.Tref = 280;
- create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_mX);
- interf_props.Tref = 350;
- create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_pX);
-
- /* Create the interface for the fluid on the left */
- interf_props.temperature = 280;
- interf_props.h = -1;
- interf_props.emissivity = 1;
- interf_props.specular_fraction = -1;
- interf_props.Tref = 280;
- create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_mX);
-
- /* Create the interface for the fluid on the right */
- interf_props.temperature = 350;
- interf_props.h = -1;
- interf_props.emissivity = 1;
- interf_props.specular_fraction = -1;
- interf_props.Tref = 350;
- create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_pX);
-
- /* Fetch pointers toward the solid and the interfaces */
- psolid_props = sdis_data_get(sdis_medium_get_data(solid));
- FOR_EACH(i, 0, INTERFACES_COUNT__) {
- pinterf_props[i] = sdis_data_get(sdis_interface_get_data(interfaces[i]));
- }
-
- create_scene_2d(dev, interfaces, &scn_2d);
- create_scene_3d(dev, interfaces, &scn_3d);
-
- /* Test picard1 with a constant Tref <=> regular linearisation */
- printf("Test Picard1 with a constant Tref of 300 K\n");
- ref.T = 314.99999999999989;
- ref.T1 = 307.64122364709766;
- ref.T2 = 322.35877635290217;
- pinterf_props[SOLID_FLUID_mX]->Tref = 300;
- 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);
- printf("\n");
-
- /* Test picard1 using T4 as a reference */
- printf("Test Picard1 using T4 as a reference\n");
- ref.T = 320.37126474482994;
- ref.T1 = 312.12650299072266;
- ref.T2 = 328.61602649893723;
- pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
- 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);
- printf("\n");
-
- /* Add volumic power */
- psolid_props->volumic_power = 1000;
-
- /* Test picard1 with a volumic power and constant Tref */
- printf("Test Picard1 with a volumic power of 1000 W/m^3 and a constant Tref "
- "of 300 K\n");
- ref.T = 324.25266420769509;
- ref.T1 = 315.80693133305368;
- ref.T2 = 330.52448403885825;
- pinterf_props[SOLID_FLUID_mX]->Tref = 300;
- 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);
- printf("\n");
-
- /* Test picard1 with a volumic power and T4 a the reference */
- printf("Test Picard1 with a volumic power of 1000 W/m^3 and T4 as a reference\n");
- ref.T = 327.95981050850446;
- ref.T1 = 318.75148773193359;
- ref.T2 = 334.99422024159708;
- pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
- 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);
- printf("\n");
-
- /* Release memory */
- OK(sdis_scene_ref_put(scn_2d));
- OK(sdis_scene_ref_put(scn_3d));
- OK(sdis_interface_ref_put(interfaces[ADIABATIC]));
- OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_mX]));
- OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_pX]));
- OK(sdis_interface_ref_put(interfaces[BOUNDARY_mX]));
- OK(sdis_interface_ref_put(interfaces[BOUNDARY_pX]));
- OK(sdis_medium_ref_put(fluid));
- OK(sdis_medium_ref_put(solid));
- OK(sdis_medium_ref_put(dummy));
- OK(sdis_device_ref_put(dev));
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHK(mem_allocated_size() == 0);
- return 0;
-}
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -139,6 +139,9 @@ 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;
@@ -261,8 +264,9 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
struct sdis_ambient_radiative_temperature trad =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
double lower[2], upper[2];
+ double t_range[2];
double u0, u1, u2, pos[2], pos1[2];
- double dst, fp, t;
+ double dst, fp;
struct context ctx;
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
@@ -270,6 +274,7 @@ 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;
@@ -304,6 +309,9 @@ 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;
@@ -369,19 +377,35 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
CHK(trad.temperature == 100);
CHK(trad.reference == 110);
- BA(sdis_scene_get_maximum_temperature(NULL, NULL));
- BA(sdis_scene_get_maximum_temperature(scn, NULL));
- BA(sdis_scene_get_maximum_temperature(NULL, &t));
- OK(sdis_scene_get_maximum_temperature(scn, &t));
- CHK(t == SDIS_SCENE_CREATE_ARGS_DEFAULT.tmax);
-
- t = -1;
- BA(sdis_scene_set_maximum_temperature(NULL, t));
- BA(sdis_scene_set_maximum_temperature(scn, t));
- t = 100;
- OK(sdis_scene_set_maximum_temperature(scn, t));
- OK(sdis_scene_get_maximum_temperature(scn, &t));
- CHK(t == 100);
+ BA(sdis_scene_get_temperature_range(NULL, NULL));
+ BA(sdis_scene_get_temperature_range(scn, NULL));
+ BA(sdis_scene_get_temperature_range(NULL, t_range));
+ OK(sdis_scene_get_temperature_range(scn, t_range));
+ CHK(t_range[0] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]);
+ CHK(t_range[1] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[1]);
+
+ t_range[0] = 1;
+ t_range[1] = 100;
+
+ BA(sdis_scene_set_temperature_range(NULL, t_range));
+ BA(sdis_scene_set_temperature_range(scn, NULL));
+ OK(sdis_scene_set_temperature_range(scn, t_range));
+ t_range[0] = -1;
+ t_range[1] = -1;
+ OK(sdis_scene_get_temperature_range(scn, t_range));
+ 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));
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -347,7 +347,8 @@ main(int argc, char** argv)
scn_args.nvertices = box_nvertices;
scn_args.trad.temperature = Trad;
scn_args.trad.reference = Trad;
- scn_args.tmax = MMAX(MMAX(Tf, Trad), Tb);
+ scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb);
+ scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb);
scn_args.context = box_interfaces;
OK(sdis_scene_create(dev, &scn_args, &box_scn));
@@ -359,7 +360,8 @@ main(int argc, char** argv)
scn_args.nvertices = square_nvertices;
scn_args.trad.temperature = Trad;
scn_args.trad.reference = Trad;
- scn_args.tmax = MMAX(MMAX(Tf, Trad), Tb);
+ scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb);
+ scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb);
scn_args.context = square_interfaces;
OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -634,7 +634,8 @@ main(int argc, char** argv)
scn_args.nvertices = npos;
scn_args.trad.temperature = 300;
scn_args.trad.reference = 300;
- scn_args.tmax = 350;
+ scn_args.t_range[0] = 300;
+ scn_args.t_range[1] = 350;
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, &scn));
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -286,6 +286,7 @@ main(int argc, char** argv)
struct ssp_rng* rng_state = NULL;
enum sdis_estimator_type type;
FILE* stream = NULL;
+ double t_range[2];
double ref;
const size_t N = 1000;
const size_t N_dump = 10;
@@ -551,8 +552,10 @@ main(int argc, char** argv)
/* Green and ambient radiative temperature */
solve_args.nrealisations = N;
trad.temperature = trad.reference = 300;
+ t_range[0] = 300;
+ t_range[1] = 300;
OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
- OK(sdis_scene_set_maximum_temperature(scn, 300));
+ OK(sdis_scene_set_temperature_range(scn, t_range));
interface_param->epsilon = 1;
interface_param->reference_temperature = 300;
@@ -569,8 +572,10 @@ main(int argc, char** argv)
/* Check same green used at different ambient radiative temperature */
trad.temperature = 600;
+ t_range[0] = 300;
+ t_range[1] = 600;
OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
- OK(sdis_scene_set_maximum_temperature(scn, 600));
+ OK(sdis_scene_set_temperature_range(scn, t_range));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_green_function_solve(green, &estimator2));
diff --git a/src/test_sdis_unstationary_atm.c b/src/test_sdis_unstationary_atm.c
@@ -89,7 +89,7 @@
#define X_PROBE (XH + 0.2 * XE)
-#define DELTA (XE/30.0)
+#define DELTA (XE/40.0)
/*******************************************************************************
* Box geometry
@@ -114,11 +114,11 @@ static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)
* triangle.
* ,3---,4---,5 ,3----4----5 ,4
* ,' | ,' | ,'/| ,'/| \ | \ | ,'/|
- * 9----10---11 / | 9' / | \ | \ | 10 / | Y
- * |', |', | / ,2 | / ,0---,1---,2 | / ,1 |
+ * 9----10---11 / | 9' / | \ | \ | 10 / | Y
+ * |', |', | / ,2 | / ,0---,1---,2 | / ,1 |
* | ',| ',|/,' |/,' | ,' | ,' |/,' o--X
- * 6----7----8' 6----7'---8' 7 /
- * Front, right Back, left and Internal Z
+ * 6----7----8' 6----7'---8' 7 /
+ * Front, right Back, left and Internal Z
* and Top faces bottom faces face */
static const size_t model3d_indices[22/*#triangles*/*3/*#indices per triangle*/] = {
0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */
@@ -421,10 +421,11 @@ solve_tbound1
size_t nreals;
size_t nfails;
enum sdis_scene_dimension dim;
- double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- double ref[sizeof(t) / sizeof(*t)]
- = { 290.046375, 289.903935, 289.840490, 289.802690, 289.777215,
- 289.759034, 289.745710, 289.735826, 289.728448, 289.722921 };
+ const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 290.046375, 289.903935, 289.840490, 289.802690, 289.777215, 289.759034,
+ 289.745710, 289.735826, 289.728448, 289.722921
+ };
const int nsimuls = sizeof(t) / sizeof(*t);
int isimul;
ASSERT(scn && rng);
@@ -475,9 +476,8 @@ solve_tbound1
printf("Elapsed time = %s\n", dump);
printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
- CHK(nfails + nreals == N);
- CHK(nfails <= N/1000);
CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
OK(sdis_estimator_ref_put(estimator));
}
@@ -498,10 +498,11 @@ solve_tbound2
size_t nreals;
size_t nfails;
enum sdis_scene_dimension dim;
- double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- double ref[sizeof(t) / sizeof(*t)]
- = { 309.08032, 309.34626, 309.46525, 309.53625, 309.58408,
- 309.618121, 309.642928, 309.661167, 309.674614, 309.684524 };
+ const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 309.08032, 309.34626, 309.46525, 309.53625, 309.58408, 309.618121,
+ 309.642928, 309.661167, 309.674614, 309.684524
+ };
const int nsimuls = sizeof(t) / sizeof(*t);
int isimul;
ASSERT(scn && rng);
@@ -555,6 +556,7 @@ solve_tbound2
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
OK(sdis_estimator_ref_put(estimator));
}
@@ -574,10 +576,11 @@ solve_tsolid
size_t nreals;
size_t nfails;
enum sdis_scene_dimension dim;
- double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- double ref[sizeof(t) / sizeof(*t)]
- = { 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030,
- 304.75041, 305.01595, 305.21193, 305.35641, 305.46271 };
+ const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030, 304.75041,
+ 305.01595, 305.21193, 305.35641, 305.46271
+ };
const int nsimuls = sizeof(t) / sizeof(*t);
int isimul;
ASSERT(scn && rng);
@@ -621,6 +624,7 @@ solve_tsolid
CHK(nfails + nreals == N);
CHK(nfails <= N / 1000);
CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
OK(sdis_estimator_ref_put(estimator));
}
@@ -640,10 +644,11 @@ solve_tfluid
size_t nfails;
enum sdis_scene_dimension dim;
double eps;
- double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- double ref[sizeof(t) / sizeof(*t)]
- = { 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899,
- 309.82141, 309.83055, 309.83728, 309.84224 };
+ const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899,
+ 309.82141, 309.83055, 309.83728, 309.84224
+ };
const int nsimuls = sizeof(t) / sizeof(*t);
int isimul;
ASSERT(scn);
@@ -834,7 +839,7 @@ main(int argc, char** argv)
model3d_interfaces[10] = interf_TA;
model3d_interfaces[11] = interf_TA;
/* Top */
- model3d_interfaces[12] = interf_adiabatic_1;
+ model3d_interfaces[12] = interf_adiabatic_1;
model3d_interfaces[13] = interf_adiabatic_1;
model3d_interfaces[14] = interf_adiabatic_2;
model3d_interfaces[15] = interf_adiabatic_2;
@@ -869,10 +874,8 @@ main(int argc, char** argv)
scn_args.context = model3d_interfaces;
scn_args.trad.temperature = TR;
scn_args.trad.reference = TR;
- scn_args.tmax = MMAX(T0_FLUID, T0_SOLID);
- scn_args.tmax = MMAX(scn_args.tmax, TA);
- scn_args.tmax = MMAX(scn_args.tmax, TG);
- scn_args.tmax = MMAX(scn_args.tmax, TR);
+ scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
+ scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
OK(sdis_scene_create(dev, &scn_args, &box_scn));
/* Create the square scene */
@@ -884,10 +887,8 @@ main(int argc, char** argv)
scn_args.context = model2d_interfaces;
scn_args.trad.temperature = TR;
scn_args.trad.reference = TR;
- scn_args.tmax = MMAX(T0_FLUID, T0_SOLID);
- scn_args.tmax = MMAX(scn_args.tmax, TA);
- scn_args.tmax = MMAX(scn_args.tmax, TG);
- scn_args.tmax = MMAX(scn_args.tmax, TR);
+ scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
+ scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
/* Release the interfaces */
