commit dca52ac7d74638d71ac388129a872166d83da91e
parent b23ec1cb998c44e4cc63ab453eba9ca9f72eae66
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 22 Jun 2018 14:48:05 +0200
Merge branch 'test_volumic_power' into develop
Diffstat:
26 files changed, 2547 insertions(+), 356 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -17,9 +17,14 @@ cmake_minimum_required(VERSION 3.0)
project(stardis C)
enable_testing()
+include(CMakeDependentOption)
+
set(SDIS_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
option(NO_TEST "Do not build tests" OFF)
+CMAKE_DEPENDENT_OPTION(ALL_TESTS
+ "Perform basic and advanced tests" OFF "NOT NO_TEST" OFF)
+
################################################################################
# Check dependencies
################################################################################
@@ -147,10 +152,23 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe_boundary)
new_test(test_sdis_volumic_power)
+ # Additionnal tests
+ build_test(test_sdis_volumic_power2)
+ build_test(test_sdis_volumic_power2_2d)
+ build_test(test_sdis_volumic_power3_2d)
+ build_test(test_sdis_volumic_power4_2d)
+
+ if(ALL_TESTS)
+ add_test(test_sdis_volumic_power2 test_sdis_volumic_power2)
+ add_test(test_sdis_volumic_power2_2d test_sdis_volumic_power2_2d)
+ add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
+ add_test(test_sdis_volumic_power4_2d test_sdis_volumic_power4_2d)
+ endif()
+
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
+ target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB})
target_link_libraries(test_sdis_solve_camera Star3DUT)
- target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB})
endif()
################################################################################
diff --git a/src/sdis.h b/src/sdis.h
@@ -17,6 +17,7 @@
#define SDIS_H
#include <rsys/rsys.h>
+#include <float.h>
/* Library symbol management */
#if defined(SDIS_SHARED_BUILD)
@@ -105,6 +106,7 @@ struct sdis_accum {
double sum_weights; /* Sum of Monte-Carlo weights */
double sum_weights_sqr; /* Sum of Monte-Carlo square weights */
size_t nweights; /* #accumulated weights */
+ size_t nfailures; /* #failures */
};
/* Monte-Carlo estimation */
@@ -137,7 +139,6 @@ struct sdis_solid_shader {
sdis_medium_getter_T thermal_conductivity; /* In W.m^-1.K^-1 */
sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */
sdis_medium_getter_T delta_solid;
- sdis_medium_getter_T delta_boundary;
/* May be NULL if there is no volumic power. One can also return
* SDIS_VOLUMIC_POWER_NONE to define that there is no volumic power at the
@@ -148,7 +149,7 @@ struct sdis_solid_shader {
* unknown for the submitted random walk vertex. */
sdis_medium_getter_T temperature;
};
-#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, NULL}
+#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL}
static const struct sdis_solid_shader SDIS_SOLID_SHADER_NULL =
SDIS_SOLID_SHADER_NULL__;
@@ -373,6 +374,10 @@ SDIS_API enum sdis_medium_type
sdis_medium_get_type
(const struct sdis_medium* medium);
+SDIS_API struct sdis_data*
+sdis_medium_get_data
+ (struct sdis_medium* medium);
+
/*******************************************************************************
* An interface is the boundary between 2 media.
******************************************************************************/
diff --git a/src/sdis_medium.c b/src/sdis_medium.c
@@ -39,7 +39,6 @@ check_solid_shader(const struct sdis_solid_shader* shader)
&& shader->thermal_conductivity
&& shader->volumic_mass
&& shader->delta_solid
- && shader->delta_boundary
&& shader->temperature;
}
@@ -205,3 +204,9 @@ sdis_medium_get_type(const struct sdis_medium* medium)
ASSERT(medium != NULL);
return medium->type;
}
+struct sdis_data*
+sdis_medium_get_data(struct sdis_medium* medium)
+{
+ ASSERT(medium);
+ return medium->data;
+}
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -94,14 +94,6 @@ solid_get_delta
}
static INLINE double
-solid_get_delta_boundary
- (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
-{
- ASSERT(mdm && mdm->type == SDIS_SOLID);
- return mdm->shader.solid.delta_boundary(vtx, mdm->data);
-}
-
-static INLINE double
solid_get_volumic_power
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -24,9 +24,6 @@
#include <rsys/double3.h>
#include <rsys/mem_allocator.h>
-#include <star/s2d.h>
-#include <star/s3d.h>
-
#include <limits.h>
/* Context used to wrap the user geometry to Star-XD. */
@@ -455,53 +452,76 @@ static INLINE res_T
scene_get_medium_2d
(const struct sdis_scene* scn,
const double pos[2],
+ struct get_medium_info* info, /* May be NULL */
const struct sdis_medium** out_medium)
{
const struct sdis_medium* medium = NULL;
size_t iprim, nprims;
size_t nfailures = 0;
const size_t max_failures = 10;
+ /* Range of the parametric coordinate into which positions are challenged */
+ const float s_range[2] = {0.25, 0.75};
+ const size_t s_nsteps = 3; /* #challenges per primitive into the range */
+ float s_step;
res_T res = RES_OK;
ASSERT(scn && pos);
+ s_step = (s_range[1] - s_range[0]) / (float)(s_nsteps-1);
+
S2D(scene_view_primitives_count(scn->s2d_view, &nprims));
FOR_EACH(iprim, 0, nprims) {
struct s2d_hit hit;
struct s2d_attrib attr;
struct s2d_primitive prim;
- float s;
const float range[2] = {0.f, FLT_MAX};
float N[2], P[2], dir[2], cos_N_dir;
- s = 1.f / 3.f;
-
- /* Retrieve a position onto the primitive */
- S2D(scene_view_get_primitive(scn->s2d_view, (unsigned)iprim, &prim));
- S2D(primitive_get_attrib(&prim, S2D_POSITION, s, &attr));
+ float s = s_range[0];
+
+ do {
+ /* Retrieve a position onto the primitive */
+ S2D(scene_view_get_primitive(scn->s2d_view, (unsigned)iprim, &prim));
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, s, &attr));
+
+ /* Trace a ray from the random walk vertex toward the retrieved primitive
+ * position */
+ f2_normalize(dir, f2_sub(dir, attr.value, f2_set_d2(P, pos)));
+ S2D(scene_view_trace_ray(scn->s2d_view, P, dir, range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(S2D_HIT_NONE(&hit)) {
+ ++nfailures;
+ if(nfailures < max_failures) {
+ continue;
+ } else {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Discard the hit if it is on a vertex, i.e. between 2 segments, and
+ * target a new position onto the current primitive */
+ } while((S2D_HIT_NONE(&hit) || hit_on_vertex(&hit))
+ && (s+=s_step) <= s_range[1]);
- /* Trace a ray from the randomw walk vertex toward the retrieved primitive
- * position */
- f2_normalize(dir, f2_sub(dir, attr.value, f2_set_d2(P, pos)));
- S2D(scene_view_trace_ray(scn->s2d_view, P, dir, range, NULL, &hit));
+ /* The hits of all targeted positions on the current primitive are on
+ * vertices. Challenge positions on another primitive. */
+ if(s > s_range[1]) continue;
f2_normalize(N, hit.normal);
cos_N_dir = f2_dot(N, dir);
- /* Unforeseen error. One has to intersect a primitive ! */
- if(S2D_HIT_NONE(&hit)) {
- ++nfailures;
- if(nfailures < max_failures) {
- continue;
- } else {
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
const struct sdis_interface* interf;
interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
(interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+
+ /* Register the get_medium_info */
+ if(info) {
+ f2_set(info->pos_tgt, attr.value);
+ f2_set(info->ray_org, P);
+ f2_set(info->ray_dir, dir);
+ info->hit_2d = hit;
+ }
break;
}
}
@@ -519,53 +539,75 @@ static INLINE res_T
scene_get_medium_3d
(const struct sdis_scene* scn,
const double pos[3],
+ struct get_medium_info* info,
const struct sdis_medium** out_medium)
{
const struct sdis_medium* medium = NULL;
size_t iprim, nprims;
size_t nfailures = 0;
const size_t max_failures = 10;
+ float st[3][2]; /* Position to challenge onto the primitive */
+ const size_t nsteps = sizeof(st)/sizeof(float[2]);
res_T res = RES_OK;
ASSERT(scn && pos);
+ /* Setup the position to challenge onto the primitives */
+ f2(st[0], 1.f/6.f, 5.f/12.f);
+ f2(st[1], 5.f/12.f, 1.f/6.f);
+ f2(st[2], 5.f/12.f, 5.f/12.f);
+
S3D(scene_view_primitives_count(scn->s3d_view, &nprims));
FOR_EACH(iprim, 0, nprims) {
struct s3d_hit hit;
struct s3d_attrib attr;
struct s3d_primitive prim;
- float st[2];
const float range[2] = {0.f, FLT_MAX};
float N[3], P[3], dir[3], cos_N_dir;
- st[0] = st[1] = 1.f / 3.f; /* Or MSVC will issue a warning */
-
- /* Retrieve a position onto the primitive */
- S3D(scene_view_get_primitive(scn->s3d_view, (unsigned)iprim, &prim));
- S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attr));
+ size_t istep = 0;
+
+ do {
+ /* Retrieve a position onto the primitive */
+ S3D(scene_view_get_primitive(scn->s3d_view, (unsigned)iprim, &prim));
+ S3D(primitive_get_attrib(&prim, S3D_POSITION, st[istep], &attr));
+
+ /* Trace a ray from the randomw walk vertex toward the retrieved primitive
+ * position */
+ f3_normalize(dir, f3_sub(dir, attr.value, f3_set_d3(P, pos)));
+ S3D(scene_view_trace_ray(scn->s3d_view, P, dir, range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(S3D_HIT_NONE(&hit)) {
+ ++nfailures;
+ if(nfailures < max_failures) {
+ continue;
+ } else {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Discard the hit if it is on an edge, i.e. between 2 triangles, and
+ * target a new position onto the current primitive */
+ } while((S3D_HIT_NONE(&hit) || hit_on_edge(&hit)) && ++istep < nsteps);
- /* Trace a ray from the randomw walk vertex toward the retrieved primitive
- * position */
- f3_normalize(dir, f3_sub(dir, attr.value, f3_set_d3(P, pos)));
- S3D(scene_view_trace_ray(scn->s3d_view, P, dir, range, NULL, &hit));
+ /* The hits of all targeted positions on the current primitive are on
+ * edges. Challenge positions on another primitive. */
+ if(istep >= nsteps) continue;
f3_normalize(N, hit.normal);
cos_N_dir = f3_dot(N, dir);
- /* Unforeseen error. One has to intersect a primitive ! */
- if(S3D_HIT_NONE(&hit)) {
- ++nfailures;
- if(nfailures < max_failures) {
- continue;
- } else {
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
const struct sdis_interface* interf;
interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
(interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+
+ if(info) {
+ f3_set(info->pos_tgt, attr.value);
+ f3_set(info->ray_org, P);
+ f3_set(info->ray_dir, dir);
+ info->hit_3d = hit;
+ }
break;
}
}
@@ -775,10 +817,11 @@ res_T
scene_get_medium
(const struct sdis_scene* scn,
const double pos[],
+ struct get_medium_info* info,
const struct sdis_medium** out_medium)
{
return scene_is_2d(scn)
- ? scene_get_medium_2d(scn, pos, out_medium)
- : scene_get_medium_3d(scn, pos, out_medium);
+ ? scene_get_medium_2d(scn, pos, info, out_medium)
+ : scene_get_medium_3d(scn, pos, info, out_medium);
}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -19,6 +19,20 @@
#include <rsys/dynamic_array.h>
#include <rsys/ref_count.h>
+#include <star/s2d.h>
+#include <star/s3d.h>
+
+struct get_medium_info {
+ /* Targeted position */
+ float pos_tgt[3];
+ /* Ray trace to the targeted position in order to define the current medium */
+ float ray_org[3];
+ float ray_dir[3];
+ /* Hit encouters along the ray and used to define the current medium */
+ struct s2d_hit hit_2d;
+ struct s3d_hit hit_3d;
+};
+
static INLINE void
interface_init
(struct mem_allocator* allocator,
@@ -62,6 +76,7 @@ extern LOCAL_SYM res_T
scene_get_medium
(const struct sdis_scene* scene,
const double position[],
+ struct get_medium_info* info, /* May be NULL */
const struct sdis_medium** medium);
static FINLINE int
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -98,6 +98,7 @@ solve_pixel
accum->sum_weights = sum_weights;
accum->sum_weights_sqr = sum_weights_sqr;
accum->nweights = N;
+ accum->nfailures = nrealisations - N;
exit:
return res;
@@ -210,7 +211,7 @@ sdis_solve_probe
if(res != RES_OK) goto error;
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, position, &medium);
+ res = scene_get_medium(scn, position, NULL, &medium);
if(res != RES_OK) goto error;
/* Here we go! Launch the Monte Carlo estimation */
@@ -218,7 +219,7 @@ sdis_solve_probe
#pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
for(irealisation = 0; irealisation < rcount; ++irealisation) {
res_T res_local;
- double w;
+ double w = NaN;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -361,7 +362,7 @@ sdis_solve_probe_boundary
#pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
for(irealisation = 0; irealisation < rcount; ++irealisation) {
res_T res_local;
- double w;
+ double w = NaN;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -451,7 +452,7 @@ sdis_solve_camera
}
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, cam->position, &medium);
+ res = scene_get_medium(scn, cam->position, NULL, &medium);
if(res != RES_OK) goto error;
if(medium->type != SDIS_FLUID) {
@@ -496,7 +497,7 @@ sdis_solve_camera
pix_sz[1] = 1.0 / (double)height;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static, 1/*chunki size*/)
+ #pragma omp parallel for schedule(static, 1/*chunk size*/)
for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
size_t tile_org[2] = {0, 0};
size_t tile_sz[2] = {0, 0};
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -22,13 +22,26 @@
#include "sdis_medium_c.h"
#include "sdis_scene_c.h"
+#include <rsys/float2.h>
+#include <rsys/float3.h>
#include <rsys/stretchy_array.h>
#include <star/ssp.h>
-/* Emperical scale factor to apply to the upper bound of the ray range in order
+/* Define a new result code from RES_BAD_OP saying that the bad operation is
+ * definitive, i.e. in the current state, the realisation will inevitably fail.
+ * It is thus unecessary to retry a specific section of the random walk */
+#define RES_BAD_OP_IRRECOVERABLE (-RES_BAD_OP)
+
+/* Empirical scale factor to apply to the upper bound of the ray range in order
* to handle numerical imprecisions */
-#define RAY_RANGE_MAX_SCALE 1.0001f
+#define RAY_RANGE_MAX_SCALE 1.001f
+
+/* Emperical scale factor applied to the challenged reinjection distance. If
+ * the distance to reinject is less than this adjusted value, the solver
+ * switches from 2D reinjection scheme to the 1D reinjection scheme in order to
+ * avoid numerical issues. */
+#define REINJECT_DST_MIN_SCALE 0.125f
#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
@@ -40,7 +53,22 @@ struct rwalk_context {
/* Reflect the vector V wrt the normal N. By convention V points outward the
* surface. */
static INLINE float*
-reflect(float res[3], const float V[3], const float N[3])
+reflect_2d(float res[2], const float V[2], const float N[2])
+{
+ float tmp[2];
+ float cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(f2_is_normalized(V) && f2_is_normalized(N));
+ cos_V_N = f2_dot(V, N);
+ f2_mulf(tmp, N, 2*cos_V_N);
+ f2_sub(res, tmp, V);
+ return res;
+}
+
+/* Reflect the vector V wrt the normal N. By convention V points outward the
+ * surface. */
+static INLINE float*
+reflect_3d(float res[3], const float V[3], const float N[3])
{
float tmp[3];
float cos_V_N;
@@ -163,6 +191,48 @@ XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
#endif
}
+static FINLINE 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
+}
+
/* Check that the interface fragment is consistent with the current state of
* the random walk */
static INLINE int
@@ -248,7 +318,7 @@ XD(trace_radiative_path)
FUNC_NAME,
ctx->Tarad,
SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
+ res = RES_BAD_OP;
goto error;
}
}
@@ -298,10 +368,10 @@ XD(trace_radiative_path)
res = RES_BAD_OP;
goto error;
}
- alpha = interface_side_get_specular_fraction(interf, &frag);
+ alpha = interface_side_get_specular_fraction(interf, &frag);
r = ssp_rng_canonical(rng);
if(r < alpha) { /* Sample specular part */
- reflect(dir, f3_minus(dir, dir), N);
+ reflect_3d(dir, f3_minus(dir, dir), N);
} else { /* Sample diffuse part */
ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
}
@@ -387,17 +457,27 @@ XD(solid_solid_boundary_temperature)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
+ struct sXd(hit) hit0, hit1, hit2, hit3;
+ struct sXd(hit)* hit;
const struct sdis_interface* interf = NULL;
const struct sdis_medium* solid_front = NULL;
const struct sdis_medium* solid_back = NULL;
+ const struct sdis_medium* mdm;
double lambda_front, lambda_back;
- double delta_front_boundary, delta_back_boundary;
- double delta_front_boundary_meter, delta_back_boundary_meter;
+ double delta_front, delta_back;
+ double delta_boundary_front, delta_boundary_back;
double delta_boundary;
+ double reinject_dst_front, reinject_dst_back;
+ double reinject_dst;
double proba;
double tmp;
double r;
- float pos[DIM], dir[DIM], range[2];
+ double power;
+ float range0[2], range1[2];
+ float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
+ float* dir;
+ float pos[DIM];
+ int dim = DIM;
ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
(void)frag, (void)ctx;
@@ -412,40 +492,105 @@ XD(solid_solid_boundary_temperature)
/* 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);
- delta_front_boundary = solid_get_delta_boundary(solid_front, &rwalk->vtx);
- delta_back_boundary = solid_get_delta_boundary(solid_back, &rwalk->vtx);
-
- /* Convert the delta boundary from floating point units to meters */
- delta_front_boundary_meter = fp_to_meter * delta_front_boundary;
- delta_back_boundary_meter = fp_to_meter * delta_back_boundary;
-
- /* Compute the proba to switch in solid0 or in solid1 */
- tmp = lambda_front / delta_front_boundary_meter;
- proba = tmp / (tmp + lambda_back / delta_back_boundary_meter);
+ /* 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);
+
+ /* 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);
+
+ /* Trace the sampled directions on both sides of the interface to adjust the
+ * reinjection distance of the random walk . */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range0, 0, (float)delta_boundary_front*RAY_RANGE_MAX_SCALE);
+ f2(range1, 0, (float)delta_boundary_back *RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range0, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range1, &rwalk->hit, &hit1));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir2, range0, &rwalk->hit, &hit2));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir3, range1, &rwalk->hit, &hit3));
+
+ /* Adjust the reinjection distance */
+ reinject_dst_front = MMIN(MMIN(delta_boundary_front, hit0.distance), hit2.distance);
+ reinject_dst_back = MMIN(MMIN(delta_boundary_back, hit1.distance), hit3.distance);
+
+ /* 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. */
r = ssp_rng_canonical(rng);
- fX(normalize)(dir, rwalk->hit.normal);
+ proba = (lambda_front/reinject_dst_front)
+ / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
if(r < proba) { /* Reinject in front */
- delta_boundary = delta_front_boundary;
- rwalk->mdm = solid_front;
+ dir = dir0;
+ hit = &hit0;
+ mdm = solid_front;
+ reinject_dst = reinject_dst_front;
+ delta_boundary = delta_boundary_front;
} else { /* Reinject in back */
- delta_boundary = delta_back_boundary;
- fX(minus)(dir, dir);
- rwalk->mdm = solid_back;
+ dir = dir1;
+ hit = &hit1;
+ mdm = solid_back;
+ reinject_dst = reinject_dst_back;
+ delta_boundary = delta_boundary_back;
}
- /* "Reinject" the path into the solid along the surface normal. */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
+ /* Switch in 1D reinjection scheme */
+ if(reinject_dst < delta_boundary * REINJECT_DST_MIN_SCALE) {
+ if(dir == dir0) {
+ fX(set)(dir, rwalk->hit.normal);
+ } else {
+ fX(minus)(dir, rwalk->hit.normal);
+ }
+
+ f2(range0, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir, range0, &rwalk->hit, hit));
+ reinject_dst = MMIN(delta_boundary, hit->distance),
+ dim = 1;
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance */
+ if(!SXD_HIT_NONE(hit)) {
+ reinject_dst *= 0.5;
+ *hit = SXD_HIT_NULL;
+ }
+ }
+
+ /* 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 * fp_to_meter;
+ const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
+
+ /* Reinject */
+ XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
+ if(eq_epsf(hit->distance, (float)reinject_dst, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = *hit;
+ rwalk->hit_side = fX(dot)(hit->normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
+ T->func = XD(solid_temperature);
+ rwalk->mdm = mdm;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ }
}
static void
@@ -463,6 +608,8 @@ XD(solid_fluid_boundary_temperature)
const struct sdis_medium* mdm_back = NULL;
const struct sdis_medium* solid = NULL;
const struct sdis_medium* fluid = NULL;
+ struct sXd(hit) hit0 = SXD_HIT_NULL;
+ struct sXd(hit) hit1 = SXD_HIT_NULL;
struct sdis_interface_fragment frag_fluid;
double hc;
double hr;
@@ -470,15 +617,19 @@ XD(solid_fluid_boundary_temperature)
double lambda;
double fluid_proba;
double radia_proba;
+ double delta;
double delta_boundary;
double r;
double tmp;
- float dir[DIM], pos[DIM], range[2];
+ float pos[DIM];
+ float dir0[DIM], dir1[DIM];
+ float range[2];
+ int dim = DIM;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- /* Retrieve the solid and the fluid split by the boundary */
+ /* 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);
@@ -497,7 +648,54 @@ XD(solid_fluid_boundary_temperature)
/* Fetch the solid properties */
lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta_boundary = solid_get_delta_boundary(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;
+
+ /* 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);
+ }
+
+ /* Trace dir0/dir1 to adjust the reinjection distance */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
+
+ /* Adjust the delta boundary to the hit distance */
+ tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
+
+ if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
+ delta_boundary = tmp;
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = delta_boundary / sqrt(DIM);
+ } else { /* Switch in 1D reinjection scheme. */
+ fX(set)(dir0, rwalk->hit.normal);
+ if(solid == mdm_back) fX(minus)(dir0, dir0);
+ f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ delta_boundary = MMIN(hit0.distance, delta);
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance in order to avoid infinite bounces for parallel plane */
+ if(!SXD_HIT_NONE(&hit0)) {
+ delta_boundary *= 0.5;
+ hit0 = SXD_HIT_NULL;
+ }
+
+ delta = delta_boundary;
+ dim = 1;
+ }
/* Fetch the boundary properties */
epsilon = interface_side_get_emissivity(interf, &frag_fluid);
@@ -506,8 +704,8 @@ XD(solid_fluid_boundary_temperature)
/* Compute the radiative coefficient */
hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
- /* Compute the probas to switch in solid or fluid random walk */
- tmp = lambda / (delta_boundary*fp_to_meter);
+ /* Compute the probas to switch in solid, fluid or radiative random walk */
+ tmp = lambda / (delta*fp_to_meter);
fluid_proba = hc / (tmp + hr + hc);
radia_proba = hr / (tmp + hr + hc);
/*solid_proba = tmp / (tmp + hr + hc);*/
@@ -522,20 +720,140 @@ XD(solid_fluid_boundary_temperature)
rwalk->mdm = fluid;
rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
} else { /* Solid random walk */
- rwalk->mdm = solid;
- fX(normalize)(dir, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir, dir);
+ /* 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 = delta_boundary * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
- /* "Reinject" the random walk into the solid */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+ /* Reinject */
+ XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
+ if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
+ T->func = XD(solid_temperature);
+ rwalk->mdm = solid;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ }
+ }
+}
- /* Switch in solid random walk */
+static void
+XD(solid_boundary_with_flux_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ 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)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm = NULL;
+ double lambda;
+ double delta;
+ double delta_boundary;
+ double delta_in_meter;
+ double power;
+ double tmp;
+ struct sXd(hit) hit0;
+ struct sXd(hit) hit1;
+ float pos[DIM];
+ float dir0[DIM];
+ float dir1[DIM];
+ float range[2];
+ int dim = DIM;
+ 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);
+
+ /* 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);
+ }
+
+ /* Trace dir0/dir1 to adjust the reinjection distance wrt the geometry */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
+
+ /* Adjust the delta boundary to the hit distance */
+ tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
+
+ if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
+ delta_boundary = tmp;
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = delta_boundary / sqrt(DIM);
+ } else { /* Switch in 1D reinjection scheme. */
+ fX(set)(dir0, rwalk->hit.normal);
+ if(frag->side == SDIS_BACK) fX(minus)(dir0, dir0);
+ f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ delta_boundary = MMIN(hit0.distance, delta_boundary);
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance in order to avoid infinite bounces for parallel plane */
+ if(!SXD_HIT_NONE(&hit0)) {
+ delta_boundary *= 0.5;
+ hit0 = SXD_HIT_NULL;
+ }
+
+ delta = delta_boundary;
+ dim = 1;
+ }
+
+ /* Handle the flux */
+ delta_in_meter = delta*fp_to_meter;
+ T->value += phi * delta_in_meter / lambda;
+
+ /* Handle the volumic power */
+ power = solid_get_volumic_power(mdm, &rwalk->vtx);
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ delta_in_meter = delta_boundary * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
+
+ /* Reinject into the solid */
+ XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
+ if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
T->func = XD(solid_temperature);
+ rwalk->mdm = mdm;
rwalk->hit = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
}
@@ -562,6 +880,8 @@ XD(boundary_temperature)
XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+ fX(normalize)(rwalk->hit.normal, rwalk->hit.normal);
+
/* Retrieve the current interface */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
@@ -576,40 +896,15 @@ XD(boundary_temperature)
/* Check if the boundary flux is known. Note that actually, only solid media
* can have a flux as limit condition */
mdm = interface_get_medium(interf, frag.side);
- if(sdis_medium_get_type(mdm) == SDIS_SOLID) {
+ if(sdis_medium_get_type(mdm) == SDIS_SOLID ) {
const double phi = interface_side_get_flux(interf, &frag);
-
if(phi != SDIS_FLUX_NONE) {
- double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- double delta_b = solid_get_delta_boundary(mdm, &rwalk->vtx);
- double delta_b_in_meter = delta_b * fp_to_meter;
- float pos[3];
- float range[2];
- float dir[3];
-
- /* Update the temperature */
- T->value += phi * delta_b_in_meter / lambda;
-
- /* Ensuure that the normal points toward the solid */
- fX(normalize)(dir, rwalk->hit.normal);
- if(frag.side == SDIS_BACK) fX(minus)(dir, dir);
-
- /* "Reinject" the random walk into the solid */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_b*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_b = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_b);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
- rwalk->hit = SXD_HIT_NULL;
- rwalk->hit_side = SDIS_SIDE_NULL__;
- rwalk->mdm = mdm;
+ XD(solid_boundary_with_flux_temperature)
+ (scn, fp_to_meter, ctx, &frag, phi, rwalk, rng, T);
return RES_OK;
}
}
+
mdm_front = interface_get_medium(interf, SDIS_FRONT);
mdm_back = interface_get_medium(interf, SDIS_BACK);
@@ -639,17 +934,18 @@ XD(solid_temperature)
(void)ctx;
/* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK);
if(mdm != rwalk->mdm) {
log_err(scn->dev, "%s: invalid solid random walk. "
"Unexpected medium at {%g, %g, %g}.\n",
FUNC_NAME, SPLIT3(rwalk->vtx.P));
- return RES_BAD_OP;
+ return RES_BAD_OP_IRRECOVERABLE;
}
/* Save the submitted position */
dX(set)(position_start, rwalk->vtx.P);
do { /* Solid random walk */
+ struct get_medium_info info;
struct sXd(hit) hit0, hit1;
double lambda; /* Thermal conductivity */
double rho; /* Volumic mass */
@@ -675,6 +971,7 @@ XD(solid_temperature)
lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+ power = solid_get_volumic_power(mdm, &rwalk->vtx);
#if (SDIS_SOLVE_DIMENSION == 2)
/* Sample a direction around 2PI */
@@ -696,9 +993,57 @@ XD(solid_temperature)
if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
/* Hit nothing: move along dir0 of the original delta */
delta = delta_solid;
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_in_meter = delta * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
+ T->value += tmp;
+ }
} else {
/* Hit something: move along dir0 of the minimum hit distance */
delta = MMIN(hit0.distance, hit1.distance);
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_s_in_meter = delta_solid * fp_to_meter;
+ double h;
+ double h_in_meter;
+ double cos_U_N;
+ float N[DIM];
+
+ if(delta == hit0.distance) {
+ fX(normalize)(N, hit0.normal);
+ cos_U_N = fX(dot)(dir0, N);
+ } else {
+ ASSERT(delta == hit1.distance);
+ fX(normalize)(N, hit1.normal);
+ cos_U_N = fX(dot)(dir1, N);
+ }
+
+ h = delta * fabs(cos_U_N);
+ h_in_meter = h * fp_to_meter;
+
+ /* The regular power term at wall */
+ tmp = power * h_in_meter * h_in_meter / (2.0 * lambda);
+
+ /* Add the power corrective term */
+ if(h < delta_solid) {
+ const double sin_a = h / delta_solid;
+#if DIM==2
+ /* tmp1 = sin(2a) / (PI - 2*a) */
+ const double tmp1 = sin_a * sqrt(1 - sin_a*sin_a)/acos(sin_a);
+ tmp += -(power*delta_s_in_meter*delta_s_in_meter)/(4.0*lambda) * tmp1;
+#else
+ const double tmp1 = (sin_a*sin_a*sin_a - sin_a)/ (1-sin_a);
+ tmp += (power*delta_s_in_meter*delta_s_in_meter)/(6*lambda) * tmp1;
+#endif
+
+ } else if (h == delta_solid) {
+ tmp += -(delta_s_in_meter*delta_s_in_meter*power)/(2.0*DIM*lambda);
+ }
+ T->value += tmp;
+ }
}
/* Sample the time */
@@ -714,14 +1059,6 @@ XD(solid_temperature)
return RES_OK;
}
- /* Add the volumic power density to the measured temperature */
- power = solid_get_volumic_power(mdm, &rwalk->vtx);
- if(power != SDIS_VOLUMIC_POWER_NONE) {
- const double delta_in_meter = delta * fp_to_meter;
- tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += tmp;
- }
-
/* Define if the random walk hits something along dir0 */
if(hit0.distance > delta) {
rwalk->hit = SXD_HIT_NULL;
@@ -736,7 +1073,7 @@ XD(solid_temperature)
/* Fetch the current medium */
if(SXD_HIT_NONE(&rwalk->hit)) {
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK);
} else {
const struct sdis_interface* interf;
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
@@ -747,15 +1084,28 @@ XD(solid_temperature)
if(mdm != rwalk->mdm) {
log_err(scn->dev,
"%s: inconsistent medium during the solid random walk.\n", FUNC_NAME);
- if(DIM == 2) {
- log_err(scn->dev,
- " start position: %g %g; current position: %g %g\n",
- SPLIT2(position_start), SPLIT2(rwalk->vtx.P));
- } else {
- log_err(scn->dev,
- " start position: %g %g %g; current position: %g %g %g\n",
- SPLIT3(position_start), SPLIT3(rwalk->vtx.P));
+#if DIM == 2
+ #define VEC_STR "%g %g"
+ #define VEC_SPLIT SPLIT2
+#else
+ #define VEC_STR "%g %g %g"
+ #define VEC_SPLIT SPLIT3
+#endif
+ log_err(scn->dev,
+ " start position: " VEC_STR "; current position: " VEC_STR "\n",
+ VEC_SPLIT(position_start), VEC_SPLIT(rwalk->vtx.P));
+ if(SXD_HIT_NONE(&rwalk->hit)) {
+ float hit_pos[DIM];
+ fX(mulf)(hit_pos, info.ray_dir, info.XD(hit).distance);
+ fX(add)(hit_pos, info.ray_org, hit_pos);
+ log_err(scn->dev, " ray org: " VEC_STR "; ray dir: " VEC_STR "\n",
+ VEC_SPLIT(info.ray_org), VEC_SPLIT(info.ray_dir));
+ log_err(scn->dev, " targeted point: " VEC_STR "\n",
+ VEC_SPLIT(info.pos_tgt));
+ log_err(scn->dev, " hit pos: " VEC_STR "\n", VEC_SPLIT(hit_pos));
}
+#undef VEC_STR
+#undef VEC_SPLIT
return RES_BAD_OP;
}
@@ -777,18 +1127,37 @@ XD(compute_temperature)
struct XD(temperature)* T)
{
#ifndef NDEBUG
- struct XD(temperature)* stack = NULL;
+ struct entry {
+ struct XD(temperature) temperature;
+ struct XD(rwalk) rwalk;
+ }* stack = NULL;
size_t istack = 0;
#endif
+ const size_t max_fails = 10;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
do {
+ /* Save the current random walk state */
+ const struct XD(rwalk) rwalk_bkp = *rwalk;
+ const struct XD(temperature) T_bkp = *T;
+
+ size_t nfails = 0;
+
#ifndef NDEBUG
- sa_push(stack, *T);
+ struct entry e;
+ e.temperature = *T;
+ e.rwalk = *rwalk;
+ sa_push(stack, e);
++istack;
#endif
- res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
+
+ /* Reject the current step if a BAD_OP occurs and retry up to "max_fails"
+ * times */
+ do {
+ res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
+ if(res == RES_BAD_OP) { *rwalk = rwalk_bkp; *T = T_bkp; }
+ } while(res == RES_BAD_OP && ++nfails < max_fails);
if(res != RES_OK) goto error;
} while(!T->done);
@@ -797,7 +1166,7 @@ exit:
#ifndef NDEBUG
sa_release(stack);
#endif
- return res;
+ return res == RES_BAD_OP_IRRECOVERABLE ? RES_BAD_OP : res;
error:
goto exit;
}
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -168,14 +168,6 @@ solid_get_delta
return 1.0/10.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx != NULL); (void)data;
- return 2.1/10.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -303,7 +295,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -316,7 +307,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -399,12 +389,13 @@ main(int argc, char** argv)
double ref, u;
size_t nreals = 0;
size_t nfails = 0;
+ const size_t N = 10000;
pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[2] = ssp_rng_uniform_double(rng, -0.9, 0.9);
- CHK(sdis_solve_probe(scn, 10000, pos, INF, 1, -1, Tref, &estimator) == RES_OK);
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1, -1, Tref, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
@@ -413,7 +404,10 @@ main(int argc, char** argv)
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -149,14 +149,6 @@ solid_get_delta
return 1.0/10.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx != NULL); (void)data;
- return 2.1/10.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -317,7 +309,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -330,7 +321,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -403,6 +393,7 @@ main(int argc, char** argv)
double ref, u;
size_t nreals = 0;
size_t nfails = 0;
+ const size_t N = 10000;
pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
@@ -416,8 +407,11 @@ main(int argc, char** argv)
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE) == 1);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
}
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -150,15 +150,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -235,7 +226,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -310,25 +300,27 @@ main(int argc, char** argv)
CHK(sdis_solve_probe(box_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*3));
/* Solve in 2D */
CHK(sdis_solve_probe(square_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2.0));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*3));
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
CHK(sdis_scene_ref_put(square_scn) == RES_OK);
diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c
@@ -20,6 +20,7 @@ int
main(int argc, char** argv)
{
struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
@@ -40,6 +41,9 @@ main(int argc, char** argv)
CHK(sdis_fluid_create(NULL, &fluid_shader, NULL, &fluid) == RES_BAD_ARG);
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+ CHK(sdis_medium_get_type(fluid) == SDIS_FLUID);
+ CHK(sdis_medium_get_data(fluid) == NULL);
+
CHK(sdis_medium_ref_get(NULL) == RES_BAD_ARG);
CHK(sdis_medium_ref_get(fluid) == RES_OK);
CHK(sdis_medium_ref_put(NULL) == RES_BAD_ARG);
@@ -61,15 +65,19 @@ main(int argc, char** argv)
CHK(sdis_fluid_create
(dev, &SDIS_FLUID_SHADER_NULL, NULL, &fluid) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, NULL, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, NULL, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, &solid_shader, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, NULL, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, NULL, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, &solid_shader, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+ CHK(sdis_data_create(dev, 4, 16, NULL, &data) == RES_OK);
+ CHK(sdis_solid_create(NULL, NULL, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, NULL, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, &solid_shader, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, &solid_shader, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, NULL, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, NULL, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, &solid_shader, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_medium_get_type(solid) == SDIS_SOLID);
+ CHK(sdis_medium_get_data(solid) == data);
CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
solid_shader.calorific_capacity = NULL;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
@@ -87,10 +95,6 @@ main(int argc, char** argv)
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
solid_shader.delta_solid = DUMMY_SOLID_SHADER.delta_solid;
- solid_shader.delta_boundary = NULL;
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
- solid_shader.delta_boundary = DUMMY_SOLID_SHADER.delta_boundary;
-
solid_shader.temperature = NULL;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
solid_shader.temperature = DUMMY_SOLID_SHADER.temperature;
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -30,7 +30,7 @@
#define UNKOWN_TEMPERATURE -1
#define IMG_WIDTH 640
#define IMG_HEIGHT 480
-#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */
+#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */
/*
* The scene is composed of a solid cube whose temperature is unknown. The
@@ -222,14 +222,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->delta * 2.1;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -308,7 +300,6 @@ create_solid
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
/* Create the solid medium */
@@ -451,7 +442,7 @@ dump_image(const struct sdis_accum_buffer* buf)
double Tmax = -DBL_MAX;
double Tmin = DBL_MAX;
double norm;
- size_t ix, iy;
+ size_t i, ix, iy;
CHK(buf != NULL);
CHK(sdis_accum_buffer_get_layout(buf, &layout) == RES_OK);
@@ -459,8 +450,17 @@ dump_image(const struct sdis_accum_buffer* buf)
temps = mem_alloc(layout.width*layout.height*sizeof(double));
CHK(temps != NULL);
- /* Compute the per pixel temperature */
CHK(sdis_accum_buffer_map(buf, &accums) == RES_OK);
+
+ /* Check the results validity */
+ FOR_EACH(i, 0, layout.height * layout.width) {
+ CHK(accums[i].nweights + accums[i].nfailures == SPP);
+ CHK(accums[i].nfailures <= SPP/100);
+ CHK(accums[i].sum_weights >= 0);
+ CHK(accums[i].sum_weights_sqr >= 0);
+ }
+
+ /* Compute the per pixel temperature */
FOR_EACH(iy, 0, layout.height) {
const struct sdis_accum* row_accums = accums + iy * layout.width;
double* row = temps + iy * layout.width;
@@ -605,6 +605,11 @@ main(int argc, char** argv)
/* Create the accum buffer */
CHK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf) == RES_OK);
+#if 0
+ dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
+ exit(0);
+#endif
+
/* Launch the simulation */
CHK(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
sdis_accum_buffer_write, buf) == RES_OK);
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -126,14 +126,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->delta * 2.1;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -230,7 +222,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -285,8 +276,6 @@ main(int argc, char** argv)
CHK(sdis_estimator_get_failure_count(NULL, &nfails) == RES_BAD_ARG);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, NULL) == RES_BAD_ARG);
CHK(sdis_estimator_get_temperature(NULL, &T) == RES_BAD_ARG);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
@@ -294,6 +283,10 @@ main(int argc, char** argv)
ref = 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, T.SE));
CHK(sdis_estimator_ref_get(NULL) == RES_BAD_ARG);
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -114,15 +114,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -180,7 +171,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -191,7 +182,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -260,12 +250,12 @@ main(int argc, char** argv)
ref = 350 * pos[2] + (1-pos[2]) * 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Release data */
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -111,15 +111,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -177,7 +168,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -188,7 +179,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -258,12 +248,12 @@ main(int argc, char** argv)
ref = 350 * pos[0] + (1-pos[0]) * 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*2));
/* Release data */
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -136,15 +136,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -207,7 +198,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -218,7 +209,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -316,11 +306,11 @@ main(int argc, char** argv)
ref = 350 * pos[2] + (1-pos[2]) * 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 2*T.SE));
/* Release data */
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -133,15 +133,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -202,7 +193,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -213,7 +204,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -309,11 +299,11 @@ main(int argc, char** argv)
ref = 350 * pos[0] + (1-pos[0]) * 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 3*T.SE));
/* Release data */
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -108,14 +108,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)vtx, (void)data;
- return 2.1/20.0;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -159,7 +151,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
@@ -170,7 +162,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -202,13 +193,15 @@ main(int argc, char** argv)
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
ref = 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, T.SE));
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c
@@ -161,15 +161,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -238,7 +229,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
@@ -249,7 +240,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
@@ -335,38 +325,34 @@ main(int argc, char** argv)
CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG);
CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
+ ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
+
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
-
CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
-
- ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
-
printf("Boundary temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
uv[0] = 0.5;
iprim = 3;
CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
-
CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
-
printf("Boundary temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
#undef SOLVE
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -69,7 +69,7 @@ static const double square_vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
};
static const size_t square_nvertices = sizeof(square_vertices)/sizeof(double[2]);
-static const size_t square_indices[4/*#triangles*/*2/*#indices per segment*/]= {
+static const size_t square_indices[4/*#segments*/*2/*#indices per segment*/]= {
0, 1, /* Bottom */
1, 2, /* Left */
2, 3, /* Top */
@@ -104,7 +104,6 @@ static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
dummy_medium_getter,
dummy_medium_getter,
dummy_medium_getter,
- dummy_medium_getter,
dummy_medium_getter
};
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -49,6 +49,7 @@
#define T0 320
#define LAMBDA 0.1
#define P0 10
+#define DELTA 1.0/20.0
/*******************************************************************************
* Geometry 3D
@@ -114,6 +115,13 @@ square_get_interface
/*******************************************************************************
* Media
******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double delta;
+};
+
static double
fluid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
@@ -127,45 +135,32 @@ static double
solid_get_calorific_capacity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return 2.0;
+ return ((struct solid*)sdis_data_cget(data))->cp;
}
static double
solid_get_thermal_conductivity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return LAMBDA;
+ return ((struct solid*)sdis_data_cget(data))->lambda;
}
static double
solid_get_volumic_mass
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return 25.0;
+ return ((struct solid*)sdis_data_cget(data))->rho;
}
static double
solid_get_delta
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
- CHK(vtx != NULL);
- return 1.0/20.0;
-}
-
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
CHK(vtx != NULL);
- return 2.1/20.0;
+ return ((struct solid*)sdis_data_cget(data))->delta;
}
static double
@@ -222,6 +217,7 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid2 = NULL; /* For debug */
struct sdis_interface* interf_adiabatic = NULL;
struct sdis_interface* interf_T0 = NULL;
struct sdis_scene* box_scn = NULL;
@@ -233,6 +229,7 @@ main(int argc, char** argv)
struct sdis_interface* box_interfaces[12 /*#triangles*/];
struct sdis_interface* square_interfaces[4/*#segments*/];
struct interf* interf_props = NULL;
+ struct solid* solid_props = NULL;
double pos[3];
double x;
double ref;
@@ -242,21 +239,37 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
- /* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
- /* Create the solid_medium */
+ /* Setup the solid shader */
solid_shader.calorific_capacity = solid_get_calorific_capacity;
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
solid_shader.volumic_power = solid_get_volumic_power;
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Create the solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data) == RES_OK);
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ CHK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data) == RES_OK);
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = 0;
+ solid_props->cp = 0;
+ solid_props->rho = 0;
+ solid_props->delta = DELTA/4;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
/* Setup the interface shader */
interf_shader.convection_coef = interface_get_convection_coef;
@@ -280,6 +293,7 @@ main(int argc, char** argv)
/* Release the media */
CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
CHK(sdis_medium_ref_put(fluid) == RES_OK);
/* Map the interfaces to their box triangles */
@@ -324,19 +338,22 @@ main(int argc, char** argv)
printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Solve in 2D */
CHK(sdis_solve_probe(square_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2.0));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
CHK(sdis_scene_ref_put(square_scn) == RES_OK);
diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c
@@ -0,0 +1,468 @@
+/* Copyright (C) 2016-2018 |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 <rsys/math.h>
+
+#define N 10000 /* #realisations */
+#define Pw 10000 /* Volumic power */
+#define NONE -1
+#define DELTA 0.01
+#define DELTA_PSQUARE 0.01
+
+struct reference {
+ double pos[3];
+ double temperature_2d; /* In celcius */
+ double temperature_3d; /* In celcius */
+};
+
+/* Temperature in Celcius. The reference is computed by EDF with Syrthes
+ * #realisations: 100000
+ *
+ * >>> Check 1
+ * 0.85 0 = 190.29 ~ 190.198 +/- 0.572596; #failures: 46
+ * 0.65 0 = 259.95 ~ 259.730 +/- 0.678251; #failures: 73
+ * 0.45 0 = 286.33 ~ 285.287 +/- 0.693572; #failures: 74
+ * 0.25 0 = 235.44 ~ 235.672 +/- 0.710927; #failures: 61
+ * 0.05 0 = 192.33 ~ 192.464 +/- 0.693148; #failures: 70
+ *-0.15 0 = 156.82 ~ 157.526 +/- 0.668902; #failures: 43
+ *-0.35 0 = 123.26 ~ 124.234 +/- 0.634061; #failures: 31
+ *-0.55 0 = 90.250 ~ 91.0285 +/- 0.566423; #failures: 32
+ *
+ * >>> Check 2
+ * 0.85 0 = 678.170 ~ 671.302 +/- 4.03424; #failures: 186
+ * 0.65 0 = 1520.84 ~ 1523.42 +/- 5.38182; #failures: 442
+ * 0.45 0 = 1794.57 ~ 1790.60 +/- 5.44808; #failures: 528
+ * 0.25 0 = 1429.74 ~ 1419.80 +/- 5.33467; #failures: 406 */
+
+static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = {
+ -0.5,-1.0,-0.5,
+ -0.5, 1.0,-0.5,
+ 0.5, 1.0,-0.5,
+ 0.5,-1.0,-0.5,
+ -0.5,-1.0, 0.5,
+ -0.5, 1.0, 0.5,
+ 0.5, 1.0, 0.5,
+ 0.5,-1.0, 0.5,
+ -0.1, 0.4,-0.5,
+ -0.1, 0.6,-0.5,
+ 0.1, 0.6,-0.5,
+ 0.1, 0.4,-0.5,
+ -0.1, 0.4, 0.5,
+ -0.1, 0.6, 0.5,
+ 0.1, 0.6, 0.5,
+ 0.1, 0.4, 0.5
+};
+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 */
+ 1, 5, 6, 6, 2, 1, /* Cuboid top */
+ 6, 7, 3, 3, 2, 6, /* Cuboid right */
+ 0, 3, 7, 7, 4, 0, /* Cuboid bottom */
+ /* Cuboid back */
+ 0, 1, 9, 9, 8, 0,
+ 1, 2, 10, 10, 9, 1,
+ 2, 3, 11, 11, 10, 2,
+ 3, 0, 8, 8, 11, 3,
+ /* Cuboid front */
+ 5, 4, 12, 12, 13, 5,
+ 5, 13, 14, 14, 6, 5,
+ 6, 14, 15, 15, 7, 6,
+ 7, 15, 12, 12, 4, 7,
+ 8, 12, 13, 13, 9, 8, /* Cube left */
+ 9, 13, 14, 14, 10, 9, /* Cube top */
+ 14, 15, 11, 11, 10, 14, /* Cube right */
+ 8, 11, 15, 15, 12, 8, /* Cube bottom */
+ 8, 9, 10, 10, 11, 8, /* Cube back */
+ 12, 15, 14, 14, 13, 12 /* Cube front */
+};
+static const size_t ntriangles = sizeof(indices)/sizeof(size_t[3]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[itri*3+0];
+ ids[1] = indices[itri*3+1];
+ ids[2] = indices[itri*3+2];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*3+0];
+ pos[1] = vertices[ivert*3+1];
+ pos[2] = vertices[ivert*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[itri/2];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double P;
+ double T;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->T;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->P;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double pos[3] = {0,0};
+ size_t i;
+
+ FOR_EACH(i, 0, nrefs) {
+ double Tc;
+ pos[0] = refs[i].pos[0];
+ pos[1] = refs[i].pos[1];
+ pos[2] = refs[i].pos[2];
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ Tc = T.E - 273.15; /* Convert in Celcius */
+ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT3(pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+ /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+}
+
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_solid1_adiabatic = NULL;
+ struct sdis_interface* interf_solid2_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid1 = NULL;
+ struct sdis_interface* interf_solid1_fluid2 = NULL;
+ struct sdis_interface* interfaces[18 /*#rectangles*/];
+ /* In celcius. Computed by EDF with Syrthes */
+ const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
+ {{0, 0.85, 0}, 190.29, 189.13},
+ {{0, 0.65, 0}, 259.95, 247.09},
+ {{0, 0.45, 0}, 286.33, 308.42},
+ {{0, 0.25, 0}, 235.44, 233.55},
+ {{0, 0.05, 0}, 192.33, 192.30},
+ {{0,-0.15, 0}, 156.82, 156.98},
+ {{0,-0.35, 0}, 123.26, 123.43},
+ {{0,-0.55, 0}, 90.250, 90.040}
+ };
+ const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
+ {{0, 0.85}, 678.170, -1},
+ {{0, 0.65}, 1520.84, -1},
+ {{0, 0.45}, 1794.57, -1},
+ {{0, 0.25}, 1429.74, -1}
+ };
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ /* Create the fluid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 373.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 273.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 1;
+ solid_param->delta = DELTA;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 10;
+ solid_param->delta = DELTA_PSQUARE;
+ solid_param->P = Pw;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL,
+ NULL, &interf_solid1_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interfaces */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, fluid1, &interf_shader, data,
+ &interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid1/fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 5;
+ interf_param->temperature = NONE;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/fluid2 interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 10;
+ interf_param->temperature = NONE;
+ CHK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data,
+ &interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Map the interfaces to their faces */
+ interfaces[0] = interf_solid1_adiabatic;
+ interfaces[1] = interf_solid1_fluid1;
+ interfaces[2] = interf_solid1_adiabatic;
+ interfaces[3] = interf_solid1_fluid2;
+ interfaces[4] = interf_solid1_adiabatic;
+ interfaces[5] = interf_solid1_adiabatic;
+ interfaces[6] = interf_solid1_adiabatic;
+ interfaces[7] = interf_solid1_adiabatic;
+ interfaces[8] = interf_solid1_adiabatic;
+ interfaces[9] = interf_solid1_adiabatic;
+ interfaces[10] = interf_solid1_adiabatic;
+ interfaces[11] = interf_solid1_adiabatic;
+ interfaces[12] = interf_solid1_solid2;
+ interfaces[13] = interf_solid1_solid2;
+ interfaces[14] = interf_solid1_solid2;
+ interfaces[15] = interf_solid1_solid2;
+ interfaces[16] = interf_solid2_adiabatic;
+ interfaces[17] = interf_solid2_adiabatic;
+
+ /* Create the scene */
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+#if 0
+ dump_mesh(stdout, vertices, nvertices, indices, ntriangles);
+ exit(0);
+#endif
+
+ printf(">>> Check 1\n");
+ check(scn, refs1, sizeof(refs1)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 0.1;
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 2\n");
+ check(scn, refs2, sizeof(refs2)/sizeof(struct reference));
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_solid2) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c
@@ -0,0 +1,502 @@
+/* Copyright (C) 2016-2018 |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 <rsys/math.h>
+
+#define N 10000 /* #realisations */
+#define Pw 10000 /* Volumic power */
+#define NONE -1
+
+/* H delta T */
+#define Tboundary1 NONE
+#define Tboundary2 NONE
+#define DELTA 0.01
+#define Tref 286.83 /* In Celsius. Computed with Syrthes at the position 0.5 */
+
+/* Dirichlets */
+/*#define Tboundary1 373.15*/
+/*#define Tboundary2 273.15*/
+/*#define DELTA 0.01*/
+/*#define Tref 246.93*/ /* In Celsius. Computed with Syrthes at the position 0.5 */
+
+/* Temperature in Celcius. The reference is computed by EDF with Syrthes
+ * #realisations: 100000
+ *
+ * >>> Check1
+ * 0.85 = 190.29 ~ 189.322 +/- 0.566717; #failures: 51
+ * 0.65 = 259.95 ~ 259.995 +/- 0.674453; #failures: 82
+ * 0.45 = 286.33 ~ 285.928 +/- 0.691044; #failures: 76
+ * 0.25 = 235.44 ~ 234.672 +/- 0.700354; #failures: 80
+ * 0.05 = 192.33 ~ 191.977 +/- 0.690793; #failures: 64
+ *-0.15 = 156.82 ~ 155.765 +/- 0.660722; #failures: 40
+ *-0.35 = 123.26 ~ 122.973 +/- 0.621093; #failures: 29
+ *-0.55 = 90.250 ~ 90.3501 +/- 0.561255; #failures: 27
+ *
+ * >>> Check 2
+ * 0.85 = 678.170 ~ 662.616 +/- 3.97997; #failures: 221
+ * 0.65 = 1520.84 ~ 1486.35 +/- 5.25785; #failures: 474
+ * 0.45 = 1794.57 ~ 1767.21 +/- 5.36318; #failures: 584
+ * 0.25 = 1429.74 ~ 1401.39 +/- 5.25579; #failures: 465
+ *
+ * >>> Check 3
+ * 0.85 = 83.99 ~ 84.0098 +/- 0.115932; #failures: 51
+ * 0.65 = 73.90 ~ 73.9596 +/- 0.138835; #failures: 82
+ * 0.45 = 68.43 ~ 70.0292 +/- 0.144928; #failures: 76
+ * 0.25 = 60.61 ~ 61.4412 +/- 0.153980; #failures: 80
+ * 0.05 = 52.09 ~ 51.9452 +/- 0.158045; #failures: 64
+ *-0.15 = 42.75 ~ 42.9072 +/- 0.156546; #failures: 40
+ *-0.35 = 33.04 ~ 33.9338 +/- 0.149751; #failures: 29
+ *-0.55 = 24.58 ~ 24.7237 +/- 0.136441; #failures: 27 */
+
+/*
+ * _\ T1
+ * / /
+ * \__/
+ * ///+-----H1-------+///
+ * ///| |///
+ * ///| +------+ |///
+ * ///| |LAMBDA| |///
+ * ///| | Pw | |///
+ * ///| +------+ |///
+ * ///| |///
+ * ///| |///
+ * ///| LAMBDA1 |///
+ * ///| |///
+ * ///| |///
+ * ///| |///
+ * ///+-----H2-------+///
+ * _\ T2
+ * / /
+ * \__/
+ */
+
+struct reference {
+ double pos[2];
+ double temperature; /* In celcius */
+};
+
+static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
+ -0.5,-1.0,
+ -0.5, 1.0,
+ 0.5, 1.0,
+ 0.5,-1.0,
+ -0.1, 0.4,
+ -0.1, 0.6,
+ 0.1, 0.6,
+ 0.1, 0.4
+};
+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 */
+ 1, 2, /* Rectangle top */
+ 2, 3, /* Rectangle right */
+ 3, 0, /* Rectangle bottom */
+ 4, 5, /* Square left */
+ 5, 6, /* Square top */
+ 6, 7, /* Square right */
+ 7, 4 /* Square bottom */
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double P;
+ double T;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->T;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->P;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double pos[2] = {0,0};
+ size_t i;
+
+ FOR_EACH(i, 0, nrefs) {
+ double Tc;
+ pos[0] = refs[i].pos[0];
+ pos[1] = refs[i].pos[1];
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ Tc = T.E - 273.15; /* Convert in Celcius */
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), refs[i].temperature, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+ /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid1 = NULL;
+ struct sdis_interface* interf_solid1_fluid2 = NULL;
+ struct sdis_interface* interfaces[8 /*#segment*/];
+
+ /* In celcius. Computed by EDF with Syrthes */
+ const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
+ {{0, 0.85}, 190.29},
+ {{0, 0.65}, 259.95},
+ {{0, 0.45}, 286.33},
+ {{0, 0.25}, 235.44},
+ {{0, 0.05}, 192.33},
+ {{0,-0.15}, 156.82},
+ {{0,-0.35}, 123.26},
+ {{0,-0.55}, 90.250}
+ };
+ const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
+ {{0, 0.85}, 678.17},
+ {{0, 0.65}, 1520.84},
+ {{0, 0.45}, 1794.57},
+ {{0, 0.25}, 1429.74}
+ };
+ const struct reference refs3[] = { /* Lambda1=1, Lambda2=10, Pw=NONE */
+ {{0, 0.85}, 83.99},
+ {{0, 0.65}, 73.90},
+ {{0, 0.45}, 68.43},
+ {{0, 0.25}, 60.61},
+ {{0, 0.05}, 52.09},
+ {{0,-0.15}, 42.75},
+ {{0,-0.35}, 33.04},
+ {{0,-0.55}, 24.58}
+ };
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ /* Create the fluid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 373.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 273.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 1;
+ solid_param->delta = DELTA;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 10;
+ solid_param->delta = DELTA;
+ solid_param->P = Pw;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL,
+ NULL, &interf_solid1_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid1/fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 5;
+ interf_param->temperature = Tboundary1;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/fluid2 interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 10;
+ interf_param->temperature = Tboundary2;
+ CHK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data,
+ &interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_adiabatic;
+ interfaces[1] = interf_solid1_fluid1;
+ interfaces[2] = interf_adiabatic;
+ interfaces[3] = interf_solid1_fluid2;
+ interfaces[4] = interf_solid1_solid2;
+ interfaces[5] = interf_solid1_solid2;
+ interfaces[6] = interf_solid1_solid2;
+ interfaces[7] = interf_solid1_solid2;
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf(">>> Check 1\n");
+ check(scn, refs1, sizeof(refs1)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 0.1;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 2\n");
+ check(scn, refs2, sizeof(refs2)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 1;
+ data = sdis_medium_get_data(solid2);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 10;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 3\n");
+ check(scn, refs3, sizeof(refs3)/sizeof(struct reference));
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_solid2) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c
@@ -0,0 +1,465 @@
+/* Copyright (C) 2016-2018 |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 <rsys/math.h>
+
+#define Pw 10000.0 /* Volumic power */
+#define LAMBDA 10.0 /* Lambda of the middle slab */
+#define LAMBDA1 1.0 /* Lambda of the upper slab */
+#define LAMBDA2 LAMBDA1 /* Lambda of the lower slab */
+#define T1 373.15 /* Temperature of the upper fluid */
+#define T2 273.15 /* Temperature of the lower fluid */
+#define H1 5.0 /* Convection coef between the upper slab and the fluid */
+#define H2 10.0 /* Convection coef between the lower slab and the fluid */
+#define DELTA 0.01 /* Delta of the middle slab */
+#define DELTA1 0.02 /* Delta of the upper slab */
+#define DELTA2 0.07 /* Delta of the lower slab */
+#define L 0.2 /* Size of the middle slab */
+#define L1 0.4 /* Size of the upper slab */
+#define L2 1.4 /* Size of the lower slab */
+
+#define N 10000 /* #realisations */
+
+/* Analitically computed temperatures wrt the previous parameters.*/
+#define Tp1 648.6217
+#define Tp2 335.4141
+#define Ta 1199.5651
+#define Tb 1207.1122
+
+/* Fixed temperatures */
+#define UNKNOWN_TEMPERATURE -1
+#define Tsolid1_fluid UNKNOWN_TEMPERATURE /*Tp1*/
+#define Tsolid2_fluid UNKNOWN_TEMPERATURE /*Tp2*/
+#define Tsolid_solid1 UNKNOWN_TEMPERATURE /*Ta*/
+#define Tsolid_solid2 UNKNOWN_TEMPERATURE /*Tb*/
+
+#define PROBE_POS 1.8
+
+/*
+ * The 2D scene is composed of 3 stacked solid slabs whose middle slab has a
+ * volumic power. The +/-X sides of the slabs are stretched far away to
+ * simulate a 1D case. The upper and lower bounds of the "sandwich" has a
+ * convective exchange with the surrounding fluid whose temperature is known.
+ *
+ * _\ T1
+ * / /
+ * \__/
+ * ... -----H1------ ... Tp1
+ * LAMBDA1
+ *
+ * ... ------------- ... Ta
+ * LAMBDA, Pw
+ * ... ------------- ... Tb
+ *
+ * LAMBDA2
+ *
+ *
+ * ... -----H2------ ... Tp2
+ * _\ T2
+ * / /
+ * \__/
+ */
+
+static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
+ -100000.5, 0.0, /* 0 */
+ -100000.5, 1.4, /* 1 */
+ -100000.5, 1.6, /* 2 */
+ -100000.5, 2.0, /* 3 */
+ 100000.5, 2.0, /* 4 */
+ 100000.5, 1.6, /* 5 */
+ 100000.5, 1.4, /* 6 */
+ 100000.5, 0.0 /* 7 */
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
+
+static const size_t indices[10/*#segments*/*2/*#indices per segment*/]= {
+ 0, 1,
+ 1, 2,
+ 2, 3,
+ 3, 4,
+ 4, 5,
+ 5, 6,
+ 6, 7,
+ 7, 0,
+ 6, 1,
+ 2, 5
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double volumic_power;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->volumic_power;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature_lower;
+ double temperature_upper;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct fluid* fluid;
+ CHK(data != NULL && vtx != NULL);
+ fluid = sdis_data_cget(data);
+ return vtx->P[1] < 1 ? fluid->temperature_lower : fluid->temperature_upper;
+}
+
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_solid_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_adiabatic = NULL;
+ struct sdis_interface* interf_solid2_adiabatic = NULL;
+ struct sdis_interface* interf_solid_solid1 = NULL;
+ struct sdis_interface* interf_solid_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid = NULL;
+ struct sdis_interface* interf_solid2_fluid = NULL;
+ struct sdis_interface* interfaces[10/*#segment*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ double Tref;
+ double pos[2];
+ size_t nfails;
+ size_t nreals;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature_upper = T1;
+ fluid_param->temperature_lower = T2;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the medium of the upper slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA1;
+ solid_param->delta = DELTA1;
+ solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the medium of the lower slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA2;
+ solid_param->delta = DELTA2;
+ solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the medium of the middle slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA;
+ solid_param->delta = DELTA;
+ solid_param->volumic_power = Pw;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid/solid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->temperature = Tsolid_solid1;
+ CHK(sdis_interface_create(dev, solid, solid1, &interf_shader,
+ data, &interf_solid_solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->temperature = Tsolid_solid2;
+ CHK(sdis_interface_create(dev, solid, solid2, &interf_shader,
+ data, &interf_solid_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interfaces */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid, &interf_shader, data,
+ &interf_solid_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data,
+ &interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data,
+ &interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1 fluid interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H1;
+ interf_param->temperature = Tsolid1_fluid;
+ CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data,
+ &interf_solid1_fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 fluid interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H2;
+ interf_param->temperature = Tsolid2_fluid;
+ CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data,
+ &interf_solid2_fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_solid2_adiabatic;
+ interfaces[1] = interf_solid_adiabatic;
+ interfaces[2] = interf_solid1_adiabatic;
+ interfaces[3] = interf_solid1_fluid;
+ interfaces[4] = interf_solid1_adiabatic;
+ interfaces[5] = interf_solid_adiabatic;
+ interfaces[6] = interf_solid2_adiabatic;
+ interfaces[7] = interf_solid2_fluid;
+ interfaces[8] = interf_solid_solid2;
+ interfaces[9] = interf_solid_solid1;
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_solid_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_solid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_solid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_fluid) == RES_OK);
+
+ pos[0] = 0;
+ pos[1] = PROBE_POS;
+
+ if(pos[1] > 0 && pos[1] < L2) { /* Lower slab */
+ Tref = Tp2 + (Tb - Tp2) * pos[1] / L2;
+ } else if(pos[1] > L2 && pos[1] < L2 + L) { /* Middle slab */
+ Tref =
+ (Ta + Tb) / 2
+ + (Ta - Tb)/L * (pos[1] - (L2+L/2))
+ + Pw * (L*L/4.0 - pow((pos[1] - (L2+L/2)), 2)) / (2*LAMBDA);
+ } else if(pos[1] > L2 + L && pos[1] < L2 + L1 + L) {
+ Tref = Ta + (Tp1 - Ta) / L1 * (pos[1] - (L+L2));
+ } else {
+ FATAL("Unreachable code.\n");
+ }
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c
@@ -0,0 +1,371 @@
+/* Copyright (C) 2016-2018 |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 <rsys/math.h>
+
+#define Tf1 0
+#define Tf2 100
+#define Power 0 /*10000*/
+#define H1 50
+#define H2 50
+#define LAMBDA 100.0
+#define DELTA (1.0/20.0)
+#define N 10000
+
+/*
+ * The 2D scene is a solid slabs stretched along the X dimension to simulate a
+ * 1D case. The slab has a volumic power and has a convective exchange with the
+ * surrounding fluid whose temperature is fixed to Tfluid.
+ *
+ *
+ * _\ TFluid
+ * / /
+ * \__/
+ *
+ * ... -----Hboundary----- ...
+ *
+ * Lambda, Power
+ *
+ * ... -----Hboundary----- ...
+ *
+ * _\ TFluid
+ * / /
+ * \__/
+ *
+ */
+
+static const double vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
+ -10000.5,-0.5,
+ -10000.5, 0.5,
+ 10000.5, 0.5,
+ 10000.5,-0.5
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
+
+static const size_t indices[4/*#segments*/*2/*#indices per segment*/]= {
+ 0, 1,
+ 1, 2,
+ 2, 3,
+ 3, 0
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double volumic_power;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->volumic_power;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct fluid* fluid;
+ CHK(data != NULL && vtx != NULL);
+ fluid = sdis_data_cget(data);
+ return fluid->temperature;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_solid_fluid1 = NULL;
+ struct sdis_interface* interf_solid_fluid2 = NULL;
+ struct sdis_interface* interfaces[4/*#segment*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals, nfails;
+ double pos[2];
+ double Tref;
+ double a, b, x;
+ double L;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf1;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf2;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA;
+ solid_param->delta = DELTA;
+ solid_param->volumic_power = Power;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H1;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_solid_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid fluid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H2;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
+ &interf_solid_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_adiabatic;
+ interfaces[1] = interf_solid_fluid1;
+ interfaces[2] = interf_adiabatic;
+ interfaces[3] = interf_solid_fluid2;
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_fluid2) == RES_OK);
+
+ pos[0] = 0;
+ pos[1] = 0.25;
+
+ L = vertices[3] - vertices[1];
+#if 1
+ x = pos[1] + vertices[3];
+ a = (H2*Power*L + H1*H2*(Tf1 - Tf2) + H1*H2*Power*L*L/(2*LAMBDA))
+ / (LAMBDA * (H1 + H2) + H1*H2*L);
+ b = Tf2 + a * LAMBDA / H2;
+ Tref = -Power / (2*LAMBDA) * x*x + a * x + b;
+#else
+ tmp = LAMBDA / L;
+ T1 = H1 * (H2+tmp) / (tmp*(H1+H2) + H1*H2) * Tf1
+ + H2 * tmp / (tmp*(H1+H2) + H1*H2) * Tf2;
+ T2 = H1 * tmp / (tmp*(H1+H2) + H1*H2) * Tf1
+ + H2 * (H1+tmp) / (tmp*(H1+H2) + H1*H2) * Tf2;
+ Tref = T2 + (T1-T2)/L * (pos[1] + vertices[3]);
+#endif
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
