commit f7d4355e899579b82cef12cfda227281ea90ce1d
parent 6228fa33d8f0bb9eac3599d4cdf02462d9c7fa1b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 2 Nov 2021 15:19:20 +0100
Merge branch 'feature_picard1' into develop
Diffstat:
39 files changed, 2779 insertions(+), 1324 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -68,6 +68,7 @@ set(SDIS_FILES_SRC
sdis_estimator_buffer.c
sdis_green.c
sdis_heat_path.c
+ sdis_heat_path_boundary.c
sdis_interface.c
sdis_log.c
sdis_medium.c
@@ -85,7 +86,11 @@ set(SDIS_FILES_INC
sdis_estimator_c.h
sdis_green.h
sdis_heat_path.h
+ sdis_heat_path_boundary_c.h
sdis_heat_path_boundary_Xd.h
+ sdis_heat_path_boundary_Xd_fixed_flux.h
+ sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
+ sdis_heat_path_boundary_Xd_solid_solid.h
sdis_heat_path_conductive_Xd.h
sdis_heat_path_convective_Xd.h
sdis_heat_path_radiative_Xd.h
@@ -187,6 +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_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
@@ -207,8 +207,11 @@ struct sdis_interface_side_shader {
* interface or if the emissivity is 0 onto the whole interface. */
sdis_interface_getter_T emissivity; /* Overall emissivity. */
sdis_interface_getter_T specular_fraction; /* Specular part in [0,1] */
+
+ /* Reference temperature used in Picard 1 */
+ sdis_interface_getter_T reference_temperature;
};
-#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL }
+#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL, NULL }
static const struct sdis_interface_side_shader SDIS_INTERFACE_SIDE_SHADER_NULL =
SDIS_INTERFACE_SIDE_SHADER_NULL__;
@@ -357,20 +360,29 @@ typedef void
double pos[], /* Output list of vertex coordinates */
void* ctx);
+struct sdis_ambient_radiative_temperature {
+ double temperature; /* In Kelvin */
+ double reference; /* Used to linearise the radiative transfert */
+};
+#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ {-1, -1}
+static const struct sdis_ambient_radiative_temperature
+SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL =
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__;
+
struct sdis_scene_create_args {
/* Functors to retrieve the geometric description */
sdis_get_primitive_indices_T get_indices;
sdis_get_primitive_interface_T get_interface;
sdis_get_vertex_position_T get_position;
- /* Pointer toward client side sent as the last argument of the callbacks */
+ /* Pointer toward client side sent as the last argument of the callbacks */
void* context;
size_t nprimitives; /* #primitives, i.e. #segments or #triangles */
size_t nvertices; /* #vertices */
double fp_to_meter; /* Scale factor used to convert 1.0 in 1 meter */
- double trad; /* Ambiant radiative temperature */
- double tref; /* Temperature used to linearize the radiative temperature */
+ struct sdis_ambient_radiative_temperature trad; /* Ambient radiative temp */
+ double tmax; /* Max temperature used to linearize the radiative temperature */
};
#define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \
@@ -381,8 +393,8 @@ struct sdis_scene_create_args {
0, /* #primitives */ \
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
- -1.0, /* Ambient radiative temperature */ \
- -1.0 /* Reference temperature */ \
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
+ -1.0, /* Maximum temperature */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
SDIS_SCENE_CREATE_ARGS_DEFAULT__;
@@ -826,27 +838,27 @@ sdis_scene_set_fp_to_meter
SDIS_API res_T
sdis_scene_get_ambient_radiative_temperature
(const struct sdis_scene* scn,
- double* trad);
+ struct sdis_ambient_radiative_temperature* trad);
/* Set scene's ambient radiative temperature. If set negative, any sample
* ending in ambient radiative temperature will fail */
SDIS_API res_T
sdis_scene_set_ambient_radiative_temperature
(struct sdis_scene* scn,
- const double trad);
+ const struct sdis_ambient_radiative_temperature* trad);
-/* Get scene's reference temperature */
+/* Get scene's maximum temperature */
SDIS_API res_T
-sdis_scene_get_reference_temperature
+sdis_scene_get_maximum_temperature
(const struct sdis_scene* scn,
- double* tref);
+ double* tmax);
-/* Set scene's reference temperature. If set to 0, there is no radiative
- * transfert in the whole system */
+/* Set scene's maximum temperature. Must be correctly defined if there is any
+ * radiative transfert in the scene. */
SDIS_API res_T
-sdis_scene_set_reference_temperature
+sdis_scene_set_maximum_temperature
(struct sdis_scene* scn,
- const double tref);
+ const double tmax);
/* 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
@@ -25,10 +25,17 @@ struct sdis_heat_path;
struct rwalk_context {
struct green_path_handle* green_path;
struct sdis_heat_path* heat_path;
- double Tarad; /* Ambient radiative temperature */
- double Tref3; /* Reference temperature ^ 3 */
+
+ /* That is the upper bound temperature */
+ double That2; /* That^2 */
+ double That3; /* That^3 */
};
-#define RWALK_CONTEXT_NULL__ {NULL, NULL, 0, 0}
+#define RWALK_CONTEXT_NULL__ { \
+ NULL, /* Green path */ \
+ NULL, /* Heat path */ \
+ 0, /* (Temperature upper bound)^2 */ \
+ 0 /* (Temperature upper bound)^3 */ \
+}
static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
#endif /* SDIS_XD_BEGIN_H */
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -387,6 +387,8 @@ green_function_solve_path
const size_t ipath,
double* weight)
{
+ struct sdis_ambient_radiative_temperature trad =
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
const struct power_term* power_terms = NULL;
const struct flux_term* flux_terms = NULL;
const struct green_path* path = NULL;
@@ -443,7 +445,8 @@ green_function_solve_path
break;
case SDIS_GREEN_PATH_END_RADIATIVE:
SDIS(green_function_get_scene(green, &scn));
- SDIS(scene_get_ambient_radiative_temperature(scn, &end_temperature));
+ SDIS(scene_get_ambient_radiative_temperature(scn, &trad));
+ end_temperature = trad.temperature;
if(end_temperature < 0) { /* Cannot be negative if used */
res = RES_BAD_ARG;
goto error;
@@ -1549,6 +1552,15 @@ green_path_set_limit_radiative
}
res_T
+green_path_reset_limit(struct green_path_handle* handle)
+{
+ ASSERT(handle);
+ handle->path->elapsed_time = -INF;
+ handle->path->end_type = SDIS_GREEN_PATH_END_TYPES_COUNT__;
+ return RES_OK;
+}
+
+res_T
green_path_add_power_term
(struct green_path_handle* handle,
struct sdis_medium* mdm,
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -87,6 +87,10 @@ green_path_set_limit_radiative
const double elapsed_time);
extern LOCAL_SYM res_T
+green_path_reset_limit
+ (struct green_path_handle* handle);
+
+extern LOCAL_SYM res_T
green_path_add_power_term
(struct green_path_handle* path,
struct sdis_medium* mdm,
diff --git a/src/sdis_heat_path.c b/src/sdis_heat_path.c
@@ -33,12 +33,6 @@
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_conductive_Xd.h"
-/* Generate the boundary path routines */
-#define SDIS_XD_DIMENSION 2
-#include "sdis_heat_path_boundary_Xd.h"
-#define SDIS_XD_DIMENSION 3
-#include "sdis_heat_path_boundary_Xd.h"
-
/*******************************************************************************
* Exported functions
******************************************************************************/
diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h
@@ -19,6 +19,7 @@
#include "sdis.h"
#include <rsys/dynamic_array.h>
+#include <rsys/dynamic_array_size_t.h>
#include <rsys/rsys.h>
/* Forward declarations */
@@ -39,7 +40,12 @@ struct temperature_3d;
* Heat path data structure
******************************************************************************/
struct sdis_heat_path {
+ /* List of the path vertices */
struct darray_heat_vertex vertices;
+
+ /* Indices of the vertices that mark a break in the path */
+ struct darray_size_t breaks;
+
enum sdis_heat_path_flag status;
};
@@ -49,6 +55,7 @@ heat_path_init(struct mem_allocator* allocator, struct sdis_heat_path* path)
ASSERT(path);
path->status = SDIS_HEAT_PATH_NONE;
darray_heat_vertex_init(allocator, &path->vertices);
+ darray_size_t_init(allocator, &path->breaks);
}
static INLINE void
@@ -56,30 +63,46 @@ heat_path_release(struct sdis_heat_path* path)
{
ASSERT(path);
darray_heat_vertex_release(&path->vertices);
+ darray_size_t_release(&path->breaks);
}
static INLINE res_T
heat_path_copy(struct sdis_heat_path* dst, const struct sdis_heat_path* src)
{
+ res_T res = RES_OK;
ASSERT(dst && src);
dst->status = src->status;
- return darray_heat_vertex_copy(&dst->vertices, &src->vertices);
+ res = darray_heat_vertex_copy(&dst->vertices, &src->vertices);
+ if(res != RES_OK) return res;
+ res = darray_size_t_copy(&dst->breaks, &src->breaks);
+ if(res != RES_OK) return res;
+ return RES_OK;
}
static INLINE res_T
heat_path_copy_and_release(struct sdis_heat_path* dst, struct sdis_heat_path* src)
{
+ res_T res = RES_OK;
ASSERT(dst && src);
dst->status = src->status;
- return darray_heat_vertex_copy_and_release(&dst->vertices, &src->vertices);
+ res = darray_heat_vertex_copy_and_release(&dst->vertices, &src->vertices);
+ if(res != RES_OK) return res;
+ res = darray_size_t_copy_and_release(&dst->breaks, &src->breaks);
+ if(res != RES_OK) return res;
+ return RES_OK;
}
static INLINE res_T
heat_path_copy_and_clear(struct sdis_heat_path* dst, struct sdis_heat_path* src)
{
+ res_T res = RES_OK;
ASSERT(dst && src);
dst->status = src->status;
- return darray_heat_vertex_copy_and_clear(&dst->vertices, &src->vertices);
+ res = darray_heat_vertex_copy_and_clear(&dst->vertices, &src->vertices);
+ if(res != RES_OK) return res;
+ res = darray_size_t_copy_and_clear(&dst->breaks, &src->breaks);
+ if(res != RES_OK) return res;
+ return RES_OK;
}
static INLINE res_T
@@ -99,6 +122,18 @@ heat_path_get_last_vertex(struct sdis_heat_path* path)
return darray_heat_vertex_data_get(&path->vertices) + (sz-1);
}
+static INLINE res_T
+heat_path_add_break(struct sdis_heat_path* path)
+{
+ size_t id;
+ size_t sz;
+ ASSERT(path);
+ sz = darray_heat_vertex_size_get(&path->vertices);
+ if(sz == 0) return RES_OK; /* Nothing to do */
+ id = sz-1;
+ return darray_size_t_push_back(&path->breaks, &id);
+}
+
/* Generate the dynamic array of heat paths */
#define DARRAY_NAME heat_path
#define DARRAY_DATA struct sdis_heat_path
diff --git a/src/sdis_heat_path_boundary.c b/src/sdis_heat_path_boundary.c
@@ -0,0 +1,43 @@
+/* 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_heat_path.h"
+#include "sdis_heat_path_boundary_c.h"
+
+/* Generate the helper routines */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd_c.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd_c.h"
+
+/* Generate the boundary path sub-routines */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd_fixed_flux.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd_fixed_flux.h"
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd_solid_fluid_picard1.h"
+#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_solid.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd_solid_solid.h"
+
+/* Generate the boundary path routines */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd.h"
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -16,955 +16,14 @@
#include "sdis_device_c.h"
#include "sdis_green.h"
#include "sdis_heat_path.h"
+#include "sdis_heat_path_boundary_c.h"
#include "sdis_interface_c.h"
#include "sdis_medium_c.h"
-#include "sdis_scene_c.h"
#include <star/ssp.h>
#include "sdis_Xd_begin.h"
-/* Emperical scale factor applied to the challenged reinjection distance. If
- * the distance to reinject is less than this adjusted value, the solver will
- * try to discard the reinjection distance if possible in order to avoid
- * numerical issues. */
-#define REINJECT_DST_MIN_SCALE 0.125f
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-XD(sample_reinjection_dir)
- (const struct XD(rwalk)* rwalk, struct ssp_rng* rng, float dir[DIM])
-{
-#if DIM == 2
- /* The sampled directions is defined by rotating the normal around the Z axis
- * of an angle of PI/4 or -PI/4. Let the rotation matrix defined as
- * | cos(a) -sin(a) |
- * | sin(a) cos(a) |
- * with a = PI/4, dir = sqrt(2)/2 * | 1 -1 | . N
- * | 1 1 |
- * with a =-PI/4, dir = sqrt(2)/2 * | 1 1 | . N
- * |-1 1 |
- * Note that since the sampled direction is finally normalized, we can
- * discard the sqrt(2)/2 constant. */
- const uint64_t r = ssp_rng_uniform_uint64(rng, 0, 1);
- ASSERT(rwalk && dir);
- if(r) {
- dir[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
- dir[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
- } else {
- dir[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
- dir[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
- }
- f2_normalize(dir, dir);
-#else
- /* Sample a random direction around the normal whose cosine is 1/sqrt(3). To
- * do so we sample a position onto a cone whose height is 1/sqrt(2) and the
- * radius of its base is 1. */
- float frame[9];
- ASSERT(fX(is_normalized)(rwalk->hit.normal));
-
- ssp_ran_circle_uniform_float(rng, dir, NULL);
- dir[2] = (float)(1.0/sqrt(2));
-
- f33_basis(frame, rwalk->hit.normal);
- f33_mulf3(dir, frame, dir);
- f3_normalize(dir, dir);
- ASSERT(eq_epsf(f3_dot(dir, rwalk->hit.normal), (float)(1.0/sqrt(3)), 1.e-4f));
-#endif
-}
-
-#if DIM == 2
-static void
-XD(move_away_primitive_boundaries)
- (struct XD(rwalk)* rwalk,
- const double delta)
-{
- struct sXd(attrib) attr;
- float pos[DIM];
- float dir[DIM];
- float len;
- const float st = 0.5f;
- ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->hit) && delta > 0);
-
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr));
-
- fX_set_dX(pos, rwalk->vtx.P);
- fX(sub)(dir, attr.value, pos);
- len = fX(normalize)(dir, dir);
- len = MMIN(len, (float)(delta*0.1));
-
- XD(move_pos)(rwalk->vtx.P, dir, len);
-}
-#else
-/* Move the random walk away from the primitive boundaries to avoid numerical
- * issues leading to inconsistent random walks. */
-static void
-XD(move_away_primitive_boundaries)
- (struct XD(rwalk)* rwalk,
- const double delta)
-{
- struct s3d_attrib v0, v1, v2; /* Triangle vertices */
- float E[3][4]; /* 3D edge equations */
- float dst[3]; /* Distance from current position to edge equation */
- float N[3]; /* Triangle normal */
- float P[3]; /* Random walk position */
- float tmp[3];
- float min_dst, max_dst;
- float cos_a1, cos_a2;
- float len;
- int imax = 0;
- int imin = 0;
- int imid = 0;
- int i;
- ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->hit));
-
- fX_set_dX(P, rwalk->vtx.P);
-
- /* Fetch triangle vertices */
- S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 0, S3D_POSITION, &v0));
- S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 1, S3D_POSITION, &v1));
- S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 2, S3D_POSITION, &v2));
-
- /* Compute the edge vector */
- f3_sub(E[0], v1.value, v0.value);
- f3_sub(E[1], v2.value, v1.value);
- f3_sub(E[2], v0.value, v2.value);
-
- /* Compute the triangle normal */
- f3_cross(N, E[1], E[0]);
-
- /* Compute the 3D edge equation */
- f3_normalize(E[0], f3_cross(E[0], E[0], N));
- f3_normalize(E[1], f3_cross(E[1], E[1], N));
- f3_normalize(E[2], f3_cross(E[2], E[2], N));
- E[0][3] = -f3_dot(E[0], v0.value);
- E[1][3] = -f3_dot(E[1], v1.value);
- E[2][3] = -f3_dot(E[2], v2.value);
-
- /* Compute the distance from current position to the edges */
- dst[0] = f3_dot(E[0], P) + E[0][3];
- dst[1] = f3_dot(E[1], P) + E[1][3];
- dst[2] = f3_dot(E[2], P) + E[2][3];
-
- /* Retrieve the min and max distance from random walk position to triangle
- * edges */
- min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
- max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
-
- /* Sort the edges with respect to their distance to the random walk position */
- FOR_EACH(i, 0, 3) {
- if(dst[i] == min_dst) {
- imin = i;
- } else if(dst[i] == max_dst) {
- imax = i;
- } else {
- imid = i;
- }
- }
- (void)imax;
-
- /* TODO if the current position is near a vertex, one should move toward the
- * farthest edge along its normal to avoid too small displacement */
-
- /* Compute the distance `dst' from the current position to the edges to move
- * to, along the normal of the edge from which the random walk is the nearest
- *
- * +. cos(a) = d / dst => dst = d / cos_a
- * / `*.
- * / | `*.
- * / dst| a /`*.
- * / | / `*.
- * / | / d `*.
- * / |/ `*.
- * +---------o----------------+ */
- cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
- cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
- dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
- dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
-
- /* Compute the maximum displacement distance into the triangle along the
- * normal of the edge from which the random walk is the nearest */
- len = MMIN(dst[imid], dst[imax]);
- ASSERT(len != FLT_MAX);
-
- /* Define the displacement distance as the minimum between 10 percent of
- * delta and len / 2. */
- len = MMIN(len*0.5f, (float)(delta*0.1));
- XD(move_pos)(rwalk->vtx.P, E[imin], len);
-}
-#endif
-
-static res_T
-XD(select_reinjection_dir)
- (const struct sdis_scene* scn,
- const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
- struct XD(rwalk)* rwalk, /* Current random walk state */
- const float dir0[DIM], /* Challenged direction */
- const float dir1[DIM], /* Challanged direction */
- const double delta, /* Max reinjection distance */
- float reinject_dir[DIM], /* Selected direction */
- float* reinject_dst, /* Effective reinjection distance */
- int can_move, /* Define of the random wal pos can be moved or not */
- int* move_pos, /* Define if the current random walk was moved. May be NULL */
- struct sXd(hit)* reinject_hit) /* Hit along the reinjection dir */
-{
- struct sdis_interface* interf;
- struct sdis_medium* mdm0;
- struct sdis_medium* mdm1;
- struct hit_filter_data filter_data;
- struct sXd(hit) hit;
- struct sXd(hit) hit0;
- struct sXd(hit) hit1;
- double tmp[DIM];
- double dst;
- double dst0;
- double dst1;
- const double delta_adjusted = delta * RAY_RANGE_MAX_SCALE;
- const float* dir;
- const float reinject_threshold = (float)delta * REINJECT_DST_MIN_SCALE;
- float org[DIM];
- float range[2];
- enum sdis_side side;
- int iattempt = 0;
- const int MAX_ATTEMPTS = can_move ? 2 : 1;
- res_T res = RES_OK;
- ASSERT(scn && mdm && rwalk && dir0 && dir1 && delta > 0);
- ASSERT(reinject_dir && reinject_dst && reinject_hit);
-
- if(move_pos) *move_pos = 0;
-
- do {
- f2(range, 0, FLT_MAX);
- fX_set_dX(org, rwalk->vtx.P);
- filter_data.XD(hit) = rwalk->hit;
- filter_data.epsilon = delta * 0.01;
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, &filter_data, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, &filter_data, &hit1));
-
- /* Retrieve the medium at the reinjection pos along dir0 */
- if(SXD_HIT_NONE(&hit0)) {
- XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir0, (float)delta);
- res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0);
- if(res == RES_BAD_OP) { mdm0 = NULL; res = RES_OK; }
- if(res != RES_OK) goto error;
- } else {
- interf = scene_get_interface(scn, hit0.prim.prim_id);
- side = fX(dot)(dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
- mdm0 = interface_get_medium(interf, side);
- }
-
- /* Retrieve the medium at the reinjection pos along dir1 */
- if(SXD_HIT_NONE(&hit1)) {
- XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir1, (float)delta);
- res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1);
- if(res == RES_BAD_OP) { mdm1 = NULL; res = RES_OK; }
- if(res != RES_OK) goto error;
- } else {
- interf = scene_get_interface(scn, hit1.prim.prim_id);
- side = fX(dot)(dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
- mdm1 = interface_get_medium(interf, side);
- }
-
- dst0 = dst1 = -1;
- if(mdm0 == mdm) { /* Check reinjection consistency */
- if(hit0.distance <= delta_adjusted) {
- dst0 = hit0.distance;
- } else {
- dst0 = delta;
- hit0 = SXD_HIT_NULL;
- }
- }
- if(mdm1 == mdm) {/* Check reinjection consistency */
- if(hit1.distance <= delta_adjusted) {
- dst1 = hit1.distance;
- } else {
- dst1 = delta;
- hit1 = SXD_HIT_NULL;
- }
- }
-
- /* No valid reinjection. Maybe the random walk is near a sharp corner and
- * thus the ray-tracing misses the enclosure geometry. Another possibility
- * is that the random walk lies roughly on an edge. In this case, sampled
- * reinjecton dirs can intersect the primitive on the other side of the
- * edge. Normally, this primitive should be filtered by the "hit_filter"
- * function but this may be not the case due to a "threshold effect". In
- * both situations, try to slightly move away from the primitive boundaries
- * and retry to find a valid reinjection. */
- if(dst0 == -1 && dst1 == -1) {
- XD(move_away_primitive_boundaries)(rwalk, delta);
- if(move_pos) *move_pos = 1;
- }
- } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS);
-
- if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */
- log_warn(scn->dev, "%s: no valid reinjection direction at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
-
- if(dst0 == -1) {
- /* Invalid dir0 -> move along dir1 */
- dir = dir1;
- dst = dst1;
- hit = hit1;
- } else if(dst1 == -1) {
- /* Invalid dir1 -> move along dir0 */
- dir = dir0;
- dst = dst0;
- hit = hit0;
- } else if(dst0 < reinject_threshold && dst1 < reinject_threshold) {
- /* The displacement along dir0 and dir1 are both below the reinjection
- * threshold that defines a distance under which the temperature gradients
- * are ignored. Move along the direction that allows the maximum
- * displacement. */
- if(dst0 > dst1) {
- dir = dir0;
- dst = dst0;
- hit = hit0;
- } else {
- dir = dir1;
- dst = dst1;
- hit = hit1;
- }
- } else if(dst0 < reinject_threshold) {
- /* Ingore dir0 that is bellow the reinject threshold */
- dir = dir1;
- dst = dst1;
- hit = hit1;
- } else if(dst1 < reinject_threshold) {
- /* Ingore dir1 that is bellow the reinject threshold */
- dir = dir0;
- dst = dst0;
- hit = hit0;
- } else {
- /* All reinjection directions are valid. Choose the first 1 that was
- * randomly selected by the sample_reinjection_dir procedure and adjust
- * the displacement distance. */
- dir = dir0;
-
- /* Define the reinjection distance along dir0 and its corresponding hit */
- if(dst0 <= dst1) {
- dst = dst0;
- hit = hit0;
- } else {
- dst = dst1;
- hit = SXD_HIT_NULL;
- }
-
- /* If the displacement distance is too close of a boundary, move to the
- * boundary in order to avoid numerical uncertainty. */
- if(!SXD_HIT_NONE(&hit0)
- && dst0 != dst
- && eq_eps(dst0, dst, dst0*0.1)) {
- dst = dst0;
- hit = hit0;
- }
- }
-
- /* Setup output variable */
- fX(set)(reinject_dir, dir);
- *reinject_dst = (float)dst;
- *reinject_hit = hit;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-XD(select_reinjection_dir_and_check_validity)
- (const struct sdis_scene* scn,
- const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
- struct XD(rwalk)* rwalk, /* Current random walk state */
- const float dir0[DIM], /* Challenged direction */
- const float dir1[DIM], /* Challanged direction */
- const double delta, /* Max reinjection distance */
- float out_reinject_dir[DIM], /* Selected direction */
- float* out_reinject_dst, /* Effective reinjection distance */
- int can_move, /* Define of the random wal pos can be moved or not */
- int* move_pos, /* Define if the current random walk was moved. May be NULL */
- int* is_valid, /* Define if the reinjection defines a valid pos */
- struct sXd(hit)* out_reinject_hit) /* Hit along the reinjection dir */
-{
- double pos[DIM];
- struct sdis_medium* reinject_mdm;
- struct sXd(hit) reinject_hit;
- float reinject_dir[DIM];
- float reinject_dst;
- res_T res = RES_OK;
- ASSERT(is_valid && out_reinject_dir && out_reinject_dst && out_reinject_hit);
-
- /* Select a reinjection direction */
- res = XD(select_reinjection_dir)(scn, mdm, rwalk, dir0, dir1, delta,
- reinject_dir, &reinject_dst, can_move, move_pos, &reinject_hit);
- if(res != RES_OK) goto error;
-
- if(!SXD_HIT_NONE(&reinject_hit)) {
- *is_valid = 1;
- } else {
- /* Check medium consistency at the reinjection position */
- XD(move_pos)(dX(set)(pos, rwalk->vtx.P), reinject_dir, reinject_dst);
- res = scene_get_medium_in_closed_boundaries
- (scn, pos, &reinject_mdm);
- if(res == RES_BAD_OP) { reinject_mdm = NULL; res = RES_OK; }
- if(res != RES_OK) goto error;
-
- *is_valid = reinject_mdm == mdm;
- }
-
- if(*is_valid) {
- fX(set)(out_reinject_dir, reinject_dir);
- *out_reinject_dst = reinject_dst;
- *out_reinject_hit = reinject_hit;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/* Check that the interface fragment is consistent with the current state of
- * the random walk */
-static INLINE int
-XD(check_rwalk_fragment_consistency)
- (const struct XD(rwalk)* rwalk,
- const struct sdis_interface_fragment* frag)
-{
- double N[DIM];
- double uv[2] = {0, 0};
- ASSERT(rwalk && frag);
- dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
- if( SXD_HIT_NONE(&rwalk->hit)
- || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
- || !dX(eq_eps)(N, frag->Ng, 1.e-6)
- || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
- || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
- return 0;
- }
-#if (SDIS_XD_DIMENSION == 2)
- uv[0] = rwalk->hit.u;
-#else
- d2_set_f2(uv, rwalk->hit.uv);
-#endif
- return d2_eq_eps(uv, frag->uv, 1.e-6);
-}
-
-static res_T
-XD(solid_solid_boundary_path)
- (const struct sdis_scene* scn,
- const struct rwalk_context* ctx,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- struct sXd(hit) hit0, hit1;
- struct sXd(hit)* hit;
- struct XD(rwalk) rwalk_saved;
- struct sdis_interface* interf = NULL;
- struct sdis_medium* solid_front = NULL;
- struct sdis_medium* solid_back = NULL;
- struct sdis_medium* mdm;
- double lambda_front, lambda_back;
- double delta_front, delta_back;
- double delta_boundary_front, delta_boundary_back;
- double proba;
- double tmp;
- double r;
- double power;
- double tcr;
- float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
- float dir_front[DIM], dir_back[DIM];
- float* dir;
- float reinject_dst_front = 0, reinject_dst_back = 0;
- float reinject_dst;
- /* In 2D it is useless to try to resample a reinjection direction since there
- * is only one possible direction */
- const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
- int iattempt;
- int move;
- int reinjection_is_valid;
- res_T res = RES_OK;
- ASSERT(scn && ctx && frag && rwalk && rng && T);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- (void)frag, (void)ctx;
-
- /* Retrieve the current boundary media */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- solid_front = interface_get_medium(interf, SDIS_FRONT);
- solid_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(solid_front->type == SDIS_SOLID);
- ASSERT(solid_back->type == SDIS_SOLID);
-
- /* Retrieve the thermal contact resistance */
- tcr = interface_get_thermal_contact_resistance(interf, frag);
-
- /* Fetch the properties of the media */
- lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx);
- lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx);
-
- /* Note that reinjection distance is *FIXED*. It MUST ensure that the orthogonal
- * distance from the boundary to the point to challenge is equal to delta. */
- delta_front = solid_get_delta(solid_front, &rwalk->vtx);
- delta_back = solid_get_delta(solid_back, &rwalk->vtx);
- delta_boundary_front = delta_front*sqrt(DIM);
- delta_boundary_back = delta_back *sqrt(DIM);
-
- rwalk_saved = *rwalk;
- reinjection_is_valid = 0;
- iattempt = 0;
- do {
- if(iattempt != 0) *rwalk = rwalk_saved;
-
- /* Sample a reinjection direction and reflect it around the normal. Then
- * reflect them on the back side of the interface. */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
- XD(reflect)(dir2, dir0, rwalk->hit.normal);
- fX(minus)(dir1, dir0);
- fX(minus)(dir3, dir2);
-
- /* Select the reinjection direction and distance for the front side */
- res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front, rwalk,
- dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front, 1, &move,
- &reinjection_is_valid, &hit0);
- if(res != RES_OK) goto error;
- if(!reinjection_is_valid) continue;
-
- /* Select the reinjection direction and distance for the back side */
- res = XD(select_reinjection_dir_and_check_validity)(scn, solid_back, rwalk,
- dir1, dir3, delta_boundary_back, dir_back, &reinject_dst_back, 1, &move,
- &reinjection_is_valid, &hit1);
- if(res != RES_OK) goto error;
- if(!reinjection_is_valid) continue;
-
- /* If random walk was moved by the select_reinjection_dir on back side, one
- * has to rerun the select_reinjection_dir on front side at the new pos */
- if(move) {
- res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front,
- rwalk, dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front,
- 0, NULL, &reinjection_is_valid, &hit0);
- if(res != RES_OK) goto error;
- if(!reinjection_is_valid) continue;
- }
- } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
-
- /* Could not find a valid reinjection */
- if(iattempt >= MAX_ATTEMPTS) {
- *rwalk = rwalk_saved;
- log_warn(scn->dev,
- "%s: could not find a valid solid/solid reinjection at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
-
- r = ssp_rng_canonical(rng);
- if(tcr == 0) { /* No thermal contact resistance */
- /* Define the reinjection side. Note that the proba should be : Lf/Df' /
- * (Lf/Df' + Lb/Db')
- *
- * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted
- * delta of the <front|back> side, i.e. : D<f|b>' =
- * reinject_dst_<front|back> / sqrt(DIM)
- *
- * Anyway, one can avoid to compute the adjusted delta by directly using the
- * adjusted reinjection distance since the resulting proba is strictly the
- * same; sqrt(DIM) can be simplified. */
- proba = (lambda_front/reinject_dst_front)
- / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
- } else {
- const double df = reinject_dst_front/sqrt(DIM);
- const double db = reinject_dst_back/sqrt(DIM);
- const double tmp_front = lambda_front/df;
- const double tmp_back = lambda_back/db;
- const double tmp_r = tcr*tmp_front*tmp_back;
- switch(rwalk->hit_side) {
- case SDIS_BACK:
- /* When coming from the BACK side, the probability to be reinjected on
- * the FRONT side depends on the thermal contact resistance: it
- * decreases when the TCR increases (and tends to 0 when TCR -> +inf) */
- proba = (tmp_front) / (tmp_front + tmp_back + tmp_r);
- break;
- case SDIS_FRONT:
- /* Same thing when coming from the FRONT side: the probability of
- * reinjection on the FRONT side depends on the thermal contact
- * resistance: it increases when the TCR increases (and tends to 1 when
- * the TCR -> +inf) */
- proba = (tmp_front + tmp_r) / (tmp_front + tmp_back + tmp_r);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- }
-
- if(r < proba) { /* Reinject in front */
- dir = dir_front;
- hit = &hit0;
- mdm = solid_front;
- reinject_dst = reinject_dst_front;
- } else { /* Reinject in back */
- dir = dir_back;
- hit = &hit1;
- mdm = solid_back;
- reinject_dst = reinject_dst_back;
- }
-
- /* Handle the volumic power */
- power = solid_get_volumic_power(mdm, &rwalk->vtx);
- if(power != SDIS_VOLUMIC_POWER_NONE) {
- const double delta_in_meter = reinject_dst * scn->fp_to_meter;
- const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += power * tmp;
-
- if(ctx->green_path) {
- res = green_path_add_power_term(ctx->green_path, mdm, &rwalk->vtx, tmp);
- if(res != RES_OK) goto error;
- }
- }
-
- /* Time rewind */
- res = XD(time_rewind)(mdm, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T);
- if(res != RES_OK) goto error;
- if(T->done) goto exit; /* Limit condition was reached */
-
- /* Perform reinjection. */
- XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
- if(hit->distance == reinject_dst) {
- T->func = XD(boundary_path);
- rwalk->mdm = NULL;
- rwalk->hit = *hit;
- rwalk->hit_side = fX(dot)(hit->normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
- } else {
- T->func = XD(conductive_path);
- rwalk->mdm = mdm;
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
- }
-
- /* Register the new vertex against the heat path */
- res = register_heat_vertex
- (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
- if(res != RES_OK) goto error;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-XD(solid_fluid_boundary_path)
- (const struct sdis_scene* scn,
- const struct rwalk_context* ctx,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm_front = NULL;
- struct sdis_medium* mdm_back = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_medium* fluid = NULL;
- struct XD(rwalk) rwalk_saved;
- struct sXd(hit) hit = SXD_HIT_NULL;
- struct sdis_interface_fragment frag_fluid;
- double hc;
- double hr;
- double epsilon; /* Interface emissivity */
- double lambda;
- double fluid_proba;
- double radia_proba;
- double delta;
- double delta_boundary;
- double r;
- double tmp;
- float dir0[DIM], dir1[DIM];
- float reinject_dst;
- /* In 2D it is useless to try to resample a reinjection direction since there
- * is only one possible direction */
- const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
- int iattempt;
- int reinjection_is_valid = 0;
- res_T res = RES_OK;
- ASSERT(scn && rwalk && rng && T && ctx);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
-
- /* Retrieve the solid and the fluid split by the boundary */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm_front = interface_get_medium(interf, SDIS_FRONT);
- mdm_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(mdm_front->type != mdm_back->type);
-
- frag_fluid = *frag;
- if(mdm_front->type == SDIS_SOLID) {
- solid = mdm_front;
- fluid = mdm_back;
- frag_fluid.side = SDIS_BACK;
- } else {
- solid = mdm_back;
- fluid = mdm_front;
- frag_fluid.side = SDIS_FRONT;
- }
-
- /* Fetch the solid properties */
- lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta = solid_get_delta(solid, &rwalk->vtx);
-
- /* Note that the reinjection distance is *FIXED*. It MUST ensure that the
- * orthogonal distance from the boundary to the point to chalenge is equal to
- * delta. */
- delta_boundary = sqrt(DIM) * delta;
-
- rwalk_saved = *rwalk;
- reinjection_is_valid = 0;
- iattempt = 0;
- do {
- if(iattempt != 0) *rwalk = rwalk_saved;
-
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
-
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
-
- if(solid == mdm_back) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
-
- /* Select the solid reinjection direction and distance */
- res = XD(select_reinjection_dir_and_check_validity)(scn, solid, rwalk,
- dir0, dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
- &reinjection_is_valid, &hit);
- if(res != RES_OK) goto error;
-
- } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
-
- /* Could not find a valid reinjecton */
- if(iattempt >= MAX_ATTEMPTS) {
- *rwalk = rwalk_saved;
- log_warn(scn->dev,
- "%s: could not find a valid solid/fluid reinjection at {%g, %g %g}.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
-
- /* Define the orthogonal dst from the reinjection pos to the interface */
- delta = reinject_dst / sqrt(DIM);
-
- /* Fetch the boundary properties */
- epsilon = interface_side_get_emissivity(interf, &frag_fluid);
- hc = interface_get_convection_coef(interf, frag);
-
- /* Compute the radiative coefficient */
- hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
-
- /* Compute the probas to switch in solid, fluid or radiative random walk */
- tmp = lambda / (delta * scn->fp_to_meter);
- fluid_proba = hc / (tmp + hr + hc);
- radia_proba = hr / (tmp + hr + hc);
- /*solid_proba = tmp / (tmp + hr + hc);*/
-
- r = ssp_rng_canonical(rng);
- if(r < radia_proba) { /* Switch in radiative random walk */
- T->func = XD(radiative_path);
- rwalk->mdm = fluid;
- rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
- } else if(r < fluid_proba + radia_proba) { /* Switch to convective random walk */
- T->func = XD(convective_path);
- rwalk->mdm = fluid;
- rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
- } else { /* Solid random walk */
- /* Handle the volumic power */
- const double power = solid_get_volumic_power(solid, &rwalk->vtx);
- if(power != SDIS_VOLUMIC_POWER_NONE) {
- const double delta_in_meter = reinject_dst * scn->fp_to_meter;
- tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += power * tmp;
-
- if(ctx->green_path) {
- res = green_path_add_power_term(ctx->green_path, solid, &rwalk->vtx, tmp);
- if(res != RES_OK) goto error;
- }
- }
-
- /* Time rewind */
- res = XD(time_rewind)(solid, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T);
- if(res != RES_OK) goto error;
- if(T->done) goto exit; /* Limit condition was reached */
-
- /* Perform solid reinjection */
- XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
- if(hit.distance == reinject_dst) {
- T->func = XD(boundary_path);
- rwalk->mdm = NULL;
- rwalk->hit = hit;
- rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
- } else {
- T->func = XD(conductive_path);
- rwalk->mdm = solid;
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
- }
-
- /* Register the new vertex against the heat path */
- res = register_heat_vertex
- (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-XD(solid_boundary_with_flux_path)
- (const struct sdis_scene* scn,
- const struct rwalk_context* ctx,
- const struct sdis_interface_fragment* frag,
- const double phi,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- struct XD(rwalk) rwalk_saved;
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm = NULL;
- double lambda;
- double delta;
- double delta_boundary;
- double delta_in_meter;
- double power;
- double tmp;
- struct sXd(hit) hit;
- float dir0[DIM];
- float dir1[DIM];
- float reinject_dst;
- /* In 2D it is useless to try to resample a reinjection direction since there
- * is only one possible direction */
- const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
- int iattempt = 0;
- int reinjection_is_valid = 0;
- res_T res = RES_OK;
- ASSERT(frag && phi != SDIS_FLUX_NONE);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- (void)ctx;
-
- /* Fetch current interface */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- ASSERT(phi == interface_side_get_flux(interf, frag));
-
- /* Fetch incoming solid */
- mdm = interface_get_medium(interf, frag->side);
- ASSERT(mdm->type == SDIS_SOLID);
-
- /* Fetch medium properties */
- lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- delta = solid_get_delta(mdm, &rwalk->vtx);
-
- /* Compute the reinjection distance. It MUST ensure that the orthogonal
- * distance from the boundary to the point to chalenge is equal to delta. */
- delta_boundary = delta * sqrt(DIM);
-
- rwalk_saved = *rwalk;
- reinjection_is_valid = 0;
- iattempt = 0;
- do {
- if(iattempt != 0) *rwalk = rwalk_saved;
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
-
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
-
- if(frag->side == SDIS_BACK) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
-
- /* Select the reinjection direction and distance */
- res = XD(select_reinjection_dir_and_check_validity)(scn, mdm, rwalk, dir0,
- dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
- &reinjection_is_valid, &hit);
- if(res != RES_OK) goto error;
-
- } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
-
- /* Could not find a valid reinjecton */
- if(iattempt >= MAX_ATTEMPTS) {
- *rwalk = rwalk_saved;
- log_warn(scn->dev,
- "%s: could not find a valid solid/fluid with flux reinjection "
- "at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
-
- /* Define the orthogonal dst from the reinjection pos to the interface */
- delta = reinject_dst / sqrt(DIM);
-
- /* Handle the flux */
- delta_in_meter = delta * scn->fp_to_meter;
- tmp = delta_in_meter / lambda;
- T->value += phi * tmp;
- if(ctx->green_path) {
- res = green_path_add_flux_term(ctx->green_path, interf, frag, tmp);
- if(res != RES_OK) goto error;
- }
-
- /* Handle the volumic power */
- power = solid_get_volumic_power(mdm, &rwalk->vtx);
- if(power != SDIS_VOLUMIC_POWER_NONE) {
- delta_in_meter = reinject_dst * scn->fp_to_meter;
- tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += power * tmp;
- if(ctx->green_path) {
- res = green_path_add_power_term(ctx->green_path, mdm, &rwalk->vtx, tmp);
- if(res != RES_OK) goto error;
- }
- }
-
- /* Time rewind */
- res = XD(time_rewind)(mdm, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T);
- if(res != RES_OK) goto error;
- if(T->done) goto exit; /* Limit condition was reached */
-
- /* Reinject. If the reinjection move the point too close of a boundary,
- * assume that the zone is isotherm and move to the boundary. */
- XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
- if(hit.distance == reinject_dst) {
- T->func = XD(boundary_path);
- rwalk->mdm = NULL;
- rwalk->hit = hit;
- rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
- } else {
- T->func = XD(conductive_path);
- rwalk->mdm = mdm;
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
- }
-
- /* Register the new vertex against the heat path */
- res = register_heat_vertex
- (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
- if(res != RES_OK) goto error;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -1031,7 +90,7 @@ XD(boundary_path)
if(mdm_front->type == mdm_back->type) {
res = XD(solid_solid_boundary_path)(scn, ctx, &frag, rwalk, rng, T);
} else {
- res = XD(solid_fluid_boundary_path)(scn, ctx, &frag, rwalk, rng, T);
+ res = XD(solid_fluid_boundary_picard1_path)(scn, ctx, &frag, rwalk, rng, T);
}
if(res != RES_OK) goto error;
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -0,0 +1,828 @@
+/* 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_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_scene_c.h"
+
+#include <star/ssp.h>
+
+#include "sdis_Xd_begin.h"
+
+struct XD(find_reinjection_ray_args) {
+ const struct sdis_medium* solid; /* Medium into which the reinjection occurs */
+ const struct XD(rwalk)* rwalk; /* Current random walk state */
+ float dir0[DIM]; /* Challenged ray direction */
+ float dir1[DIM]; /* Challenged ray direction */
+ double distance; /* Maximum reinjection distance */
+
+ /* Define if the random walk position can be moved or not to find a valid
+ * reinjection direction */
+ int can_move;
+};
+static const struct XD(find_reinjection_ray_args)
+XD(FIND_REINJECTION_RAY_ARGS_NULL) = { NULL, NULL, {0}, {0}, 0, 0 };
+
+struct XD(reinjection_ray) {
+ double org[DIM]; /* Origin of the reinjection */
+ float dir[DIM]; /* Direction of the reinjection */
+ float dst; /* Reinjection distance along dir */
+ struct sXd(hit) hit; /* Hit along the reinjection dir */
+
+ /* Define whether or not the random walk was moved to find this reinjection
+ * ray */
+ int position_was_moved;
+};
+static const struct XD(reinjection_ray)
+XD(REINJECTION_RAY_NULL) = { {0}, {0}, 0, SXD_HIT_NULL__, 0 };
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+XD(check_find_reinjection_ray_args)
+ (const struct XD(find_reinjection_ray_args)* args)
+{
+ return args
+ && args->solid
+ && args->rwalk
+ && args->distance > 0
+ && fX(is_normalized)(args->dir0)
+ && fX(is_normalized)(args->dir1);
+}
+
+static INLINE int
+XD(check_sample_reinjection_step_args)
+ (const struct XD(sample_reinjection_step_args)* args)
+{
+ return args
+ && args->rng
+ && args->solid
+ && args->solid->type == SDIS_SOLID
+ && args->rwalk
+ && args->distance > 0
+ && (unsigned)args->side < SDIS_SIDE_NULL__;
+}
+
+static INLINE int
+XD(check_reinjection_step)(const struct XD(reinjection_step)* step)
+{
+ return step
+ && fX(is_normalized)(step->direction)
+ && step->distance > 0;
+}
+
+static INLINE int
+XD(check_solid_reinjection_args)(const struct XD(solid_reinjection_args)* args)
+{
+ return args
+ && XD(check_reinjection_step)(args->reinjection)
+ && args->rng
+ && args->rwalk
+ && args->rwalk_ctx
+ && args->T
+ && args->fp_to_meter > 0;
+}
+
+/* Check that the interface fragment is consistent with the current state of
+ * the random walk */
+static INLINE int
+XD(check_rwalk_fragment_consistency)
+ (const struct XD(rwalk)* rwalk,
+ const struct sdis_interface_fragment* frag)
+{
+ double N[DIM];
+ double uv[2] = {0, 0};
+ ASSERT(rwalk && frag);
+ dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
+ if( SXD_HIT_NONE(&rwalk->hit)
+ || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
+ || !dX(eq_eps)(N, frag->Ng, 1.e-6)
+ || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
+ || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
+ return 0;
+ }
+#if (SDIS_XD_DIMENSION == 2)
+ uv[0] = rwalk->hit.u;
+#else
+ d2_set_f2(uv, rwalk->hit.uv);
+#endif
+ return d2_eq_eps(uv, frag->uv, 1.e-6);
+}
+
+static void
+XD(sample_reinjection_dir)
+ (const struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ float dir[DIM])
+{
+#if DIM == 2
+ /* The sampled directions is defined by rotating the normal around the Z axis
+ * of an angle of PI/4 or -PI/4. Let the rotation matrix defined as
+ * | cos(a) -sin(a) |
+ * | sin(a) cos(a) |
+ * with a = PI/4, dir = sqrt(2)/2 * | 1 -1 | . N
+ * | 1 1 |
+ * with a =-PI/4, dir = sqrt(2)/2 * | 1 1 | . N
+ * |-1 1 |
+ * Note that since the sampled direction is finally normalized, we can
+ * discard the sqrt(2)/2 constant. */
+ const uint64_t r = ssp_rng_uniform_uint64(rng, 0, 1);
+ ASSERT(rwalk && rng && dir);
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ ASSERT(!rwalk->mdm);
+
+ if(r) {
+ dir[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
+ dir[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ } else {
+ dir[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ dir[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ }
+ f2_normalize(dir, dir);
+#else
+ /* Sample a random direction around the normal whose cosine is 1/sqrt(3). To
+ * do so we sample a position onto a cone whose height is 1/sqrt(2) and the
+ * radius of its base is 1. */
+ float frame[9];
+ ASSERT(rwalk && rng && dir);
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ ASSERT(!rwalk->mdm);
+ ASSERT(fX(is_normalized)(rwalk->hit.normal));
+
+ ssp_ran_circle_uniform_float(rng, dir, NULL);
+ dir[2] = (float)(1.0/sqrt(2));
+
+ f33_basis(frame, rwalk->hit.normal);
+ f33_mulf3(dir, frame, dir);
+ f3_normalize(dir, dir);
+ ASSERT(eq_epsf(f3_dot(dir, rwalk->hit.normal), (float)(1.0/sqrt(3)), 1.e-4f));
+#endif
+}
+
+
+#if DIM == 2
+static void
+XD(move_away_primitive_boundaries)
+ (const struct XD(rwalk)* rwalk,
+ const double delta,
+ double position[DIM]) /* Position to move */
+{
+ struct sXd(attrib) attr;
+ float pos[DIM];
+ float dir[DIM];
+ float len;
+ const float st = 0.5f;
+ ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->hit) && delta > 0);
+
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr));
+
+ fX_set_dX(pos, position);
+ fX(sub)(dir, attr.value, pos);
+ len = fX(normalize)(dir, dir);
+ len = MMIN(len, (float)(delta*0.1));
+
+ XD(move_pos)(position, dir, len);
+}
+#else
+/* Move the submitted position away from the primitive boundaries to avoid
+ * numerical issues leading to inconsistent random walks. */
+static void
+XD(move_away_primitive_boundaries)
+ (const struct XD(rwalk)* rwalk,
+ const double delta,
+ double position[DIM])
+{
+ struct s3d_attrib v0, v1, v2; /* Triangle vertices */
+ float E[3][4]; /* 3D edge equations */
+ float dst[3]; /* Distance from current position to edge equation */
+ float N[3]; /* Triangle normal */
+ float P[3]; /* Random walk position */
+ float tmp[3];
+ float min_dst, max_dst;
+ float cos_a1, cos_a2;
+ float len;
+ int imax = 0;
+ int imin = 0;
+ int imid = 0;
+ int i;
+ ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->hit));
+
+ fX_set_dX(P, position);
+
+ /* Fetch triangle vertices */
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the edge vector */
+ f3_sub(E[0], v1.value, v0.value);
+ f3_sub(E[1], v2.value, v1.value);
+ f3_sub(E[2], v0.value, v2.value);
+
+ /* Compute the triangle normal */
+ f3_cross(N, E[1], E[0]);
+
+ /* Compute the 3D edge equation */
+ f3_normalize(E[0], f3_cross(E[0], E[0], N));
+ f3_normalize(E[1], f3_cross(E[1], E[1], N));
+ f3_normalize(E[2], f3_cross(E[2], E[2], N));
+ E[0][3] = -f3_dot(E[0], v0.value);
+ E[1][3] = -f3_dot(E[1], v1.value);
+ E[2][3] = -f3_dot(E[2], v2.value);
+
+ /* Compute the distance from current position to the edges */
+ dst[0] = f3_dot(E[0], P) + E[0][3];
+ dst[1] = f3_dot(E[1], P) + E[1][3];
+ dst[2] = f3_dot(E[2], P) + E[2][3];
+
+ /* Retrieve the min and max distance from random walk position to triangle
+ * edges */
+ min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
+ max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
+
+ /* Sort the edges with respect to their distance to the random walk position */
+ FOR_EACH(i, 0, 3) {
+ if(dst[i] == min_dst) {
+ imin = i;
+ } else if(dst[i] == max_dst) {
+ imax = i;
+ } else {
+ imid = i;
+ }
+ }
+ (void)imax;
+
+ /* TODO if the current position is near a vertex, one should move toward the
+ * farthest edge along its normal to avoid too small displacement */
+
+ /* Compute the distance `dst' from the current position to the edges to move
+ * to, along the normal of the edge from which the random walk is the nearest
+ *
+ * +. cos(a) = d / dst => dst = d / cos_a
+ * / `*.
+ * / | `*.
+ * / dst| a /`*.
+ * / | / `*.
+ * / | / d `*.
+ * / |/ `*.
+ * +---------o----------------+ */
+ cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
+ cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
+ dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
+ dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
+
+ /* Compute the maximum displacement distance into the triangle along the
+ * normal of the edge from which the random walk is the nearest */
+ len = MMIN(dst[imid], dst[imax]);
+ ASSERT(len != FLT_MAX);
+
+ /* Define the displacement distance as the minimum between 10 percent of
+ * delta and len / 2. */
+ len = MMIN(len*0.5f, (float)(delta*0.1));
+ XD(move_pos)(position, E[imin], len);
+}
+#endif
+
+static res_T
+XD(find_reinjection_ray)
+ (const struct sdis_scene* scn,
+ const struct XD(find_reinjection_ray_args)* args,
+ struct XD(reinjection_ray)* ray)
+{
+ /* Emperical scale factor applied to the challenged reinjection distance. If
+ * the distance to reinject is less than this adjusted value, the solver will
+ * try to discard the reinjection distance if possible in order to avoid
+ * numerical issues. */
+ const float REINJECT_DST_MIN_SCALE = 0.125f;
+
+ /* # attempts to find a ray direction */
+ int MAX_ATTEMPTS = 1;
+
+ /* Physical properties */
+ struct sdis_interface* interf;
+ struct sdis_medium* mdm0;
+ struct sdis_medium* mdm1;
+
+ struct hit_filter_data filter_data;
+ struct sXd(hit) hit;
+ struct sXd(hit) hit0;
+ struct sXd(hit) hit1;
+ double tmp[DIM];
+ double dst;
+ double dst0;
+ double dst1;
+ const float* dir;
+ float reinject_threshold;
+ double dst_adjusted;
+ float org[DIM];
+ const float range[2] = {0, FLT_MAX};
+ enum sdis_side side;
+ int iattempt = 0;
+ res_T res = RES_OK;
+
+ ASSERT(scn && args && ray);
+ ASSERT(XD(check_find_reinjection_ray_args)(args));
+
+ *ray = XD(REINJECTION_RAY_NULL);
+ MAX_ATTEMPTS = args->can_move ? 2 : 1;
+
+ dst_adjusted = args->distance * RAY_RANGE_MAX_SCALE;
+ reinject_threshold = (float)args->distance * REINJECT_DST_MIN_SCALE;
+
+ dX(set)(ray->org, args->rwalk->vtx.P);
+
+ do {
+ fX_set_dX(org, ray->org);
+ filter_data.XD(hit) = args->rwalk->hit;
+ filter_data.epsilon = args->distance * 0.01;
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), org, args->dir0, range, &filter_data, &hit0));
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), org, args->dir1, range, &filter_data, &hit1));
+
+ /* Retrieve the medium at the reinjection pos along dir0 */
+ if(SXD_HIT_NONE(&hit0)) {
+ XD(move_pos)(dX(set)(tmp, ray->org), args->dir0, (float)args->distance);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0);
+ if(res == RES_BAD_OP) { mdm0 = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit0.prim.prim_id);
+ side = fX(dot)(args->dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm0 = interface_get_medium(interf, side);
+ }
+
+ /* Retrieve the medium at the reinjection pos along dir1 */
+ if(SXD_HIT_NONE(&hit1)) {
+ XD(move_pos)(dX(set)(tmp, ray->org), args->dir1, (float)args->distance);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1);
+ if(res == RES_BAD_OP) { mdm1 = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit1.prim.prim_id);
+ side = fX(dot)(args->dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm1 = interface_get_medium(interf, side);
+ }
+
+ dst0 = dst1 = -1;
+ if(mdm0 == args->solid) { /* Check reinjection consistency */
+ if(hit0.distance <= dst_adjusted) {
+ dst0 = hit0.distance;
+ } else {
+ dst0 = args->distance;
+ hit0 = SXD_HIT_NULL;
+ }
+ }
+ if(mdm1 == args->solid) { /* Check reinjection consistency */
+ if(hit1.distance <= dst_adjusted) {
+ dst1 = hit1.distance;
+ } else {
+ dst1 = args->distance;
+ hit1 = SXD_HIT_NULL;
+ }
+ }
+
+ /* No valid reinjection. Maybe the random walk is near a sharp corner and
+ * thus the ray-tracing misses the enclosure geometry. Another possibility
+ * is that the random walk lies roughly on an edge. In this case, sampled
+ * reinjection dirs can intersect the primitive on the other side of the
+ * edge. Normally, this primitive should be filtered by the "hit_filter"
+ * function but this may be not the case due to a "threshold effect". In
+ * both situations, try to slightly move away from the primitive boundaries
+ * and retry to find a valid reinjection. */
+ if(dst0 == -1 && dst1 == -1) {
+ XD(move_away_primitive_boundaries)(args->rwalk, args->distance, ray->org);
+ ray->position_was_moved = 1;
+ }
+ } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS);
+
+ if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */
+#if DIM == 2
+ log_warn(scn->dev, "%s: no valid reinjection direction at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(ray->org));
+#else
+ log_warn(scn->dev, "%s: no valid reinjection direction at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(ray->org));
+#endif
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ if(dst0 == -1) {
+ /* Invalid dir0 -> move along dir1 */
+ dir = args->dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 == -1) {
+ /* Invalid dir1 -> move along dir0 */
+ dir = args->dir0;
+ dst = dst0;
+ hit = hit0;
+ } else if(dst0 < reinject_threshold && dst1 < reinject_threshold) {
+ /* The displacement along dir0 and dir1 are both below the reinjection
+ * threshold that defines a distance under which the temperature gradients
+ * are ignored. Move along the direction that allows the maximum
+ * displacement. */
+ if(dst0 > dst1) {
+ dir = args->dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dir = args->dir1;
+ dst = dst1;
+ hit = hit1;
+ }
+ } else if(dst0 < reinject_threshold) {
+ /* Ingore dir0 that is bellow the reinject threshold */
+ dir = args->dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 < reinject_threshold) {
+ /* Ingore dir1 that is bellow the reinject threshold */
+ dir = args->dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ /* All reinjection directions are valid. Choose the first 1 that was
+ * randomly selected by the sample_reinjection_dir procedure and adjust
+ * the displacement distance. */
+ dir = args->dir0;
+
+ /* Define the reinjection distance along dir0 and its corresponding hit */
+ if(dst0 <= dst1) {
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dst = dst1;
+ hit = SXD_HIT_NULL;
+ }
+
+ /* If the displacement distance is too close of a boundary, move to the
+ * boundary in order to avoid numerical uncertainty. */
+ if(!SXD_HIT_NONE(&hit0)
+ && dst0 != dst
+ && eq_eps(dst0, dst, dst0*0.1)) {
+ dst = dst0;
+ hit = hit0;
+ }
+ }
+
+ /* Setup the ray */
+ fX(set)(ray->dir, dir);
+ ray->dst = (float)dst;
+ ray->hit = hit;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(find_reinjection_ray_and_check_validity)
+ (const struct sdis_scene* scn,
+ const struct XD(find_reinjection_ray_args)* args,
+ struct XD(reinjection_ray)* ray)
+{
+ double pos[DIM];
+ struct sdis_medium* reinject_mdm;
+ res_T res = RES_OK;
+
+ ASSERT(scn && args && ray);
+ ASSERT(XD(check_find_reinjection_ray_args)(args));
+
+ /* Select a reinjection direction */
+ res = XD(find_reinjection_ray)(scn, args, ray);
+ if(res != RES_OK) goto error;
+
+ if(SXD_HIT_NONE(&ray->hit)) {
+ /* Check medium consistency at the reinjection position */
+ XD(move_pos)(dX(set)(pos, ray->org), ray->dir, (float)ray->dst);
+ res = scene_get_medium_in_closed_boundaries(scn, pos, &reinject_mdm);
+ if(res != RES_OK) goto error;
+
+ if(reinject_mdm != args->solid) {
+ res = RES_BAD_OP;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+XD(sample_reinjection_step_solid_fluid)
+ (const struct sdis_scene* scn,
+ const struct XD(sample_reinjection_step_args)* args,
+ struct XD(reinjection_step)* step)
+{
+ /* Input/output data of the function finding a valid reinjection ray */
+ struct XD(find_reinjection_ray_args) find_reinject_ray_args =
+ XD(FIND_REINJECTION_RAY_ARGS_NULL);
+ struct XD(reinjection_ray) ray = XD(REINJECTION_RAY_NULL);
+
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+
+ /* Miscellaneous variables */
+ float dir0[DIM]; /* Sampled direction */
+ float dir1[DIM]; /* Sampled direction reflected */
+ int iattempt = 0; /* #attempts to find a reinjection dir */
+ res_T res = RES_OK;
+
+ /* Pre-conditions */
+ ASSERT(scn && args && step);
+ ASSERT(XD(check_sample_reinjection_step_args)(args));
+
+ iattempt = 0;
+ do {
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(args->rwalk, args->rng, dir0);
+
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, args->rwalk->hit.normal);
+
+ /* Flip the sampled directions if one wants to reinject to back side */
+ if(args->side == SDIS_BACK) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
+ }
+
+ /* Find the reinjection step */
+ find_reinject_ray_args.solid = args->solid;
+ find_reinject_ray_args.rwalk = args->rwalk;
+ find_reinject_ray_args.distance = args->distance;
+ find_reinject_ray_args.can_move = 1;
+ fX(set)(find_reinject_ray_args.dir0, dir0);
+ fX(set)(find_reinject_ray_args.dir1, dir1);
+ res = XD(find_reinjection_ray_and_check_validity)
+ (scn, &find_reinject_ray_args, &ray);
+ if(res == RES_BAD_OP) continue; /* Cannot find a valid reinjection ray. Retry */
+ if(res != RES_OK) goto error;
+
+ } while(res != RES_OK && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjecton step */
+ if(iattempt >= MAX_ATTEMPTS) {
+ log_warn(scn->dev,
+ "%s: could not find a valid reinjection step at `%g %g %g'.\n",
+ FUNC_NAME, SPLIT3(args->rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ /* Setup the reinjection step */
+ step->hit = ray.hit;
+ step->distance = ray.dst;
+ fX(set)(step->direction, ray.dir);
+
+ /* Update the random walk position if necessary */
+ if(ray.position_was_moved) {
+ dX(set)(args->rwalk->vtx.P, ray.org);
+ }
+
+ /* Post-conditions */
+ ASSERT(dX(eq)(args->rwalk->vtx.P, ray.org));
+ ASSERT(XD(check_reinjection_step)(step));
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+XD(sample_reinjection_step_solid_solid)
+ (const struct sdis_scene* scn,
+ const struct XD(sample_reinjection_step_args)* args_frt,
+ const struct XD(sample_reinjection_step_args)* args_bck,
+ struct XD(reinjection_step)* step_frt,
+ struct XD(reinjection_step)* step_bck)
+{
+ /* Input/output data of the function finding a valid reinjection ray */
+ struct XD(find_reinjection_ray_args) find_reinject_ray_frt_args =
+ XD(FIND_REINJECTION_RAY_ARGS_NULL);
+ struct XD(find_reinjection_ray_args) find_reinject_ray_bck_args =
+ XD(FIND_REINJECTION_RAY_ARGS_NULL);
+ struct XD(reinjection_ray) ray_frt = XD(REINJECTION_RAY_NULL);
+ struct XD(reinjection_ray) ray_bck = XD(REINJECTION_RAY_NULL);
+
+ /* Initial random walk position used as a backup */
+ double rwalk_pos_backup[DIM];
+
+ /* Variables shared by the two side */
+ struct XD(rwalk)* rwalk = NULL;
+ struct ssp_rng* rng = NULL;
+
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+
+ float dir_frt_samp[DIM]; /* Sampled direction */
+ float dir_frt_refl[DIM]; /* Sampled direction reflected */
+ float dir_bck_samp[DIM]; /* Negated sampled direction */
+ float dir_bck_refl[DIM]; /* Negated sampled direction reflected */
+ int iattempt = 0; /* #attempts to find a reinjection dir */
+ res_T res = RES_OK;
+
+ /* Pre-conditions */
+ ASSERT(scn && args_frt && args_bck && step_frt && step_bck);
+ ASSERT(XD(check_sample_reinjection_step_args)(args_frt));
+ ASSERT(XD(check_sample_reinjection_step_args)(args_bck));
+ ASSERT(args_frt->side == SDIS_FRONT);
+ ASSERT(args_bck->side == SDIS_BACK);
+
+ rng = args_frt->rng;
+ rwalk = args_frt->rwalk;
+ ASSERT(args_bck->rng == rng);
+ ASSERT(args_bck->rwalk == rwalk);
+
+ dX(set)(rwalk_pos_backup, rwalk->vtx.P);
+ iattempt = 0;
+ do {
+ /* Restore random walk pos */
+ if(iattempt != 0) dX(set)(rwalk->vtx.P, rwalk_pos_backup);
+
+ /* Sample a reinjection direction and reflect it around the normal. Then
+ * reflect them on the back side of the interface. */
+ XD(sample_reinjection_dir)(rwalk, rng, dir_frt_samp);
+ XD(reflect)(dir_frt_refl, dir_frt_samp, rwalk->hit.normal);
+ fX(minus)(dir_bck_samp, dir_frt_samp);
+ fX(minus)(dir_bck_refl, dir_frt_refl);
+
+ /* Find the reinjection ray for the front side */
+ find_reinject_ray_frt_args.solid = args_frt->solid;
+ find_reinject_ray_frt_args.rwalk = args_frt->rwalk;
+ find_reinject_ray_frt_args.distance = args_frt->distance;
+ find_reinject_ray_frt_args.can_move = 1;
+ fX(set)(find_reinject_ray_frt_args.dir0, dir_frt_samp);
+ fX(set)(find_reinject_ray_frt_args.dir1, dir_frt_refl);
+ res = XD(find_reinjection_ray_and_check_validity)
+ (scn, &find_reinject_ray_frt_args, &ray_frt);
+ if(res == RES_BAD_OP) continue;
+ if(res != RES_OK) goto error;
+
+ /* Update the random walk position if necessary */
+ if(ray_frt.position_was_moved) dX(set)(rwalk->vtx.P, ray_frt.org);
+
+ /* Select the reinjection direction and distance for the back side */
+ find_reinject_ray_bck_args.solid = args_bck->solid;
+ find_reinject_ray_bck_args.rwalk = args_bck->rwalk;
+ find_reinject_ray_bck_args.distance = args_bck->distance;
+ find_reinject_ray_bck_args.can_move = 1;
+ fX(set)(find_reinject_ray_bck_args.dir0, dir_bck_samp);
+ fX(set)(find_reinject_ray_bck_args.dir1, dir_bck_refl);
+ res = XD(find_reinjection_ray_and_check_validity)
+ (scn, &find_reinject_ray_bck_args, &ray_bck);
+ if(res == RES_BAD_OP) continue;
+ if(res != RES_OK) goto error;
+
+ /* Update the random walk position if necessary */
+ if(ray_bck.position_was_moved) dX(set)(rwalk->vtx.P, ray_bck.org);
+
+ /* If random walk was moved to find a valid rinjection ray on back side,
+ * one has to find a valid reinjection ob front side from the new pos */
+ if(ray_bck.position_was_moved) {
+ find_reinject_ray_frt_args.can_move = 0;
+ res = XD(find_reinjection_ray_and_check_validity)
+ (scn, &find_reinject_ray_frt_args, &ray_frt);
+ if(res == RES_BAD_OP) continue;
+ if(res != RES_OK) goto error;
+
+ /* Update the random walk position if necessary */
+ if(ray_frt.position_was_moved) dX(set)(rwalk->vtx.P, ray_frt.org);
+ }
+ } while(res != RES_OK && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjection */
+ if(iattempt >= MAX_ATTEMPTS) {
+ dX(set)(rwalk->vtx.P, rwalk_pos_backup); /* Restore random walk pos */
+ log_warn(scn->dev,
+ "%s: could not find a valid solid/solid reinjection at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ /* Setup the front and back reinjection steps */
+ step_frt->hit = ray_frt.hit;
+ step_bck->hit = ray_bck.hit;
+ step_frt->distance = ray_frt.dst;
+ step_bck->distance = ray_bck.dst;
+ fX(set)(step_frt->direction, ray_frt.dir);
+ fX(set)(step_bck->direction, ray_bck.dir);
+
+ /* Post-conditions */
+ ASSERT(XD(check_reinjection_step)(step_frt));
+ ASSERT(XD(check_reinjection_step)(step_bck));
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+XD(solid_reinjection)
+ (struct sdis_medium* solid,
+ struct XD(solid_reinjection_args)* args)
+{
+ double power;
+ double lambda;
+ double reinject_dst_m; /* Reinjection distance in meters */
+ res_T res = RES_OK;
+ ASSERT(solid && XD(check_solid_reinjection_args)(args));
+
+ reinject_dst_m = args->reinjection->distance * args->fp_to_meter;
+
+ /* Fetch solid properties */
+ lambda = solid_get_thermal_conductivity(solid, &args->rwalk->vtx);
+ power = solid_get_volumic_power(solid, &args->rwalk->vtx);
+
+ /* Handle the volumic power */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double reinject_dst_m_sqr = reinject_dst_m * reinject_dst_m;
+ const double power_term = reinject_dst_m_sqr / (2.0 * DIM * lambda);
+
+ args->T->value += power * power_term;
+
+ /* Update the green */
+ if(args->rwalk_ctx->green_path) {
+ res = green_path_add_power_term
+ (args->rwalk_ctx->green_path, solid, &args->rwalk->vtx, power_term);
+ if(res != RES_OK) goto error;
+ }
+ }
+
+ /* Time rewind */
+ res = XD(time_rewind)
+ (solid, args->rng, reinject_dst_m, args->rwalk_ctx, args->rwalk, args->T);
+ if(res != RES_OK) goto error;
+
+ /* Test if a limit condition was reached */
+ if(args->T->done) goto exit;
+
+ /* Move the random walk to the reinjection position */
+ XD(move_pos)
+ (args->rwalk->vtx.P,
+ args->reinjection->direction,
+ args->reinjection->distance);
+
+ /* The random walk is in the solid */
+ if(args->reinjection->hit.distance != args->reinjection->distance) {
+ args->T->func = XD(conductive_path);
+ args->rwalk->mdm = solid;
+ args->rwalk->hit = SXD_HIT_NULL;
+ args->rwalk->hit_side = SDIS_SIDE_NULL__;
+
+ /* The random walk is at a boundary */
+ } else {
+ args->T->func = XD(boundary_path);
+ args->rwalk->mdm = NULL;
+ args->rwalk->hit = args->reinjection->hit;
+ if(fX(dot)(args->reinjection->hit.normal, args->reinjection->direction) < 0) {
+ args->rwalk->hit_side = SDIS_FRONT;
+ } else {
+ args->rwalk->hit_side = SDIS_BACK;
+ }
+ }
+
+ /* Register the new vertex against the heat path */
+ res = register_heat_vertex
+ (args->rwalk_ctx->heat_path,
+ &args->rwalk->vtx,
+ args->T->value,
+ SDIS_HEAT_VERTEX_CONDUCTION);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_fixed_flux.h b/src/sdis_heat_path_boundary_Xd_fixed_flux.h
@@ -0,0 +1,121 @@
+/* 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_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_scene_c.h"
+
+#include <star/ssp.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Boundary path between a solid and a fluid with a fixed flux
+ ******************************************************************************/
+res_T
+XD(solid_boundary_with_flux_path)
+ (struct sdis_scene* scn,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ const double phi,
+ 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);
+
+ /* Reinjection arguments */
+ struct XD(solid_reinjection_args) solid_reinject_args =
+ XD(SOLID_REINJECTION_ARGS_NULL);
+
+ /* Data attached to the boundary */
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+
+ /* Miscellaneous terms */
+ double lambda; /* Solid conductivity */
+ double delta_boundary; /* Orthogonal reinjection dst at the boundary */
+ double delta; /* Orthogonal fitted reinjection dst at the boundary */
+ double delta_m; /* Delta in meters */
+ double flux_term;
+ enum sdis_side solid_side = SDIS_SIDE_NULL__;
+ res_T res = RES_OK;
+
+ ASSERT(frag && phi != SDIS_FLUX_NONE);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ (void)ctx;
+
+ /* Retrieve the solid split by the interface */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ solid = interface_get_medium(interf, frag->side);
+ solid_side = frag->side;
+ ASSERT(phi == interface_side_get_flux(interf, frag));
+ ASSERT(solid->type == SDIS_SOLID);
+
+ /* Fetch the solid properties */
+ lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
+ delta = solid_get_delta(solid, &rwalk->vtx);
+
+ /* 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 = delta * sqrt(DIM);
+
+ /* 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 boundary to the reinjection position */
+ delta = reinject_step.distance / sqrt(DIM);
+ delta_m = delta * scn->fp_to_meter;
+
+ /* Handle the flux */
+ flux_term = delta_m / lambda;
+ T->value += phi * flux_term;
+ if(ctx->green_path) {
+ res = green_path_add_flux_term(ctx->green_path, interf, frag, flux_term);
+ if(res != RES_OK) goto error;
+ }
+
+ /* 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;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -25,176 +25,238 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+XD(check_Tref)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ const double Tref,
+ const char* func_name)
+{
+ ASSERT(scn && pos && func_name);
+
+#if DIM == 2
+ #define STR_VECX "%g %g"
+ #define SPLITX SPLIT2
+#else
+ #define STR_VECX "%g %g %g"
+ #define SPLITX SPLIT3
+#endif
+ if(Tref < 0) {
+ log_err(scn->dev,
+ "%s: invalid reference temperature `%gK' at the position `"STR_VECX"'.\n",
+ func_name, Tref, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+ if(Tref > scn->tmax) {
+ log_err(scn->dev,
+ "%s: invalid maximum temperature `%gK'. The reference temperature `%gK'"
+ "at the position `"STR_VECX"' is greater than this temperature.\n",
+ func_name, scn->tmax, Tref, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+#undef STR_VECX
+#undef SPLITX
+
+ return RES_OK;
+}
+
+static INLINE res_T
+XD(rwalk_get_Tref)
+ (const struct sdis_scene* scn,
+ const struct XD(rwalk)* rwalk,
+ const struct XD(temperature)* T,
+ double* out_Tref)
+{
+ double Tref = -1;
+ res_T res = RES_OK;
+ ASSERT(rwalk && T && out_Tref);
+
+ if(T->done) {
+ /* The path reaches a limit condition, i.e. it goes to the infinity and
+ * fetches the ambient radiative temperature. We do not use the limit
+ * conditions as the reference temperature to make the sampled paths
+ * independant of them. */
+ ASSERT(T->value == scn->trad.temperature);
+ Tref = scn->trad.reference;
+ } else {
+ struct sdis_interface_fragment frag;
+ struct sdis_interface* interf = NULL;
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+
+ /* Fetch the interface where the random walk ends */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ ASSERT(rwalk->hit_side!=SDIS_FRONT || interf->medium_front->type==SDIS_FLUID);
+ ASSERT(rwalk->hit_side!=SDIS_BACK || interf->medium_back->type==SDIS_FLUID);
+
+ /* Fragment on the fluid side of the boundary onto which the rwalk ends */
+ XD(setup_interface_fragment)
+ (&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+
+ Tref = interface_side_get_reference_temperature(interf, &frag);
+ }
+
+ res = XD(check_Tref)(scn, rwalk->vtx.P, Tref, FUNC_NAME);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_Tref = Tref;
+ return res;
+error:
+ Tref = -1;
+ goto exit;
+}
+
+/*******************************************************************************
* Boundary path between a solid and a fluid
******************************************************************************/
res_T
XD(solid_fluid_boundary_picard1_path)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const 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);
+
+ /* Temperature and random walk state of the sampled radiative path */
+ struct XD(temperature) T_s;
+ struct XD(rwalk) rwalk_s;
+
+ /* Fragment on the fluid side of the boundary */
+ struct sdis_interface_fragment frag_fluid;
+
+ /* Data attached to the boundary */
struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm_front = NULL;
- struct sdis_medium* mdm_back = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
- struct XD(rwalk) rwalk_saved;
- struct sXd(hit) hit = SXD_HIT_NULL;
- struct sdis_interface_fragment frag_fluid;
- double hc;
- double hr;
+
+ 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 */
+
double epsilon; /* Interface emissivity */
- double lambda;
- double fluid_proba;
- double radia_proba;
- double delta;
- double delta_boundary;
+ double Tref; /* Reference temperature */
+ double Tref_s; /* Reference temperature of the sampled radiative path */
+ double lambda; /* Solid conductivity */
+ double delta_boundary; /* Orthogonal reinjection dst at the boundary */
+ double delta; /* Orthogonal fitted reinjection dst at the boundary */
+
double r;
- double tmp;
- float dir0[DIM], dir1[DIM];
- float reinject_dst;
- /* In 2D it is useless to try to resample a reinjection direction since there
- * is only one possible direction */
- const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
- int iattempt;
- int reinjection_is_valid = 0;
+ enum sdis_heat_vertex_type current_vertex_type;
+ 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));
/* Retrieve the solid and the fluid split by the boundary */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm_front = interface_get_medium(interf, SDIS_FRONT);
- mdm_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(mdm_front->type != mdm_back->type);
-
- /* Setup the fluid side fragment */
- frag_fluid = *frag;
- if(mdm_front->type == SDIS_SOLID) {
- solid = mdm_front;
- fluid = mdm_back;
- frag_fluid.side = SDIS_BACK;
- } else {
- solid = mdm_back;
- fluid = mdm_front;
- frag_fluid.side = SDIS_FRONT;
+ 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);
}
- /* Fetch the boundary properties */
- epsilon = interface_side_get_emissivity(interf, &frag_fluid);
- Tref = interface_side_get_reference_temperature(interf, &frag_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);
- /* Note that the reinjection distance is *FIXED*. It MUST ensure that the
- * orthogonal distance from the boundary to the point to chalenge is equal to
- * delta. */
- delta_boundary = sqrt(DIM) * delta;
-
- rwalk_saved = *rwalk;
- reinjection_is_valid = 0;
- iattempt = 0;
- do {
- if(iattempt != 0) *rwalk = rwalk_saved;
-
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
-
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
-
- if(solid == mdm_back) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
+ /* Fetch the boundary emissivity */
+ epsilon = interface_side_get_emissivity(interf, &frag_fluid);
- /* Select the solid reinjection direction and distance */
- res = XD(select_reinjection_dir_and_check_validity)(scn, solid, rwalk,
- dir0, dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
- &reinjection_is_valid, &hit);
+ if(epsilon <= 0) {
+ Tref = 0;
+ } else {
+ /* Check the Tref */
+ Tref = interface_side_get_reference_temperature(interf, &frag_fluid);
+ res = XD(check_Tref)(scn, frag_fluid.P, Tref, FUNC_NAME);
if(res != RES_OK) goto error;
+ }
- } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+ /* 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;
- /* Could not find a valid reinjecton */
- if(iattempt >= MAX_ATTEMPTS) {
- *rwalk = rwalk_saved;
- log_warn(scn->dev,
- "%s: could not find a valid solid/fluid reinjection at {%g, %g %g}.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
+ /* 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_dst / sqrt(DIM);
+ 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_rad_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
-
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
+
/* Compute a global upper bound coefficient */
- h_hat = h_conv + h_cond + h_rad_hat;
+ h_hat = h_conv + h_cond + h_radi_hat;
/* Compute the probas to switch in solid, fluid or radiative random walk */
p_conv = h_conv / h_hat;
p_cond = h_cond / h_hat;
+ /* Fetch the current heat vertex type. This will be used below to restart the
+ * registration of the heat path geometry after a null collision */
+ if(ctx->heat_path) {
+ current_vertex_type = heat_path_get_last_vertex(ctx->heat_path)->type;
+ }
+
/* Null collision */
for(;;) {
- r = ssp_rng_canonical(rng);
+ double h_radi; /* Radiative coefficient */
+ double p_radi; /* Radiative proba */
+
+ r = ssp_rng_canonical(rng);
/* Switch in convective path */
if(r < p_conv) {
T->func = XD(convective_path);
rwalk->mdm = fluid;
- rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit_side = fluid_side;
break;
}
/* Switch in conductive path */
if(r < p_conv + p_cond) {
- /* Handle the volumic power */
- const double power = solid_get_volumic_power(solid, &rwalk->vtx);
- if(power != SDIS_VOLUMIC_POWER_NONE) {
- const double delta_in_meter = reinject_dst * scn->fp_to_meter;
- tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += power * tmp;
-
- if(ctx->green_path) {
- res = green_path_add_power_term(ctx->green_path, solid, &rwalk->vtx, tmp);
- if(res != RES_OK) goto error;
- }
- }
-
- /* Time rewind */
- res = XD(time_rewind)(solid, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T);
- if(res != RES_OK) goto error;
- if(T->done) goto exit; /* Limit condition was reached */
-
- /* Perform solid reinjection */
- XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
- if(hit.distance == reinject_dst) {
- T->func = XD(boundary_path);
- rwalk->mdm = NULL;
- rwalk->hit = hit;
- rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
- } else {
- T->func = XD(conductive_path);
- rwalk->mdm = solid;
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
- }
-
- /* Register the new vertex against the heat path */
- res = register_heat_vertex
- (ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
+ 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;
}
@@ -202,47 +264,51 @@ XD(solid_fluid_boundary_picard1_path)
/* From there, we know the path is either a radiative path or a
* null-collision */
- /* Trace a candidate radiative path and get the Tref at its end.
+ /* Sample a radiative path and get the Tref at its end.
* TODO handle the registration of the path geometry */
- T_candidate = *T;
- rwalk_candidate = *rwalk;
- res = XD(radiative_path)(scn, ctx, &rwalk_candidate, rng, T_candidate);
+ 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;
/* Get the Tref at the end of the candidate radiative path */
- if(T_candidate->done) {
- Tref_candidate = T_candidate->value;
- } else {
- ASSERT(!SXD_HIT_NONE(rwalk_candidate->hit));
- XD(setup_interface_fragment)
- (&frag_candidate, &rwalk_candidate->vtx, &rwalk_candidate->hit,
- rwalk_candidate->hit_side);
- interf_candidate = scene_get_interface
- (scn, rwalk_candidate->hit.prim.prim_id);
-
- Tref_candidate = interface_side_get_reference_temperature(interf, f&rag);
- }
+ res = XD(rwalk_get_Tref)(scn, &rwalk_s, &T_s, &Tref_s);
+ if(res != RES_OK) goto error;
- if(Tref_candidate < 0) {
- log_err(scn->dev,
- "%s: invalid reference temperature `%gK' at the position `%g %g %g'.\n",
- FUNC_NAME, Tref_candidate, SPLIT3(rwalk_candidate->vtx.P));
- res = RES_BAD_OP_IRRECOVERABLE;
- goto error;
- }
+ h_radi = BOLTZMANN_CONSTANT * epsilon *
+ ( Tref*Tref*Tref
+ + Tref*Tref * Tref_s
+ + Tref * Tref_s*Tref_s
+ + Tref_s*Tref_s*Tref_s);
- h_rad = BOLTZMANN_CONSTANT
- * epsilon
- * ( Tref*Tref*Tref
- + Tref*Tref * Tref_candidate
- + Tref* Tref_candidate*Tref_candidate
- + Tref_candidate*Tref_candidate*Tref_candidate);
-
- p_rad = h_rad / h_hat;
- if(r < p_conv + p_cond + p_rad) { /* Radiative path */
- *rwalk = *rwalk_candidate;
- *T = *T_candidate;
+ p_radi = h_radi / h_hat;
+ if(r < p_conv + p_cond + p_radi) { /* Radiative path */
+ *rwalk = rwalk_s;
+ *T = T_s;
break;
+
+ /* Null collision: the sampled path is rejected. */
+ } else {
+
+ if(ctx->green_path) {
+ /* The limit condition of the green path could be set by the rejected
+ * sampled radiative path. Reset this limit condition. */
+ green_path_reset_limit(ctx->green_path);
+ }
+
+ if(ctx->heat_path) {
+ /* Add a break into the heat path geometry and restart it from the
+ * current position. The sampled radiative path becomes a branch of the
+ * current sampled path */
+ res = heat_path_add_break(ctx->heat_path);
+ if(res != RES_OK) goto error;
+
+ res = register_heat_vertex
+ (ctx->heat_path, &rwalk->vtx, T->value, current_vertex_type);
+ if(res != RES_OK) goto error;
+ }
}
/* Null-collision, looping at the beginning */
diff --git a/src/sdis_heat_path_boundary_Xd_solid_solid.h b/src/sdis_heat_path_boundary_Xd_solid_solid.h
@@ -0,0 +1,175 @@
+/* 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_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_scene_c.h"
+
+#include <star/ssp.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Boundary path between a solid and a fluid
+ ******************************************************************************/
+res_T
+XD(solid_solid_boundary_path)
+ (struct sdis_scene* scn,
+ const 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_frt_args =
+ XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
+ struct XD(sample_reinjection_step_args) samp_reinject_step_bck_args =
+ XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
+ struct XD(reinjection_step) reinject_step_frt = XD(REINJECTION_STEP_NULL);
+ struct XD(reinjection_step) reinject_step_bck = XD(REINJECTION_STEP_NULL);
+ struct XD(reinjection_step)* reinject_step = NULL;
+
+ /* Reinjection arguments */
+ struct XD(solid_reinjection_args) solid_reinject_args =
+ XD(SOLID_REINJECTION_ARGS_NULL);
+
+ /* Data attached to the boundary */
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid_frt = NULL;
+ struct sdis_medium* solid_bck = NULL;
+ struct sdis_medium* solid = NULL;
+
+ double lambda_frt;
+ double lambda_bck;
+ double delta_boundary_frt;
+ double delta_boundary_bck;
+
+ double proba;
+ double r;
+ double tcr;
+
+ res_T res = RES_OK;
+ ASSERT(scn && ctx && frag && rwalk && rng && T);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ (void)frag, (void)ctx;
+
+ /* Retrieve the two solids split by the boundary */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ solid_frt = interface_get_medium(interf, SDIS_FRONT);
+ solid_bck = interface_get_medium(interf, SDIS_BACK);
+ ASSERT(solid_frt->type == SDIS_SOLID);
+ ASSERT(solid_bck->type == SDIS_SOLID);
+
+ /* Retrieve the thermal contact resistance */
+ tcr = interface_get_thermal_contact_resistance(interf, frag);
+
+ /* Fetch the properties of the media */
+ lambda_frt = solid_get_thermal_conductivity(solid_frt, &rwalk->vtx);
+ lambda_bck = solid_get_thermal_conductivity(solid_bck, &rwalk->vtx);
+
+ /* 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_frt = solid_get_delta(solid_frt, &rwalk->vtx) * sqrt(DIM);
+ delta_boundary_bck = solid_get_delta(solid_bck, &rwalk->vtx) * sqrt(DIM);
+
+ /* Sample a front/back reinjection steps */
+ samp_reinject_step_frt_args.rng = rng;
+ samp_reinject_step_bck_args.rng = rng;
+ samp_reinject_step_frt_args.solid = solid_frt;
+ samp_reinject_step_bck_args.solid = solid_bck;
+ samp_reinject_step_frt_args.rwalk = rwalk;
+ samp_reinject_step_bck_args.rwalk = rwalk;
+ samp_reinject_step_frt_args.distance = delta_boundary_frt;
+ samp_reinject_step_bck_args.distance = delta_boundary_bck;
+ samp_reinject_step_frt_args.side = SDIS_FRONT;
+ samp_reinject_step_bck_args.side = SDIS_BACK;
+ res = XD(sample_reinjection_step_solid_solid)
+ (scn,
+ &samp_reinject_step_frt_args,
+ &samp_reinject_step_bck_args,
+ &reinject_step_frt,
+ &reinject_step_bck);
+ if(res != RES_OK) goto error;
+
+ r = ssp_rng_canonical(rng);
+ if(tcr == 0) { /* No thermal contact resistance */
+ /* Define the reinjection side. Note that the proba should be : Lf/Df' /
+ * (Lf/Df' + Lb/Db')
+ *
+ * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted
+ * delta of the <front|back> side, i.e. : D<f|b>' =
+ * reinject_dst_<front|back> / sqrt(DIM)
+ *
+ * Anyway, one can avoid to compute the adjusted delta by directly using the
+ * adjusted reinjection distance since the resulting proba is strictly the
+ * same; sqrt(DIM) can be simplified. */
+ const double tmp_frt = lambda_frt / reinject_step_frt.distance;
+ const double tmp_bck = lambda_bck / reinject_step_bck.distance;
+ proba = tmp_frt / (tmp_frt + tmp_bck);
+ } else {
+ const double delta_frt = reinject_step_frt.distance/sqrt(DIM);
+ const double delta_bck = reinject_step_bck.distance/sqrt(DIM);
+ const double tmp_frt = lambda_frt/delta_frt;
+ const double tmp_bck = lambda_bck/delta_bck;
+ const double tmp_tcr = tcr*tmp_frt*tmp_bck;
+ switch(rwalk->hit_side) {
+ case SDIS_BACK:
+ /* When coming from the BACK side, the probability to be reinjected on
+ * the FRONT side depends on the thermal contact resistance: it
+ * decreases when the TCR increases (and tends to 0 when TCR -> +inf) */
+ proba = tmp_frt / (tmp_frt + tmp_bck + tmp_tcr);
+ break;
+ case SDIS_FRONT:
+ /* Same thing when coming from the FRONT side: the probability of
+ * reinjection on the FRONT side depends on the thermal contact
+ * resistance: it increases when the TCR increases (and tends to 1 when
+ * the TCR -> +inf) */
+ proba = (tmp_frt + tmp_tcr) / (tmp_frt + tmp_bck + tmp_tcr);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+
+ if(r < proba) { /* Reinject in front */
+ reinject_step = &reinject_step_frt;
+ solid = solid_frt;
+ } else { /* Reinject in back */
+ reinject_step = &reinject_step_bck;
+ solid = solid_bck;
+ }
+
+ /* Perform the reinjection into the solid */
+ solid_reinject_args.reinjection = reinject_step;
+ solid_reinject_args.rng = rng;
+ solid_reinject_args.rwalk = rwalk;
+ solid_reinject_args.rwalk_ctx = ctx;
+ 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;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -0,0 +1,219 @@
+/* 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/>. */
+
+#ifndef SDIS_HEAT_PATH_BOUNDARY_C_H
+#define SDIS_HEAT_PATH_BOUNDARY_C_H
+
+#include <star/s2d.h>
+#include <star/s3d.h>
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct rwalk_2d;
+struct rwalk_3d;
+struct sdis_scene;
+struct sdis_medium;
+
+/*******************************************************************************
+ * Sample a reinjection step
+ ******************************************************************************/
+struct sample_reinjection_step_args_2d {
+ struct ssp_rng* rng; /* Random number generator to use */
+ const struct sdis_medium* solid; /* Solid in which to reinject */
+ struct rwalk_2d* rwalk; /* Current state of the random walk */
+ double distance; /* Maximum Reinjection distance */
+ enum sdis_side side; /* Side of the boundary to re-inject */
+};
+
+struct sample_reinjection_step_args_3d {
+ struct ssp_rng* rng; /* Random number generator to use */
+ const struct sdis_medium* solid; /* Medium in which to reinject */
+ struct rwalk_3d* rwalk; /* Current random walk state */
+ double distance; /* Maximum Reinjection distance */
+ enum sdis_side side; /* Side of the boundary to re-inject */
+};
+
+struct reinjection_step_2d {
+ struct s2d_hit hit; /* Intersection along the reinjection direction */
+ float direction[2]; /* Reinjection direction */
+ float distance; /* Reinjection distance */
+};
+
+struct reinjection_step_3d {
+ struct s3d_hit hit; /* Intersection along the reinjection direction */
+ float direction[3]; /* Reinjection direction */
+ float distance; /* Reinjection distance */
+};
+
+#define SAMPLE_REINJECTION_STEP_ARGS_NULL___2d \
+ {NULL, NULL, NULL, -1, SDIS_SIDE_NULL__}
+#define SAMPLE_REINJECTION_STEP_ARGS_NULL___3d \
+ {NULL, NULL, NULL, -1, SDIS_SIDE_NULL__}
+static const struct sample_reinjection_step_args_2d
+SAMPLE_REINJECTION_STEP_ARGS_NULL_2d = SAMPLE_REINJECTION_STEP_ARGS_NULL___2d;
+static const struct sample_reinjection_step_args_3d
+SAMPLE_REINJECTION_STEP_ARGS_NULL_3d = SAMPLE_REINJECTION_STEP_ARGS_NULL___3d;
+
+#define REINJECTION_STEP_NULL___2d {S2D_HIT_NULL__, {0,0}, 0}
+#define REINJECTION_STEP_NULL___3d {S3D_HIT_NULL__, {0,0,0}, 0}
+static const struct reinjection_step_2d
+REINJECTION_STEP_NULL_2d = REINJECTION_STEP_NULL___2d;
+static const struct reinjection_step_3d
+REINJECTION_STEP_NULL_3d = REINJECTION_STEP_NULL___3d;
+
+extern LOCAL_SYM res_T
+sample_reinjection_step_solid_fluid_2d
+ (const struct sdis_scene* scn,
+ const struct sample_reinjection_step_args_2d* args,
+ struct reinjection_step_2d* step);
+
+extern LOCAL_SYM res_T
+sample_reinjection_step_solid_fluid_3d
+ (const struct sdis_scene* scn,
+ const struct sample_reinjection_step_args_3d* args,
+ struct reinjection_step_3d *step);
+
+extern LOCAL_SYM res_T
+sample_reinjection_step_solid_solid_2d
+ (const struct sdis_scene* scn,
+ const struct sample_reinjection_step_args_2d* args_front,
+ const struct sample_reinjection_step_args_2d* args_back,
+ struct reinjection_step_2d* step_front,
+ struct reinjection_step_2d* step_back);
+
+extern LOCAL_SYM res_T
+sample_reinjection_step_solid_solid_3d
+ (const struct sdis_scene* scn,
+ const struct sample_reinjection_step_args_3d* args_front,
+ const struct sample_reinjection_step_args_3d* args_back,
+ struct reinjection_step_3d* step_front,
+ struct reinjection_step_3d* step_back);
+
+/*******************************************************************************
+ * Reinject the random walk into a solid
+ ******************************************************************************/
+struct solid_reinjection_args_2d {
+ const struct reinjection_step_2d* reinjection; /* Reinjection to do */
+ const 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;
+ double fp_to_meter;
+};
+
+struct solid_reinjection_args_3d {
+ const struct reinjection_step_3d* reinjection; /* Reinjection to do */
+ const 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;
+ double fp_to_meter;
+};
+
+#define SOLID_REINJECTION_ARGS_NULL___2d {NULL,NULL,NULL,NULL,NULL,0}
+#define SOLID_REINJECTION_ARGS_NULL___3d {NULL,NULL,NULL,NULL,NULL,0}
+static const struct solid_reinjection_args_2d SOLID_REINJECTION_ARGS_NULL_2d =
+ SOLID_REINJECTION_ARGS_NULL___2d;
+static const struct solid_reinjection_args_3d SOLID_REINJECTION_ARGS_NULL_3d =
+ SOLID_REINJECTION_ARGS_NULL___3d;
+
+extern LOCAL_SYM res_T
+solid_reinjection_2d
+ (struct sdis_medium* solid,
+ struct solid_reinjection_args_2d* args);
+
+extern LOCAL_SYM res_T
+solid_reinjection_3d
+ (struct sdis_medium* solid,
+ struct solid_reinjection_args_3d* args);
+
+/*******************************************************************************
+ * Boundary sub-paths
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+solid_boundary_with_flux_path_2d
+ (struct sdis_scene* scn,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ const double phi,
+ struct rwalk_2d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_2d* T);
+
+extern LOCAL_SYM res_T
+solid_boundary_with_flux_path_3d
+ (struct sdis_scene* scn,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ const double phi,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T);
+
+extern LOCAL_SYM res_T
+solid_fluid_boundary_path_2d
+ (struct sdis_scene* scn,
+ const 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
+ (struct sdis_scene* scn,
+ const 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
+ (struct sdis_scene* scn,
+ const 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
+ (struct sdis_scene* scn,
+ const 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_solid_boundary_path_2d
+ (struct sdis_scene* scn,
+ const 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_solid_boundary_path_3d
+ (struct sdis_scene* scn,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T);
+
+#endif /* SDIS_HEAT_PATH_BOUNDARY_C_H */
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -74,8 +74,8 @@ XD(trace_radiative_path)
#endif
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
rwalk->hit_side = SDIS_SIDE_NULL__;
- if(ctx->Tarad >= 0) {
- T->value += ctx->Tarad;
+ if(scn->trad.temperature >= 0) {
+ T->value += scn->trad.temperature;
T->done = 1;
if(ctx->green_path) {
@@ -104,7 +104,7 @@ XD(trace_radiative_path)
"such temperature, one has to setup a valid ambient radiative "
"temperature, i.e. it must be greater or equal to 0.\n",
FUNC_NAME,
- ctx->Tarad,
+ scn->trad.temperature,
SPLIT3(rwalk->vtx.P));
res = RES_BAD_OP;
goto error;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -74,11 +74,13 @@ check_interface_shader
FOR_EACH(i, 0, 2) {
switch(type[i]) {
case SDIS_SOLID:
- if(shaders[i]->emissivity || shaders[i]->specular_fraction) {
+ if(shaders[i]->emissivity
+ || shaders[i]->specular_fraction
+ || shaders[i]->reference_temperature) {
log_warn(dev,
- "%s: the interface side toward a solid can neither have the "
- "emissivity nor the specular_fraction properties. The shader's "
- " pointer functions for these attributes should be NULL.\n",
+ "%s: the interface side toward a solid cannot have an emissivity, "
+ "a specular_fraction or a reference temperature. The shader's "
+ "pointer functions for these attributes should be NULL.\n",
caller_name);
}
break;
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -167,5 +167,21 @@ interface_side_get_specular_fraction
? shader->specular_fraction(frag, interf->data) : 0;
}
+static INLINE double
+interface_side_get_reference_temperature
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ const struct sdis_interface_side_shader* shader;
+ ASSERT(interf && frag);
+ switch(frag->side) {
+ case SDIS_FRONT: shader = &interf->shader.front; break;
+ case SDIS_BACK: shader = &interf->shader.back; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return shader->reference_temperature
+ ? shader->reference_temperature(frag, interf->data) : -1;
+}
+
#endif /* SDIS_INTERFACE_C_H */
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -46,11 +46,10 @@ ray_realisation_3d
rwalk.hit_side = SDIS_SIDE_NULL__;
rwalk.mdm = medium;
- ctx.Tarad = scn->ambient_radiative_temperature;
- ctx.Tref3 = scn->reference_temperature * scn->reference_temperature
- * scn->reference_temperature;
ctx.heat_path = heat_path;
-
+ ctx.That2 = scn->tmax * scn->tmax;
+ ctx.That3 = scn->tmax * ctx.That2;
+
f3_set_d3(dir, direction);
/* Register the starting position against the heat path */
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -181,11 +181,8 @@ XD(probe_realisation)
ctx.green_path = green_path;
ctx.heat_path = heat_path;
- ctx.Tarad = scn->ambient_radiative_temperature;
- ctx.Tref3 =
- scn->reference_temperature
- * scn->reference_temperature
- * scn->reference_temperature;
+ ctx.That2 = scn->tmax * scn->tmax;
+ ctx.That3 = scn->tmax * ctx.That2;
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
@@ -259,9 +256,8 @@ XD(boundary_realisation)
ctx.green_path = green_path;
ctx.heat_path = heat_path;
- ctx.Tarad = scn->ambient_radiative_temperature;
- ctx.Tref3 = scn->reference_temperature * scn->reference_temperature
- * scn->reference_temperature;
+ ctx.That2 = scn->tmax * scn->tmax;
+ ctx.That3 = scn->tmax * ctx.That2;
res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
@@ -300,8 +296,9 @@ XD(boundary_flux_realisation)
#endif
double P[SDIS_XD_DIMENSION];
float N[SDIS_XD_DIMENSION];
- const double Tr3 = scn->reference_temperature * scn->reference_temperature
- * scn->reference_temperature;
+ const double That = scn->tmax;
+ const double That2 = That * That;
+ const double That3 = That * That2;
const enum sdis_side fluid_side =
(solid_side == SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT;
res_T res = RES_OK;
@@ -336,8 +333,8 @@ XD(boundary_flux_realisation)
rwalk.mdm = (Mdm); \
rwalk.hit.prim = prim; \
SET_PARAM(rwalk.hit, st); \
- ctx.Tarad = scn->ambient_radiative_temperature; \
- ctx.Tref3 = Tr3; \
+ ctx.That2 = That2; \
+ ctx.That3 = That3; \
dX(set)(rwalk.vtx.P, P); \
fX(set)(rwalk.hit.normal, N); \
T = XD(TEMPERATURE_NULL); \
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -203,40 +203,36 @@ sdis_scene_set_fp_to_meter
res_T
sdis_scene_get_ambient_radiative_temperature
(const struct sdis_scene* scn,
- double* trad)
+ struct sdis_ambient_radiative_temperature* trad)
{
if(!scn || !trad) return RES_BAD_ARG;
- *trad = scn->ambient_radiative_temperature;
+ *trad = scn->trad;
return RES_OK;
}
res_T
-sdis_scene_set_reference_temperature
+sdis_scene_set_ambient_radiative_temperature
(struct sdis_scene* scn,
- const double tref)
+ const struct sdis_ambient_radiative_temperature* trad)
{
- if(!scn || tref < 0) return RES_BAD_ARG;
- scn->reference_temperature = tref;
+ if(!scn) return RES_BAD_ARG;
+ scn->trad = *trad;
return RES_OK;
}
res_T
-sdis_scene_get_reference_temperature
- (const struct sdis_scene* scn,
- double* tref)
+sdis_scene_get_maximum_temperature(const struct sdis_scene* scn, double* tmax)
{
- if(!scn || !tref) return RES_BAD_ARG;
- *tref = scn->reference_temperature;
+ if(!scn || !tmax) return RES_BAD_ARG;
+ *tmax = scn->tmax;
return RES_OK;
}
res_T
-sdis_scene_set_ambient_radiative_temperature
- (struct sdis_scene* scn,
- const double trad)
+sdis_scene_set_maximum_temperature(struct sdis_scene* scn, const double tmax)
{
- if(!scn) return RES_BAD_ARG;
- scn->ambient_radiative_temperature = trad;
+ if(!scn || tmax < 0) return RES_BAD_ARG;
+ scn->tmax = tmax;
return RES_OK;
}
@@ -474,7 +470,8 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
} else {
S3D(scene_view_primitives_count(scn->s3d_view, &nprims));
}
- WRITE(&scn->reference_temperature, 1);
+ WRITE(&scn->trad.reference, 1);
+ WRITE(&scn->tmax, 1);
WRITE(&scn->fp_to_meter, 1);
FOR_EACH(iprim, 0, nprims) {
struct sdis_interface* interf = NULL;
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -912,8 +912,8 @@ XD(scene_create)
SDIS(device_ref_get(dev));
scn->dev = dev;
scn->fp_to_meter = args->fp_to_meter;
- scn->ambient_radiative_temperature = args->trad;
- scn->reference_temperature = args->tref;
+ scn->trad = args->trad;
+ scn->tmax = args->tmax;
scn->outer_enclosure_id = UINT_MAX;
darray_interf_init(dev->allocator, &scn->interfaces);
darray_medium_init(dev->allocator, &scn->media);
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -209,8 +209,8 @@ struct sdis_scene {
unsigned outer_enclosure_id;
double fp_to_meter;
- double ambient_radiative_temperature; /* In Kelvin */
- double reference_temperature;
+ struct sdis_ambient_radiative_temperature trad;
+ double tmax; /* Maximum temperature of the system (In Kelvin) */
ref_T ref;
struct sdis_device* dev;
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -573,8 +573,7 @@ XD(solve_boundary_flux)
const struct sdis_medium *fmd, *bmd;
enum sdis_side solid_side, fluid_side;
struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
- const double Tref = scn->reference_temperature;
- double epsilon, hc, hr, imposed_flux, imposed_temp;
+ double epsilon, hc, hr, imposed_flux, imposed_temp, Tref;
size_t iprim;
double uv[DIM - 1];
float st[DIM - 1];
@@ -638,6 +637,7 @@ XD(solve_boundary_flux)
/* Fetch interface parameters */
epsilon = interface_side_get_emissivity(interf, &frag);
+ Tref = interface_side_get_reference_temperature(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
frag.side = solid_side;
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -480,7 +480,7 @@ XD(solve_probe_boundary_flux)
double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
- const double Tref = scn->reference_temperature;
+ double Tref = -1;
size_t n;
int pcent;
res_T res_simul = RES_OK;
@@ -496,6 +496,7 @@ XD(solve_probe_boundary_flux)
frag.time = time;
frag.side = fluid_side;
epsilon = interface_side_get_emissivity(interf, &frag);
+ Tref = interface_side_get_reference_temperature(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
frag.side = solid_side;
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -69,7 +69,7 @@ static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = {
-1.5, 1.0, 1.0,
1.5, 1.0, 1.0,
};
-static const size_t nvertices = sizeof(vertices) / (3*sizeof(double));
+static const size_t nvertices = sizeof(vertices) / (sizeof(double)*3);
static const size_t indices[32/*#triangles*/*3/*#indices per triangle*/] = {
0, 2, 1, 1, 2, 3, /* Solid back face */
@@ -91,7 +91,7 @@ static const size_t indices[32/*#triangles*/*3/*#indices per triangle*/] = {
3, 7, 11, 11, 7, 15, /* Right fluid top face */
1, 9, 5, 5, 9, 13 /* Right fluid bottom face */
};
-static const size_t ntriangles = sizeof(indices) / (3*sizeof(size_t));
+static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3);
static void
get_indices(const size_t itri, size_t ids[3], void* ctx)
@@ -193,6 +193,7 @@ struct interfac {
double convection_coef;
double emissivity;
double specular_fraction;
+ double Tref;
};
static double
@@ -227,6 +228,14 @@ interface_get_specular_fraction
return ((const struct interfac*)sdis_data_cget(data))->specular_fraction;
}
+static double
+interface_get_Tref
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interfac*)sdis_data_cget(data))->Tref;
+}
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -251,10 +260,12 @@ create_interface
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->convection_coef);
@@ -337,6 +348,7 @@ main(int argc, char** argv)
interf.convection_coef = -1;
interf.emissivity = -1;
interf.specular_fraction = -1;
+ interf.Tref = Tref;
create_interface(dev, solid, solid2, &interf, interfaces+0);
/* Create the interface that emits radiative heat from the solid */
@@ -344,6 +356,7 @@ main(int argc, char** argv)
interf.convection_coef = 0;
interf.emissivity = emissivity;
interf.specular_fraction = 1;
+ interf.Tref = Tref;
create_interface(dev, solid, fluid, &interf, interfaces+1);
/* Create the interface that forces the radiative heat to bounce */
@@ -351,6 +364,7 @@ main(int argc, char** argv)
interf.convection_coef = 0;
interf.emissivity = 0;
interf.specular_fraction = 1;
+ interf.Tref = Tref;
create_interface(dev, fluid, solid2, &interf, interfaces+2);
/* Create the interface with a limit condition of T0 Kelvin */
@@ -358,6 +372,7 @@ main(int argc, char** argv)
interf.convection_coef = 0;
interf.emissivity = 1;
interf.specular_fraction = 1;
+ interf.Tref = T0;
create_interface(dev, fluid, solid2, &interf, interfaces+3);
/* Create the interface with a limit condition of T1 Kelvin */
@@ -365,6 +380,7 @@ main(int argc, char** argv)
interf.convection_coef = 0;
interf.emissivity = 1;
interf.specular_fraction = 1;
+ interf.Tref = T1;
create_interface(dev, fluid, solid2, &interf, interfaces+4);
/* Setup the per primitive interface of the solid medium */
@@ -398,7 +414,7 @@ main(int argc, char** argv)
scn_args.get_position = get_position;
scn_args.nprimitives = ntriangles;
scn_args.nvertices = nvertices;
- scn_args.tref = Tref;
+ scn_args.tmax = MMAX(T0, T1);
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, &scn));
@@ -476,6 +492,7 @@ 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)));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -58,7 +58,7 @@ static const double vertices[8/*#vertices*/*2/*#coords par vertex*/] = {
-1.5, 1.0,
1.5, 1.0
};
-static const size_t nvertices = sizeof(vertices) / (2*sizeof(double));
+static const size_t nvertices = sizeof(vertices) / (sizeof(double)*2);
static const size_t indices[10/*#segments*/*2/*#indices per segment*/] = {
0, 1, /* Solid bottom segment */
@@ -74,7 +74,7 @@ static const size_t indices[10/*#segments*/*2/*#indices per segment*/] = {
3, 7, /* Right fluid top segment */
7, 4 /* Right fluid right segment */
};
-static const size_t nsegments = sizeof(indices) / (2*sizeof(size_t));
+static const size_t nsegments = sizeof(indices) / (sizeof(size_t)*2);
static void
get_indices(const size_t iseg, size_t ids[2], void* ctx)
@@ -158,11 +158,12 @@ struct interfac {
double temperature;
double emissivity;
double specular_fraction;
+ double reference_temperature;
} front, back;
};
static const struct interfac INTERFACE_NULL = {
- 0, {-1, -1, -1}, {-1, -1, -1}
+ 0, {-1, -1, -1, -1}, {-1, -1, -1, -1}
};
static double
@@ -223,6 +224,22 @@ interface_get_specular_fraction
return f;
}
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interfac* interf;
+ double T = -1;
+ CHK(data != NULL && frag != NULL);
+ interf = sdis_data_cget(data);
+ switch(frag->side) {
+ case SDIS_FRONT: T = interf->front.reference_temperature; break;
+ case SDIS_BACK: T = interf->back.reference_temperature; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return T;
+}
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -250,10 +267,12 @@ create_interface
if(type_f == SDIS_FLUID) {
shader.front.emissivity = interface_get_emissivity;
shader.front.specular_fraction = interface_get_specular_fraction;
+ shader.front.reference_temperature = interface_get_reference_temperature;
}
if(type_b == SDIS_FLUID) {
shader.back.emissivity = interface_get_emissivity;
shader.back.specular_fraction = interface_get_specular_fraction;
+ shader.back.reference_temperature = interface_get_reference_temperature;
}
shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef);
@@ -336,6 +355,7 @@ main(int argc, char** argv)
interf.back.temperature = UNKNOWN_TEMPERATURE;
interf.back.emissivity = emissivity;
interf.back.specular_fraction = -1; /* Should not be fetched */
+ interf.back.reference_temperature = Tref;
create_interface(dev, solid, fluid, &interf, interfaces+1);
/* Create the interface that forces the radiative heat to bounce */
@@ -343,6 +363,7 @@ main(int argc, char** argv)
interf.front.temperature = UNKNOWN_TEMPERATURE;
interf.front.emissivity = 0;
interf.front.specular_fraction = 1;
+ interf.front.reference_temperature = Tref;
create_interface(dev, fluid, solid2, &interf, interfaces+2);
/* Create the interface with a limit condition of T0 Kelvin */
@@ -350,6 +371,7 @@ main(int argc, char** argv)
interf.front.temperature = T0;
interf.front.emissivity = 1;
interf.front.specular_fraction = 1;
+ interf.front.reference_temperature = T0;
create_interface(dev, fluid, solid2, &interf, interfaces+3);
/* Create the interface with a limit condition of T1 Kelvin */
@@ -357,6 +379,7 @@ main(int argc, char** argv)
interf.front.temperature = T1;
interf.front.emissivity = 1;
interf.front.specular_fraction = 1;
+ interf.front.reference_temperature = T1;
create_interface(dev, fluid, solid2, &interf, interfaces+4);
/* Setup the per primitive interface of the solid medium */
@@ -384,8 +407,8 @@ main(int argc, char** argv)
scn_args.get_position = get_position;
scn_args.nprimitives = nsegments;
scn_args.nvertices = nvertices;
- scn_args.tref = Tref;
scn_args.context = &geom;
+ scn_args.tmax = MMAX(T0, T1);
OK(sdis_scene_2d_create(dev, &scn_args, &scn));
hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
diff --git a/src/test_sdis_contact_resistance.h b/src/test_sdis_contact_resistance.h
@@ -37,8 +37,7 @@
/*******************************************************************************
* Box geometry
******************************************************************************/
-static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/]
-= {
+static const double model3d_vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
0, 0, 0,
X0, 0, 0,
L, 0, 0,
@@ -52,7 +51,7 @@ static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/]
X0, L, L,
L, L, L
};
-static const size_t model3d_nvertices = sizeof(model3d_vertices) / (3*sizeof(double));
+static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)*3);
/* The following array lists the indices toward the 3D vertices of each
* triangle.
@@ -64,8 +63,7 @@ static const size_t model3d_nvertices = sizeof(model3d_vertices) / (3*sizeof(dou
* 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*/]
-= {
+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 */
0, 6, 3, 3, 6, 9, /* -X */
6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */
@@ -74,7 +72,7 @@ static const size_t model3d_indices[22/*#triangles*/ * 3/*#indices per triangle*
0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */
4, 10, 7, 7, 1, 4 /* Inside */
};
-static const size_t model3d_ntriangles = sizeof(model3d_indices) / (3*sizeof(size_t));
+static const size_t model3d_ntriangles = sizeof(model3d_indices)/(sizeof(size_t)*3);
static INLINE void
model3d_get_indices(const size_t itri, size_t ids[3], void* context)
@@ -110,7 +108,7 @@ model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* co
/*******************************************************************************
* Square geometry
******************************************************************************/
-static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = {
+static const double model2d_vertices[6/*#vertices*/*2/*#coords per vertex*/] = {
L, 0,
X0, 0,
0, 0,
@@ -118,7 +116,7 @@ static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] =
X0, L,
L, L
};
-static const size_t model2d_nvertices = sizeof(model2d_vertices) / (2*sizeof(double));
+static const size_t model2d_nvertices = sizeof(model2d_vertices)/(sizeof(double)*2);
static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = {
0, 1, 1, 2, /* Bottom */
@@ -127,8 +125,7 @@ static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/]
5, 0, /* Right */
4, 1 /* Inside */
};
-static const size_t model2d_nsegments = sizeof(model2d_indices) / (2*sizeof(size_t));
-
+static const size_t model2d_nsegments = sizeof(model2d_indices) / (sizeof(size_t)*2);
static INLINE void
model2d_get_indices(const size_t iseg, size_t ids[2], void* context)
diff --git a/src/test_sdis_picard1.c b/src/test_sdis_picard1.c
@@ -0,0 +1,722 @@
+/* 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
@@ -258,6 +258,8 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
double duplicated_vertices[] = { 0, 0, 0, 0 };
struct sdis_scene* scn = NULL;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_ambient_radiative_temperature trad =
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
double lower[2], upper[2];
double u0, u1, u2, pos[2], pos1[2];
double dst, fp, t;
@@ -353,28 +355,32 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
- BA(sdis_scene_get_ambient_radiative_temperature(NULL, &t));
- OK(sdis_scene_get_ambient_radiative_temperature(scn, &t));
- CHK(t == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad);
-
- t = 100;
- BA(sdis_scene_set_ambient_radiative_temperature(NULL, t));
- OK(sdis_scene_set_ambient_radiative_temperature(scn, t));
- OK(sdis_scene_get_ambient_radiative_temperature(scn, &t));
- CHK(t == 100);
-
- BA(sdis_scene_get_reference_temperature(NULL, NULL));
- BA(sdis_scene_get_reference_temperature(scn, NULL));
- BA(sdis_scene_get_reference_temperature(NULL, &t));
- OK(sdis_scene_get_reference_temperature(scn, &t));
- CHK(t == SDIS_SCENE_CREATE_ARGS_DEFAULT.tref);
+ BA(sdis_scene_get_ambient_radiative_temperature(NULL, &trad));
+ OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
+ CHK(trad.temperature == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature);
+ CHK(trad.reference == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference);
+
+ trad.temperature = 100;
+ trad.reference = 110;
+ BA(sdis_scene_set_ambient_radiative_temperature(NULL, &trad));
+ OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
+ trad = SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
+ OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
+ 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_reference_temperature(NULL, t));
- BA(sdis_scene_set_reference_temperature(scn, t));
+ BA(sdis_scene_set_maximum_temperature(NULL, t));
+ BA(sdis_scene_set_maximum_temperature(scn, t));
t = 100;
- OK(sdis_scene_set_reference_temperature(scn, t));
- OK(sdis_scene_get_reference_temperature(scn, &t));
+ OK(sdis_scene_set_maximum_temperature(scn, t));
+ OK(sdis_scene_get_maximum_temperature(scn, &t));
CHK(t == 100);
BA(sdis_scene_get_boundary_position(NULL, 1, &u0, pos));
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -142,6 +142,7 @@ struct interf {
double temperature;
double emissivity;
double hc;
+ double reference_temperature;
};
static double
@@ -171,6 +172,15 @@ interface_get_convection_coef
return interf->hc;
}
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->reference_temperature;
+}
+
/*******************************************************************************
* Helper function
******************************************************************************/
@@ -289,7 +299,9 @@ main(int argc, char** argv)
interf_props->hc = H;
interf_props->temperature = Tb;
interf_props->emissivity = EPSILON;
+ interf_props->reference_temperature = Tb;
interf_shader.back.emissivity = interface_get_emissivity;
+ interf_shader.back.reference_temperature = interface_get_reference_temperature;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_Tb));
interf_shader.back.emissivity = NULL;
@@ -301,7 +313,9 @@ main(int argc, char** argv)
interf_props->hc = H;
interf_props->temperature = UNKNOWN_TEMPERATURE;
interf_props->emissivity = EPSILON;
+ interf_props->reference_temperature = Tref;
interf_shader.back.emissivity = interface_get_emissivity;
+ interf_shader.back.reference_temperature = interface_get_reference_temperature;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_H));
interf_shader.back.emissivity = NULL;
@@ -331,8 +345,9 @@ main(int argc, char** argv)
scn_args.get_position = box_get_position;
scn_args.nprimitives = box_ntriangles;
scn_args.nvertices = box_nvertices;
- scn_args.trad = Trad;
- scn_args.tref = Tref;
+ scn_args.trad.temperature = Trad;
+ scn_args.trad.reference = Trad;
+ scn_args.tmax = MMAX(MMAX(Tf, Trad), Tb);
scn_args.context = box_interfaces;
OK(sdis_scene_create(dev, &scn_args, &box_scn));
@@ -342,8 +357,9 @@ main(int argc, char** argv)
scn_args.get_position = square_get_position;
scn_args.nprimitives = square_nsegments;
scn_args.nvertices = square_nvertices;
- scn_args.trad = Trad;
- scn_args.tref = Tref;
+ scn_args.trad.temperature = Trad;
+ scn_args.trad.reference = Trad;
+ scn_args.tmax = 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
@@ -237,8 +237,11 @@ struct interf {
double epsilon;
double specular_fraction;
double temperature;
+ double reference_temperature;
+};
+static const struct interf INTERF_NULL = {
+ 0, 0, 0, UNKOWN_TEMPERATURE, UNKOWN_TEMPERATURE
};
-static const struct interf INTERF_NULL = {0, 0, 0, UNKOWN_TEMPERATURE};
static double
interface_get_convection_coef
@@ -272,6 +275,14 @@ interface_get_temperature
return ((const struct interf*)sdis_data_cget(data))->temperature;
}
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->reference_temperature;
+}
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -372,10 +383,14 @@ create_interface
if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) {
interface_shader.front.emissivity = interface_get_emissivity;
interface_shader.front.specular_fraction = interface_get_specular_fraction;
+ interface_shader.front.reference_temperature =
+ interface_get_reference_temperature;
}
if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) {
interface_shader.back.emissivity = interface_get_emissivity;
interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ interface_shader.back.reference_temperature =
+ interface_get_reference_temperature;
}
/* Create the interface */
OK(sdis_interface_create
@@ -542,6 +557,8 @@ main(int argc, char** argv)
struct sdis_scene* scn = NULL;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
struct sdis_solve_camera_args solve_args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
+ struct sdis_ambient_radiative_temperature trad =
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
struct ssp_rng* rng_state = NULL;
struct fluid fluid_param = FLUID_NULL;
struct solid solid_param = SOLID_NULL;
@@ -590,6 +607,7 @@ main(int argc, char** argv)
interface_param.epsilon = 1;
interface_param.specular_fraction = 1;
interface_param.temperature = UNKOWN_TEMPERATURE;
+ interface_param.reference_temperature = 300;
create_interface(dev, fluid1, solid, &interface_param, &interf1);
/* Setup the cube geometry */
@@ -614,8 +632,9 @@ main(int argc, char** argv)
scn_args.get_position = geometry_get_position;
scn_args.nprimitives = ntris;
scn_args.nvertices = npos;
- scn_args.trad = 300;
- scn_args.tref = 300;
+ scn_args.trad.temperature = 300;
+ scn_args.trad.reference = 300;
+ scn_args.tmax = 350;
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, &scn));
@@ -650,9 +669,12 @@ main(int argc, char** argv)
solve_args.cam = NULL;
BA(sdis_solve_camera(scn, &solve_args, &buf));
solve_args.cam = cam;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, -1));
+ OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
+ trad.temperature = -1;
+ OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
BA(sdis_solve_camera(scn, &solve_args, &buf));
- OK(sdis_scene_set_ambient_radiative_temperature(scn, 300));
+ trad.temperature = 300;
+ OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
solve_args.time_range[0] = solve_args.time_range[1] = -1;
BA(sdis_solve_camera(scn, &solve_args, &buf));
solve_args.time_range[0] = 1;
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -141,6 +141,7 @@ struct interf {
double hc;
double epsilon;
double specular_fraction;
+ double reference_temperature;
};
static double
@@ -167,6 +168,14 @@ interface_get_specular_fraction
return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
}
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->reference_temperature;
+}
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -267,6 +276,8 @@ main(int argc, char** argv)
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_ambient_radiative_temperature trad =
+ SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
struct dump_path_context dump_ctx = DUMP_PATH_CONTEXT_NULL;
struct context ctx;
struct fluid* fluid_param;
@@ -284,7 +295,7 @@ main(int argc, char** argv)
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
- OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
/* Create the fluid medium */
OK(sdis_data_create
@@ -324,6 +335,7 @@ main(int argc, char** argv)
interface_shader.back.temperature = NULL;
interface_shader.back.emissivity = interface_get_emissivity;
interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ interface_shader.back.reference_temperature = interface_get_reference_temperature;
OK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &interf));
OK(sdis_data_ref_put(data));
@@ -538,10 +550,12 @@ main(int argc, char** argv)
/* Green and ambient radiative temperature */
solve_args.nrealisations = N;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, 300));
- OK(sdis_scene_set_reference_temperature(scn, 300));
+ trad.temperature = trad.reference = 300;
+ OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
+ OK(sdis_scene_set_maximum_temperature(scn, 300));
interface_param->epsilon = 1;
+ interface_param->reference_temperature = 300;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_solve_probe_green_function(scn, &solve_args, &green));
@@ -554,7 +568,9 @@ main(int argc, char** argv)
OK(sdis_estimator_ref_put(estimator2));
/* Check same green used at different ambient radiative temperature */
- OK(sdis_scene_set_ambient_radiative_temperature(scn, 600));
+ trad.temperature = 600;
+ OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
+ OK(sdis_scene_set_maximum_temperature(scn, 600));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_green_function_solve(green, &estimator2));
diff --git a/src/test_sdis_transcient.c b/src/test_sdis_transcient.c
@@ -47,7 +47,7 @@ static const double vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
1.0, 0.0, 1.0,
1.0, 1.0, 1.0
};
-static const size_t nvertices = sizeof(vertices) / (3*sizeof(double));
+static const size_t nvertices = sizeof(vertices) / (sizeof(double)*3);
/* The following array lists the indices toward the 3D vertices of each
* triangle.
@@ -67,7 +67,7 @@ static const size_t indices[22/*#triangles*/*3/*#indices per triangle*/] = {
0, 2, 1, 1, 2, 3, 1, 3, 8, 8, 3, 9, /* Z min */
4, 5, 6, 6, 5, 7, 5,10, 7, 7,10,11 /* Z max */
};
-static const size_t ntriangles = sizeof(indices) / (3*sizeof(size_t));
+static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3);
/*******************************************************************************
* Box geometry functions
diff --git a/src/test_sdis_unstationary_atm.c b/src/test_sdis_unstationary_atm.c
@@ -23,6 +23,20 @@
#include <star/ssp.h>
/*
+ * The physical configuration is the following: a slab of fluid with known
+ * thermophysical properties but unknown temperature is located between a
+ * "ground" and a slab of solid, with also a unknown temperature profile. On
+ * the other side of the solid slab, is a "atmosphere" with known temperature,
+ * and known radiative temperature.
+ *
+ * Solving the system means: finding the temperature of the ground, of the
+ * fluid, of the boundaries, and also the temperature inside the solid, at
+ * various locations (the 1D slab is discretized in order to obtain the
+ * reference)
+ *
+ * The reference for this system comes from a numerical method and is not
+ * analytic. Thus the compliance test MC VS reference is not the usual |MC -
+ * ref| <= 3*sigma but is |MC -ref| <= (Tmax -Tmin) * 0.01.
*
* 3D 2D
*
@@ -59,6 +73,10 @@
#define HA 400
#define TR 260
+#define TMAX (MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), MMAX(TG, TA)), TR))
+#define TMIN (MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), MMIN(TG, TA)), TR))
+#define EPS ((TMAX-TMIN)*0.01)
+
/* hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon
* Tref = (hr / (4 * 5.6696e-8 * epsilon)) ^ 1/3, hr = 6 */
#define TREF 297.974852286
@@ -73,11 +91,10 @@
#define DELTA (XE/30.0)
- /*******************************************************************************
- * Box geometry
- ******************************************************************************/
-static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/]
-= {
+/*******************************************************************************
+ * Box geometry
+ ******************************************************************************/
+static const double model3d_vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
0, 0, 0,
XH, 0, 0,
XHpE, 0, 0,
@@ -91,7 +108,7 @@ static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/]
XH, XHpE, XHpE,
XHpE, XHpE, XHpE
};
-static const size_t model3d_nvertices = sizeof(model3d_vertices) / (3*sizeof(double));
+static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)*3);
/* The following array lists the indices toward the 3D vertices of each
* triangle.
@@ -103,8 +120,7 @@ static const size_t model3d_nvertices = sizeof(model3d_vertices) / (3*sizeof(dou
* 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*/]
-= {
+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 */
0, 6, 3, 3, 6, 9, /* -X */
6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */
@@ -113,7 +129,7 @@ static const size_t model3d_indices[22/*#triangles*/ * 3/*#indices per triangle*
0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */
4, 10, 7, 7, 1, 4 /* Inside */
};
-static const size_t model3d_ntriangles = sizeof(model3d_indices) / (3*sizeof(size_t));
+static const size_t model3d_ntriangles = sizeof(model3d_indices)/(sizeof(size_t)*3);
static INLINE void
model3d_get_indices(const size_t itri, size_t ids[3], void* context)
@@ -157,7 +173,7 @@ static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] =
XH, XHpE,
XHpE, XHpE
};
-static const size_t model2d_nvertices = sizeof(model2d_vertices) / (2*sizeof(double));
+static const size_t model2d_nvertices = sizeof(model2d_vertices)/(sizeof(double)*2);
static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = {
0, 1, 1, 2, /* Bottom */
@@ -166,8 +182,7 @@ static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/]
5, 0, /* Right */
4, 1 /* Inside */
};
-static const size_t model2d_nsegments = sizeof(model2d_indices) / (2*sizeof(size_t));
-
+static const size_t model2d_nsegments = sizeof(model2d_indices)/(sizeof(size_t)*2);
static INLINE void
model2d_get_indices(const size_t iseg, size_t ids[2], void* context)
@@ -309,6 +324,7 @@ struct interf {
double temperature;
double emissivity;
double h;
+ double Tref;
};
static double
@@ -341,10 +357,56 @@ interface_get_emissivity
return interf->emissivity;
}
+static double
+interface_get_Tref
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf;
+ CHK(frag && data);
+ interf = sdis_data_cget(data);
+ return interf->Tref;
+}
+
/*******************************************************************************
* Helper functions
******************************************************************************/
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;
+ shader.convection_coef_upper_bound = interf->h;
+ }
+ if(sdis_medium_get_type(front) == SDIS_FLUID) {
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.reference_temperature = interface_get_Tref;
+ }
+ if(sdis_medium_get_type(back) == SDIS_FLUID) {
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.reference_temperature = interface_get_Tref;
+ }
+
+ OK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ *((struct interf*)sdis_data_get(data)) = *interf;
+
+ OK(sdis_interface_create(dev, front, back, &shader, data, out_interf));
+ OK(sdis_data_ref_put(data));
+}
+
+static void
solve_tbound1
(struct sdis_scene* scn,
struct ssp_rng* rng)
@@ -415,7 +477,7 @@ solve_tbound1
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- CHK(t[isimul] == 0 || eq_eps(T.E, ref[isimul], T.SE * 4));
+ CHK(eq_eps(T.E, ref[isimul], EPS));
OK(sdis_estimator_ref_put(estimator));
}
@@ -492,7 +554,7 @@ solve_tbound2
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- CHK(eq_eps(T.E, ref[isimul], T.SE * 4));
+ CHK(eq_eps(T.E, ref[isimul], EPS));
OK(sdis_estimator_ref_put(estimator));
}
@@ -558,7 +620,7 @@ solve_tsolid
CHK(nfails + nreals == N);
CHK(nfails <= N / 1000);
- CHK(eq_eps(T.E, ref[isimul], T.SE * 4));
+ CHK(eq_eps(T.E, ref[isimul], EPS));
OK(sdis_estimator_ref_put(estimator));
}
@@ -577,6 +639,7 @@ solve_tfluid
size_t nreals;
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,
@@ -621,7 +684,9 @@ solve_tfluid
CHK(nfails + nreals == N);
CHK(nfails <= N / 1000);
- CHK(eq_eps(T.E, ref[isimul], T.SE * 4));
+
+ eps = EPS;
+ CHK(eq_eps(T.E, ref[isimul], eps));
OK(sdis_estimator_ref_put(estimator));
}
@@ -650,10 +715,9 @@ main(int argc, char** argv)
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_interface* model3d_interfaces[22 /*#triangles*/];
struct sdis_interface* model2d_interfaces[7/*#segments*/];
- struct interf* interf_props = NULL;
+ struct interf interf_props;
struct solid* solid_props = NULL;
struct fluid* fluid_props = NULL;
struct ssp_rng* rng = NULL;
@@ -718,66 +782,34 @@ main(int argc, char** argv)
OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid_A));
OK(sdis_data_ref_put(data));
- /* Setup the interface shader */
- interf_shader.front.temperature = interface_get_temperature;
- interf_shader.back.temperature = interface_get_temperature;
- interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.convection_coef_upper_bound = 0;
-
/* Create the adiabatic interfaces */
- OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
- interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- interf_props->h = 0;
- OK(sdis_interface_create
- (dev, fluid, dummy_solid, &interf_shader, data, &interf_adiabatic_1));
- OK(sdis_data_ref_put(data));
-
- interf_shader.convection_coef = NULL;
- OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
- interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- interf_props->h = 0;
- OK(sdis_interface_create
- (dev, solid, dummy_solid, &interf_shader, data, &interf_adiabatic_2));
- OK(sdis_data_ref_put(data));
+ interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.h = 0;
+ interf_props.emissivity = 0;
+ interf_props.Tref = TREF;
+ create_interface(dev, fluid, dummy_solid, &interf_props, &interf_adiabatic_1);
+ create_interface(dev, solid, dummy_solid, &interf_props, &interf_adiabatic_2);
/* Create the P interface */
- interf_shader.convection_coef_upper_bound = HC;
- interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.front.emissivity = interface_get_emissivity;
- OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
- interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- interf_props->h = HC;
- interf_props->emissivity = 1;
- OK(sdis_interface_create
- (dev, fluid, solid, &interf_shader, data, &interf_P));
- OK(sdis_data_ref_put(data));
+ interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.h = HC;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TREF;
+ create_interface(dev, fluid, solid, &interf_props, &interf_P);
/* Create the TG interface */
- interf_shader.convection_coef_upper_bound = HG;
- OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
- interf_props = sdis_data_get(data);
- interf_props->temperature = TG;
- interf_props->h = HG;
- interf_props->emissivity = 1;
- OK(sdis_interface_create
- (dev, fluid, dummy_solid, &interf_shader, data, &interf_TG));
- OK(sdis_data_ref_put(data));
+ interf_props.temperature = TG;
+ interf_props.h = HG;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TG;
+ create_interface(dev, fluid, dummy_solid, &interf_props, &interf_TG);
/* Create the TA interface */
- interf_shader.convection_coef_upper_bound = HA;
- interf_shader.front.emissivity = NULL;
- interf_shader.back.emissivity = interface_get_emissivity;
- OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
- interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- interf_props->h = HA;
- interf_props->emissivity = 1;
- OK(sdis_interface_create
- (dev, solid, fluid_A, &interf_shader, data, &interf_TA));
- OK(sdis_data_ref_put(data));
+ interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.h = HA;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TREF;
+ create_interface(dev, solid, fluid_A, &interf_props, &interf_TA);
/* Release the media */
OK(sdis_medium_ref_put(solid));
@@ -835,8 +867,12 @@ main(int argc, char** argv)
scn_args.nprimitives = model3d_ntriangles;
scn_args.nvertices = model3d_nvertices;
scn_args.context = model3d_interfaces;
- scn_args.trad = TR;
- scn_args.tref = TREF;
+ 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);
OK(sdis_scene_create(dev, &scn_args, &box_scn));
/* Create the square scene */
@@ -846,8 +882,12 @@ main(int argc, char** argv)
scn_args.nprimitives = model2d_nsegments;
scn_args.nvertices = model2d_nvertices;
scn_args.context = model2d_interfaces;
- scn_args.trad = TR;
- scn_args.tref = TREF;
+ 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);
OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
/* Release the interfaces */
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -125,6 +125,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_get_limit_point(NULL, &pt));
BA(sdis_green_path_get_limit_point(path, NULL));
if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) {
+ struct sdis_ambient_radiative_temperature trad;
struct sdis_green_function* green;
struct sdis_scene* scn;
BO(sdis_green_path_get_limit_point(path, &pt));
@@ -140,8 +141,9 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
- BA(sdis_scene_get_ambient_radiative_temperature(NULL, &temp));
- OK(sdis_scene_get_ambient_radiative_temperature(scn, &temp));
+ BA(sdis_scene_get_ambient_radiative_temperature(NULL, &trad));
+ OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
+ temp = trad.temperature;
} else {
OK(sdis_green_path_get_limit_point(path, &pt));
switch(pt.type) {
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -41,7 +41,7 @@ static const double box_vertices[8/*#vertices*/*3/*#coords per vertex*/] = {
0.0, 1.0, 1.0,
1.0, 1.0, 1.0
};
-static const size_t box_nvertices = sizeof(box_vertices) / (3*sizeof(double));
+static const size_t box_nvertices = sizeof(box_vertices) / (sizeof(double)*3);
/* The following array lists the indices toward the 3D vertices of each
* triangle.
@@ -61,7 +61,7 @@ static const size_t box_indices[12/*#triangles*/*3/*#indices per triangle*/] = {
2, 6, 7, 7, 3, 2, /* +Y */
0, 1, 5, 5, 4, 0 /* -Y */
};
-static const size_t box_ntriangles = sizeof(box_indices) / (3*sizeof(size_t));
+static const size_t box_ntriangles = sizeof(box_indices) / (sizeof(size_t)*3);
static INLINE void
box_get_indices(const size_t itri, size_t ids[3], void* context)
@@ -103,7 +103,7 @@ static const double square_vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
0.0, 1.0,
1.0, 1.0
};
-static const size_t square_nvertices = sizeof(square_vertices)/(2*sizeof(double));
+static const size_t square_nvertices = sizeof(square_vertices)/(sizeof(double)*2);
static const size_t square_indices[4/*#segments*/*2/*#indices per segment*/]= {
0, 1, /* Bottom */
@@ -111,7 +111,7 @@ static const size_t square_indices[4/*#segments*/*2/*#indices per segment*/]= {
2, 3, /* Top */
3, 0 /* Right */
};
-static const size_t square_nsegments = sizeof(square_indices)/(2*sizeof(size_t));
+static const size_t square_nsegments = sizeof(square_indices)/(sizeof(size_t)*2);
static INLINE void
square_get_indices(const size_t iseg, size_t ids[2], void* context)
@@ -165,35 +165,36 @@ dummy_interface_getter
}
static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
- dummy_medium_getter,
- dummy_medium_getter,
- dummy_medium_getter,
- dummy_medium_getter,
- dummy_medium_getter,
- dummy_medium_getter,
- 0
+ dummy_medium_getter, /* Calorific capacity */
+ dummy_medium_getter, /* Thermal conductivity */
+ dummy_medium_getter, /* Volumic mass */
+ dummy_medium_getter, /* Delta */
+ dummy_medium_getter, /* Volumic power */
+ dummy_medium_getter, /* Temperature */
+ 0 /* Initial time */
};
static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
- dummy_medium_getter,
- dummy_medium_getter,
- dummy_medium_getter,
- 0
+ dummy_medium_getter, /* Calorific capacity */
+ dummy_medium_getter, /* Volumic mass */
+ dummy_medium_getter, /* Temperature */
+ 0 /* Initial time */
};
#define DUMMY_INTERFACE_SIDE_SHADER__ { \
- dummy_interface_getter, \
- dummy_interface_getter, \
- dummy_interface_getter, \
- dummy_interface_getter \
+ dummy_interface_getter, /* Temperature */ \
+ dummy_interface_getter, /* Flux */ \
+ dummy_interface_getter, /* Emissivity */ \
+ dummy_interface_getter, /* Specular fraction */ \
+ dummy_interface_getter /* Reference temperature */ \
}
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
- dummy_interface_getter,
- 0,
- dummy_interface_getter,
- DUMMY_INTERFACE_SIDE_SHADER__,
- DUMMY_INTERFACE_SIDE_SHADER__
+ dummy_interface_getter, /* Convection coef */
+ 0, /* Upper bound of the convection coef */
+ dummy_interface_getter, /* Thermal contact resistance */
+ DUMMY_INTERFACE_SIDE_SHADER__, /* Front side */
+ DUMMY_INTERFACE_SIDE_SHADER__ /* Back side */
};
/*******************************************************************************
diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c
@@ -66,7 +66,7 @@ static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = {
0.1, 0.6, 0.5,
0.1, 0.4, 0.5
};
-static const size_t nvertices = sizeof(vertices)/(3*sizeof(double));
+static const size_t nvertices = sizeof(vertices)/(sizeof(double)*3);
static const size_t indices[36/*#triangles*/*3/*#indices per triangle*/]= {
0, 4, 5, 5, 1, 0, /* Cuboid left */
@@ -90,7 +90,7 @@ static const size_t indices[36/*#triangles*/*3/*#indices per triangle*/]= {
8, 9, 10, 10, 11, 8, /* Cube back */
12, 15, 14, 14, 13, 12 /* Cube front */
};
-static const size_t ntriangles = sizeof(indices)/(3*sizeof(size_t));
+static const size_t ntriangles = sizeof(indices)/(sizeof(size_t)*3);
/*******************************************************************************
* Geometry
diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c
@@ -99,7 +99,7 @@ static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
0.1, 0.6,
0.1, 0.4
};
-static const size_t nvertices = sizeof(vertices)/(2*sizeof(double));
+static const size_t nvertices = sizeof(vertices)/(sizeof(double)*2);
static const size_t indices[8/*#segments*/*2/*#indices per segment*/]= {
0, 1, /* Rectangle left */
@@ -111,7 +111,7 @@ static const size_t indices[8/*#segments*/*2/*#indices per segment*/]= {
6, 7, /* Square right */
7, 4 /* Square bottom */
};
-static const size_t nsegments = sizeof(indices)/(2*sizeof(size_t));
+static const size_t nsegments = sizeof(indices)/(sizeof(size_t)*2);
/*******************************************************************************
* Geometry
diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c
@@ -84,7 +84,7 @@ static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
100000.5, 1.4, /* 6 */
100000.5, 0.0 /* 7 */
};
-static const size_t nvertices = sizeof(vertices)/(2*sizeof(double));
+static const size_t nvertices = sizeof(vertices)/(sizeof(double)*2);
static const size_t indices[10/*#segments*/*2/*#indices per segment*/]= {
0, 1,
@@ -98,7 +98,7 @@ static const size_t indices[10/*#segments*/*2/*#indices per segment*/]= {
6, 1,
2, 5
};
-static const size_t nsegments = sizeof(indices)/(2*sizeof(size_t));
+static const size_t nsegments = sizeof(indices)/(sizeof(size_t)*2);
/*******************************************************************************
* Geometry
