commit b5e6bb5111d0326a86159f6790df2b13d5e52096
parent 5139c69ba3a0985bfc89497a020befc691968c01
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 3 Jul 2019 16:15:35 +0200
Merge branch 'release_0.8'
Diffstat:
49 files changed, 3016 insertions(+), 1150 deletions(-)
diff --git a/README.md b/README.md
@@ -25,6 +25,19 @@ variable the install directories of its dependencies.
## Release notes
+### Version 0.8
+
+- Drastically improve the robustness of the solver~: far less realisations are
+ now rejected.
+- Add the estimation of the time spent per realisation estimate. Add the
+ `sdis_estimator_get_realisation_time` function that returns this estimate.
+- Add the `sdis_estimator_buffer` API~: it manages a two dimensional array of
+ regular estimators and provides global estimations over the whole estimators
+ saved into the buffer.
+- Update the signature of the `sdis_solve_camera` function~: it now returns a
+ `sdis_estimator_buffer`. It now also supports time integration as well as
+ heat paths registration.
+
### Version 0.7
#### Add Green function support
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -16,9 +16,6 @@
cmake_minimum_required(VERSION 3.0)
project(stardis C)
enable_testing()
-if(POLICY CMP0054)
- cmake_policy(SET CMP0054 OLD) # Disable OpenMP CMake warning
-endif()
include(CMakeDependentOption)
@@ -32,12 +29,12 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS
# Check dependencies
################################################################################
find_package(RCMake 0.4 REQUIRED)
-find_package(Star2D 0.3 REQUIRED)
+find_package(Star2D 0.3.1 REQUIRED)
find_package(Star3D 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarEnc 0.2.2 REQUIRED)
find_package(StarEnc2D 0.2.2 REQUIRED)
-find_package(RSys 0.6.1 REQUIRED)
+find_package(RSys 0.8.1 REQUIRED)
find_package(OpenMP 2.0 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
@@ -58,16 +55,16 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP StarEnc StarEnc2D)
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 7)
-set(VERSION_PATCH 1)
+set(VERSION_MINOR 8)
+set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SDIS_FILES_SRC
- sdis_accum_buffer.c
sdis_camera.c
sdis_data.c
sdis_device.c
sdis_estimator.c
+ sdis_estimator_buffer.c
sdis_green.c
sdis_heat_path.c
sdis_interface.c
@@ -169,7 +166,6 @@ if(NOT NO_TEST)
register_test(${_name} ${_name})
endfunction()
- new_test(test_sdis_accum_buffer)
new_test(test_sdis_camera)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
@@ -181,6 +177,7 @@ if(NOT NO_TEST)
new_test(test_sdis_interface)
new_test(test_sdis_medium)
new_test(test_sdis_scene)
+ new_test(test_sdis_solid_random_walk_robustness)
new_test(test_sdis_solve_camera)
new_test(test_sdis_solve_probe)
new_test(test_sdis_solve_probe2)
@@ -193,7 +190,7 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_medium)
new_test(test_sdis_solve_medium_2d)
new_test(test_sdis_volumic_power)
- new_test(test_sdis_volumic_power4_2d)
+ new_test(test_sdis_volumic_power4)
# Additionnal tests
build_test(test_sdis_volumic_power2)
@@ -206,6 +203,7 @@ if(NOT NO_TEST)
add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
endif()
+ target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT)
target_link_libraries(test_sdis_solve_medium Star3DUT)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB})
diff --git a/src/sdis.h b/src/sdis.h
@@ -56,11 +56,11 @@ struct senc_descriptor;
* get or release a reference on the data, i.e. they increment or decrement the
* reference counter, respectively. When this counter reaches 0, the object is
* silently destroyed and cannot be used anymore. */
-struct sdis_accum_buffer;
struct sdis_camera;
struct sdis_data;
struct sdis_device;
struct sdis_estimator;
+struct sdis_estimator_buffer;
struct sdis_green_function;
struct sdis_interface;
struct sdis_medium;
@@ -137,14 +137,6 @@ struct sdis_interface_fragment {
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
-/* Monte-Carlo accumulator */
-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 */
struct sdis_mc {
double E; /* Expected value */
@@ -244,22 +236,6 @@ struct sdis_interface_shader {
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
-struct sdis_accum_buffer_layout {
- size_t width;
- size_t height;
-};
-#define SDIS_ACCUM_BUFFER_LAYOUT_NULL__ {0, 0}
-static const struct sdis_accum_buffer_layout SDIS_ACCUM_BUFFER_LAYOUT_NULL =
- SDIS_ACCUM_BUFFER_LAYOUT_NULL__;
-
-/* Functor use to write accumulations performed by sdis_solve_camera */
-typedef res_T
-(*sdis_write_accums_T)
- (void* context, /* User data */
- const size_t origin[2], /* Coordinates of the 1st accumulation in image plane */
- const size_t naccums[2], /* #accumulations in X and Y */
- const struct sdis_accum* accums); /* List of row ordered accumulations */
-
/* Vertex of heat path v*/
struct sdis_heat_vertex {
double P[3];
@@ -418,51 +394,47 @@ sdis_camera_look_at
const double up[3]);
/*******************************************************************************
- * A buffer of accumulations is a 2D array whose each cell stores an
- * Monte-Carlo accumulation, i.e. a sum of MC weights, the sum of their square
- * and the overall number of summed weights (see struct sdis_accum)
- ******************************************************************************/
+ * An estimator buffer is 2D array of estimators
+******************************************************************************/
SDIS_API res_T
-sdis_accum_buffer_create
- (struct sdis_device* dev,
- const size_t width, /* #cells in X */
- const size_t height, /* #cells in Y */
- struct sdis_accum_buffer** buf);
+sdis_estimator_buffer_ref_get
+ (struct sdis_estimator_buffer* buf);
SDIS_API res_T
-sdis_accum_buffer_ref_get
- (struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_ref_put
+ (struct sdis_estimator_buffer* buf);
SDIS_API res_T
-sdis_accum_buffer_ref_put
- (struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_get_definition
+ (const struct sdis_estimator_buffer* buf,
+ size_t definition[2]);
SDIS_API res_T
-sdis_accum_buffer_get_layout
- (const struct sdis_accum_buffer* buf,
- struct sdis_accum_buffer_layout* layout);
+sdis_estimator_buffer_at
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y,
+ const struct sdis_estimator** estimator);
+
+SDIS_API res_T
+sdis_estimator_buffer_get_realisation_count
+ (const struct sdis_estimator_buffer* buf,
+ size_t* nrealisations); /* Successful ones */
-/* Get a read only pointer toward the memory space of the accum buffer */
SDIS_API res_T
-sdis_accum_buffer_map
- (const struct sdis_accum_buffer* buf,
- const struct sdis_accum** accums);
+sdis_estimator_buffer_get_failure_count
+ (const struct sdis_estimator_buffer* buf,
+ size_t* nfailures);
SDIS_API res_T
-sdis_accum_buffer_unmap
- (const struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_get_temperature
+ (const struct sdis_estimator_buffer* buf,
+ struct sdis_mc* temperature);
-/* Helper function that matches the `sdis_write_accums_T' functor type. One can
- * send this function directly to the sdis_solve_camera function, to fill the
- * accum buffer with the estimation of the radiative temperature that reaches
- * each pixel of an image whose definition matches the definition of the accum
- * buffer. */
SDIS_API res_T
-sdis_accum_buffer_write
- (void* buf, /* User data */
- const size_t origin[2], /* Coordinates of the 1st accum in image plane */
- const size_t naccum[2], /* #accum in X and Y */
- const struct sdis_accum* accums); /* List of row ordered accum */
+sdis_estimator_buffer_get_realisation_time
+ (const struct sdis_estimator_buffer* buf,
+ struct sdis_mc* time);
/*******************************************************************************
* A medium encapsulates the properties of either a fluid or a solid.
@@ -722,6 +694,11 @@ sdis_estimator_get_temperature
struct sdis_mc* temperature);
SDIS_API res_T
+sdis_estimator_get_realisation_time
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* time);
+
+SDIS_API res_T
sdis_estimator_get_convective_flux
(const struct sdis_estimator* estimator,
struct sdis_mc* flux);
@@ -800,6 +777,18 @@ sdis_green_path_get_limit_point
(struct sdis_green_path* path,
struct sdis_point* pt);
+/* Retrieve the number of "power terms" associated to a path. */
+SDIS_API res_T
+sdis_green_function_get_power_terms_count
+ (const struct sdis_green_path* path,
+ size_t* nterms);
+
+/* Retrieve the number of "flux terms" associated to a path. */
+SDIS_API res_T
+sdis_green_function_get_flux_terms_count
+ (const struct sdis_green_path* path,
+ size_t* nterms);
+
/* Iterate over all "power terms" associated to the path. Multiply each term
* by the power of their associated medium, that is assumed to be constant in
* time and space, gives the medium power registered along the path. */
@@ -914,15 +903,15 @@ SDIS_API res_T
sdis_solve_camera
(struct sdis_scene* scn,
const struct sdis_camera* cam, /* Point of view */
- const double time, /* Observation time */
+ const double time_range[2], /* Observation time */
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
const size_t width, /* Image definition in in X */
const size_t height, /* Image definition in Y */
const size_t spp, /* #samples per pixel */
- sdis_write_accums_T writer,
- void* writer_data);
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_estimator_buffer** buf);
SDIS_API res_T
sdis_solve_medium
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -22,14 +22,6 @@
struct green_path_handle;
struct sdis_heat_path;
-struct accum {
- double sum; /* Sum of MC weights */
- double sum2; /* Sum of square MC weights */
- size_t count; /* #accumulated MC weights */
-};
-#define ACCUM_NULL__ {0,0,0}
-static const struct accum ACCUM_NULL = ACCUM_NULL__;
-
struct rwalk_context {
struct green_path_handle* green_path;
struct sdis_heat_path* heat_path;
@@ -39,24 +31,6 @@ struct rwalk_context {
#define RWALK_CONTEXT_NULL__ {NULL, NULL, 0, 0}
static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
-static INLINE void
-sum_accums
- (const struct accum accums[],
- const size_t naccums,
- struct accum* accum)
-{
- struct accum acc = ACCUM_NULL;
- size_t i;
- ASSERT(accums && naccums && accum);
-
- FOR_EACH(i, 0, naccums) {
- acc.sum += accums[i].sum;
- acc.sum2 += accums[i].sum2;
- acc.count += accums[i].count;
- }
- *accum = acc;
-}
-
#endif /* SDIS_XD_BEGIN_H */
#ifdef SDIS_XD_BEGIN_H__
diff --git a/src/sdis_accum_buffer.c b/src/sdis_accum_buffer.c
@@ -1,160 +0,0 @@
-/* Copyright (C) 2016-2019 |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 "sdis_device_c.h"
-
-struct sdis_accum_buffer {
- struct sdis_accum* accums;
- size_t width;
- size_t height;
-
- ref_T ref;
- struct sdis_device* dev;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-accum_buffer_release(ref_T* ref)
-{
- struct sdis_accum_buffer* buf;
- struct sdis_device* dev;
- ASSERT(ref);
- buf = CONTAINER_OF(ref, struct sdis_accum_buffer, ref);
- dev = buf->dev;
- if(buf->accums) MEM_RM(dev->allocator, buf->accums);
- MEM_RM(dev->allocator, buf);
- SDIS(device_ref_put(dev));
-}
-
-/*******************************************************************************
- * Export functions
- ******************************************************************************/
-res_T
-sdis_accum_buffer_create
- (struct sdis_device* dev,
- const size_t width,
- const size_t height,
- struct sdis_accum_buffer** out_buf)
-{
- struct sdis_accum_buffer* buf = NULL;
- size_t sz;
- res_T res = RES_OK;
-
- if(!dev || !width || !height || !out_buf) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- buf = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_accum_buffer));
- if(!buf) {
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&buf->ref);
- SDIS(device_ref_get(dev));
- buf->dev = dev;
- buf->width = width;
- buf->height = height;
-
- sz = buf->width * buf->height * sizeof(struct sdis_accum);
- buf->accums = MEM_ALLOC_ALIGNED(dev->allocator, sz, 16);
- if(!buf->accums) {
- log_err(dev, "%s: could not allocate a buffer accumulations of %lux%lu.\n",
- FUNC_NAME, (unsigned long)width, (unsigned long)height);
- res = RES_MEM_ERR;
- goto error;
- }
-
-exit:
- if(out_buf) *out_buf = buf;
- return res;
-error:
- if(buf) {
- SDIS(accum_buffer_ref_put(buf));
- buf = NULL;
- }
- goto exit;
-}
-
-res_T
-sdis_accum_buffer_ref_get(struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- ref_get(&buf->ref);
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_ref_put(struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- ref_put(&buf->ref, accum_buffer_release);
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_get_layout
- (const struct sdis_accum_buffer* buf,
- struct sdis_accum_buffer_layout* layout)
-{
- if(!buf || !layout) return RES_BAD_ARG;
- layout->width = buf->width;
- layout->height = buf->height;
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_map
- (const struct sdis_accum_buffer* buf, const struct sdis_accum** accums)
-{
- if(!buf || !accums) return RES_BAD_ARG;
- *accums = buf->accums;
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_unmap(const struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- /* Do nothing */
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_write
- (void* buffer,
- const size_t origin[2],
- const size_t size[2],
- const struct sdis_accum* accums)
-{
- struct sdis_accum_buffer* buf = buffer;
- const struct sdis_accum* src_row = accums;
- size_t dst_iy;
-
- if(UNLIKELY(!buffer || !origin || !size || !accums)) return RES_BAD_ARG;
- if(UNLIKELY((origin[0] + size[0]) > buf->width)) return RES_BAD_ARG;
- if(UNLIKELY((origin[1] + size[1]) > buf->height)) return RES_BAD_ARG;
-
- FOR_EACH(dst_iy, origin[1], origin[1] + size[1]) {
- const size_t dst_irow = dst_iy*buf->width + origin[0];
- memcpy(buf->accums + dst_irow, src_row, size[0] * sizeof(struct sdis_accum));
- src_row += size[0];
- }
- return RES_OK;
-}
-
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -54,7 +54,6 @@ device_release(ref_T* ref)
ASSERT(flist_name_is_empty(&dev->media_names));
flist_name_release(&dev->interfaces_names);
flist_name_release(&dev->media_names);
- darray_tile_release(&dev->tiles);
MEM_RM(dev->allocator, dev);
}
@@ -98,14 +97,6 @@ sdis_device_create
ref_init(&dev->ref);
flist_name_init(allocator, &dev->interfaces_names);
flist_name_init(allocator, &dev->media_names);
- darray_tile_init(allocator, &dev->tiles);
-
- res = darray_tile_resize(&dev->tiles, dev->nthreads);
- if(res != RES_OK) {
- log_err(dev,
- "%s: could not allocate the per thread buffer of estimations.\n", FUNC_NAME);
- goto error;
- }
res = s2d_device_create(log, allocator, 0, &dev->s2d);
if(res != RES_OK) {
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -26,18 +26,6 @@ struct name { FITEM; };
#define FITEM_TYPE name
#include <rsys/free_list.h>
-#define DARRAY_NAME accum
-#define DARRAY_DATA struct sdis_accum
-#include <rsys/dynamic_array.h>
-
-#define DARRAY_NAME tile
-#define DARRAY_DATA struct darray_accum
-#define DARRAY_FUNCTOR_INIT darray_accum_init
-#define DARRAY_FUNCTOR_RELEASE darray_accum_release
-#define DARRAY_FUNCTOR_COPY darray_accum_copy
-#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_accum_copy_and_release
-#include <rsys/dynamic_array.h>
-
struct sdis_device {
struct logger* logger;
struct mem_allocator* allocator;
@@ -46,7 +34,6 @@ struct sdis_device {
struct flist_name interfaces_names;
struct flist_name media_names;
- struct darray_tile tiles;
struct s2d_device* s2d;
struct s3d_device* s3d;
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -82,12 +82,28 @@ sdis_estimator_get_failure_count
return RES_OK;
}
+#define SETUP_MC(Mc, Acc) { \
+ (Mc)->E = (Acc)->sum / (double)(Acc)->count; \
+ (Mc)->V = (Acc)->sum2 / (double)(Acc)->count - (Mc)->E*(Mc)->E; \
+ (Mc)->V = MMAX((Mc)->V, 0); \
+ (Mc)->SE = sqrt((Mc)->V / (double)(Acc)->count); \
+} (void)0
+
res_T
sdis_estimator_get_temperature
(const struct sdis_estimator* estimator, struct sdis_mc* mc)
{
if(!estimator || !mc) return RES_BAD_ARG;
- *mc = estimator->temperature;
+ SETUP_MC(mc, &estimator->temperature);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_get_realisation_time
+ (const struct sdis_estimator* estimator, struct sdis_mc* mc)
+{
+ if(!estimator || !mc) return RES_BAD_ARG;
+ SETUP_MC(mc, &estimator->realisation_time);
return RES_OK;
}
@@ -98,7 +114,7 @@ sdis_estimator_get_convective_flux
if(!estimator || !flux ||estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_CONVECTIVE];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_CONVECTIVE]);
return RES_OK;
}
@@ -109,7 +125,7 @@ sdis_estimator_get_radiative_flux
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_RADIATIVE];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_RADIATIVE]);
return RES_OK;
}
@@ -120,10 +136,12 @@ sdis_estimator_get_total_flux
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_TOTAL];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_TOTAL]);
return RES_OK;
}
+#undef SETUP_MC
+
res_T
sdis_estimator_get_paths_count
(const struct sdis_estimator* estimator, size_t* npaths)
diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c
@@ -0,0 +1,244 @@
+/* Copyright (C) 2016-2019 |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 "sdis_device_c.h"
+#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
+
+struct sdis_estimator_buffer {
+ struct sdis_estimator** estimators; /* Row major per pixe lestimators */
+ size_t width;
+ size_t height;
+
+ struct accum temperature;
+ struct accum realisation_time;
+ size_t nrealisations; /* #successes */
+ size_t nfailures;
+
+ ref_T ref;
+ struct sdis_device* dev;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+estimator_buffer_release(ref_T* ref)
+{
+ struct sdis_estimator_buffer* buf;
+ struct sdis_device* dev;
+ size_t i;
+ ASSERT(ref);
+
+ buf = CONTAINER_OF(ref, struct sdis_estimator_buffer, ref);
+ dev = buf->dev;
+
+ if(buf->estimators) {
+ FOR_EACH(i, 0, buf->width*buf->height) {
+ if(buf->estimators[i]) SDIS(estimator_ref_put(buf->estimators[i]));
+ }
+ MEM_RM(dev->allocator, buf->estimators);
+ }
+
+ MEM_RM(dev->allocator, buf);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_estimator_buffer_ref_get(struct sdis_estimator_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_get(&buf->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_ref_put(struct sdis_estimator_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_put(&buf->ref, estimator_buffer_release);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_definition
+ (const struct sdis_estimator_buffer* buf, size_t definition[2])
+{
+ if(!buf || !definition) return RES_BAD_ARG;
+ definition[0] = buf->width;
+ definition[1] = buf->height;
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_at
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y,
+ const struct sdis_estimator** estimator)
+{
+ if(!buf || x >= buf->width || y >= buf->height || !estimator)
+ return RES_BAD_ARG;
+ *estimator = estimator_buffer_grab(buf, x, y);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_realisation_count
+ (const struct sdis_estimator_buffer* buf, size_t* nrealisations)
+{
+ if(!buf || !nrealisations) return RES_BAD_ARG;
+ *nrealisations = buf->nrealisations;
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_failure_count
+ (const struct sdis_estimator_buffer* buf, size_t* nfailures)
+{
+ if(!buf || !nfailures) return RES_BAD_ARG;
+ *nfailures = buf->nfailures;
+ return RES_OK;
+}
+
+#define SETUP_MC(Mc, Acc) { \
+ (Mc)->E = (Acc)->sum / (double)(Acc)->count; \
+ (Mc)->V = (Acc)->sum2 / (double)(Acc)->count - (Mc)->E*(Mc)->E; \
+ (Mc)->V = MMAX((Mc)->V, 0); \
+ (Mc)->SE = sqrt((Mc)->V / (double)(Acc)->count); \
+} (void)0
+
+res_T
+sdis_estimator_buffer_get_temperature
+ (const struct sdis_estimator_buffer* buf, struct sdis_mc* mc)
+{
+ if(!buf || !mc) return RES_BAD_ARG;
+ SETUP_MC(mc, &buf->temperature);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_realisation_time
+(const struct sdis_estimator_buffer* buf, struct sdis_mc* mc)
+{
+ if(!buf || !mc) return RES_BAD_ARG;
+ SETUP_MC(mc, &buf->realisation_time);
+ return RES_OK;
+}
+
+#undef SETUP_MC
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+estimator_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_estimator_buffer** out_buf)
+{
+ struct sdis_estimator_buffer* buf = NULL;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!dev || !width || !height || !out_buf) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ buf = MEM_CALLOC(dev->allocator, 1, sizeof(*buf));
+ if(!buf) { res = RES_MEM_ERR; goto error; }
+
+ ref_init(&buf->ref);
+ SDIS(device_ref_get(dev));
+ buf->dev = dev;
+ buf->width = width;
+ buf->height = height;
+
+ buf->estimators = MEM_CALLOC
+ (dev->allocator, width*height, sizeof(*buf->estimators));
+ if(!buf->estimators) {
+ log_err(dev, "%s: could not allocate a buffer of estimators of %lux%lu.\n",
+ FUNC_NAME, (unsigned long)width, (unsigned long)height);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(i, 0, width*height) {
+ res = estimator_create(dev, SDIS_ESTIMATOR_TEMPERATURE, buf->estimators+i);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ if(out_buf) *out_buf = buf;
+ return res;
+error:
+ if(buf) {
+ SDIS(estimator_buffer_ref_put(buf));
+ buf = NULL;
+ }
+ goto exit;
+}
+
+struct sdis_estimator*
+estimator_buffer_grab
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y)
+{
+ ASSERT(x < buf->width && y < buf->height);
+ return buf->estimators[y*buf->width + x];
+}
+
+void
+estimator_buffer_setup_realisations_count
+ (struct sdis_estimator_buffer* buf,
+ const size_t nrealisations,
+ const size_t nsuccesses)
+{
+ ASSERT(buf && nrealisations && nsuccesses && nsuccesses<=nrealisations);
+ buf->nrealisations = nsuccesses;
+ buf->nfailures = nrealisations - nsuccesses;
+}
+
+void
+estimator_buffer_setup_temperature
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2)
+{
+ ASSERT(buf && buf->nrealisations);
+ buf->temperature.sum = sum;
+ buf->temperature.sum2 = sum2;
+ buf->temperature.count = buf->nrealisations;
+}
+
+void
+estimator_buffer_setup_realisation_time
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2)
+{
+ ASSERT(buf && buf->nrealisations);
+ buf->realisation_time.sum = sum;
+ buf->realisation_time.sum2 = sum2;
+ buf->realisation_time.count = buf->nrealisations;
+}
+
diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h
@@ -0,0 +1,58 @@
+/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_ESTIMATOR_BUFFER_C_H
+#define SDIS_ESTIMATOR_BUFFER_C_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration */
+struct sdis_device;
+struct sdis_estimator;
+struct sdis_estimator_buffer;
+
+extern LOCAL_SYM res_T
+estimator_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_estimator_buffer** buf);
+
+extern LOCAL_SYM struct sdis_estimator*
+estimator_buffer_grab
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_realisations_count
+ (struct sdis_estimator_buffer* buf,
+ const size_t nrealisations,
+ const size_t nsuccesses);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_temperature
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_realisation_time
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2);
+
+#endif /* SDIS_ESTIMATOR_BUFFER_C_H */
+
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -17,6 +17,7 @@
#define SDIS_ESTIMATOR_C_H
#include "sdis_heat_path.h"
+#include "sdis_misc.h"
#include <rsys/math.h>
#include <rsys/ref_count.h>
@@ -34,9 +35,10 @@ enum flux_name {
};
struct sdis_estimator {
- struct sdis_mc temperature;
- struct sdis_mc fluxes[FLUX_NAMES_COUNT__];
- size_t nrealisations;
+ struct accum temperature;
+ struct accum realisation_time;
+ struct accum fluxes[FLUX_NAMES_COUNT__];
+ size_t nrealisations; /* #successes */
size_t nfailures;
struct mutex* mutex;
@@ -82,13 +84,22 @@ estimator_setup_temperature
const double sum,
const double sum2)
{
- double N;
ASSERT(estim && estim->nrealisations);
- N = (double)estim->nrealisations;
- estim->temperature.E = sum/N;
- estim->temperature.V = sum2/N - estim->temperature.E*estim->temperature.E;
- estim->temperature.V = MMAX(estim->temperature.V, 0);
- estim->temperature.SE = sqrt(estim->temperature.V/N);
+ estim->temperature.sum = sum;
+ estim->temperature.sum2 = sum2;
+ estim->temperature.count = estim->nrealisations;
+}
+
+static INLINE void
+estimator_setup_realisation_time
+ (struct sdis_estimator* estim,
+ const double sum,
+ const double sum2)
+{
+ ASSERT(estim && estim->nrealisations);
+ estim->realisation_time.sum = sum;
+ estim->realisation_time.sum2 = sum2;
+ estim->realisation_time.count = estim->nrealisations;
}
static INLINE void
@@ -98,13 +109,10 @@ estimator_setup_flux
const double sum,
const double sum2)
{
- double N;
ASSERT(estim && (unsigned)name < FLUX_NAMES_COUNT__ && estim->nrealisations);
- N = (double)estim->nrealisations;
- estim->fluxes[name].E = sum/N;
- estim->fluxes[name].V = sum2/N - estim->fluxes[name].E*estim->fluxes[name].E;
- estim->fluxes[name].V = MMAX(estim->fluxes[name].V, 0);
- estim->fluxes[name].SE = sqrt(estim->fluxes[name].V/N);
+ estim->fluxes[name].sum = sum;
+ estim->fluxes[name].sum2 = sum2;
+ estim->fluxes[name].count = estim->nrealisations;
}
#endif /* SDIS_PROBE_ESTIMATOR_C_H */
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -192,6 +192,8 @@ struct sdis_green_function {
size_t npaths_valid;
size_t npaths_invalid;
+ struct accum realisation_time; /* Time per realisation */
+
struct ssp_rng_type rng_type;
FILE* rng_state;
@@ -489,6 +491,8 @@ sdis_green_function_solve
/* Setup the estimated temperature */
estimator_setup_realisations_count(estimator, npaths, N);
estimator_setup_temperature(estimator, accum, accum2);
+ estimator_setup_realisation_time
+ (estimator, green->realisation_time.sum, green->realisation_time.sum2);
exit:
if(rng) SSP(rng_ref_put(rng));
@@ -595,6 +599,60 @@ error:
}
res_T
+sdis_green_function_get_power_terms_count
+ (const struct sdis_green_path* path_handle,
+ size_t* nterms)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !nterms) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ green = path_handle->green__; (void)green;
+ ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
+
+ path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+
+ *nterms = darray_power_term_size_get(&path->power_terms);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+sdis_green_function_get_flux_terms_count
+ (const struct sdis_green_path* path_handle,
+ size_t* nterms)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !nterms) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ green = path_handle->green__; (void)green;
+ ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
+
+ path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+
+ *nterms = darray_flux_term_size_get(&path->flux_terms);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
sdis_green_path_for_each_power_term
(struct sdis_green_path* path_handle,
sdis_process_medium_power_term_T func,
@@ -781,7 +839,7 @@ green_function_redux_and_clear
{
size_t i;
res_T res = RES_OK;
- ASSERT(dst && greens && ngreens);
+ ASSERT(dst && greens);
FOR_EACH(i, 0, ngreens) {
res = green_function_merge_and_clear(dst, greens[i]);
@@ -797,12 +855,13 @@ error:
res_T
green_function_finalize
(struct sdis_green_function* green,
- struct ssp_rng_proxy* proxy)
+ struct ssp_rng_proxy* proxy,
+ const struct accum* time)
{
size_t i, n;
res_T res = RES_OK;
- if(!green || !proxy) {
+ if(!green || !proxy || !time) {
res = RES_BAD_ARG;
goto error;
}
@@ -821,6 +880,9 @@ green_function_finalize
}
green->npaths_invalid = n - green->npaths_valid;
+ ASSERT(green->npaths_valid == time->count);
+ green->realisation_time = *time;
+
exit:
return res;
error:
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -19,6 +19,7 @@
#include <rsys/rsys.h>
/* Forward declaration */
+struct accum;
struct sdis_green_function;
struct ssp_rng_proxy;
struct green_path;
@@ -53,7 +54,8 @@ green_function_redux_and_clear
extern LOCAL_SYM res_T
green_function_finalize
(struct sdis_green_function* green,
- struct ssp_rng_proxy* rng_proxy); /* Proxy RNG used to estimate the function */
+ struct ssp_rng_proxy* rng_proxy, /* Proxy RNG used to estimate the function */
+ const struct accum* time); /* Accumulator of the realisation time */
extern LOCAL_SYM res_T
green_function_create_path
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -25,15 +25,15 @@
#include "sdis_Xd_begin.h"
/* Emperical scale factor applied to the challenged reinjection distance. If
- * the distance to reinject is less than this adjusted value, the solver
- * switches from 2D reinjection scheme to the 1D reinjection scheme in order to
- * avoid numerical issues. */
+ * the distance to reinject is less than this adjusted value, the solver will
+ * try to discard the reinjection distance if possible in order to avoid
+ * numerical issues. */
#define REINJECT_DST_MIN_SCALE 0.125f
/*******************************************************************************
* Helper functions
******************************************************************************/
-static FINLINE void
+static void
XD(sample_reinjection_dir)
(const struct XD(rwalk)* rwalk, struct ssp_rng* rng, float dir[DIM])
{
@@ -75,6 +75,358 @@ XD(sample_reinjection_dir)
#endif
}
+#if DIM == 2
+static void
+XD(move_away_primitive_boundaries)
+ (struct XD(rwalk)* rwalk,
+ const double delta)
+{
+ struct sXd(attrib) attr;
+ float pos[DIM];
+ float dir[DIM];
+ float len;
+ const float st = 0.5f;
+ ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->hit) && delta > 0);
+
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr));
+
+ fX_set_dX(pos, rwalk->vtx.P);
+ fX(sub)(dir, attr.value, pos);
+ len = fX(normalize)(dir, dir);
+ len = MMIN(len, (float)(delta*0.1));
+
+ XD(move_pos)(rwalk->vtx.P, dir, len);
+}
+#else
+/* Move the random walk away from the primitive boundaries to avoid numerical
+ * issues leading to inconsistent random walks. */
+static void
+XD(move_away_primitive_boundaries)
+ (struct XD(rwalk)* rwalk,
+ const double delta)
+{
+ struct s3d_attrib v0, v1, v2; /* Triangle vertices */
+ float E[3][4]; /* 3D edge equations */
+ float dst[3]; /* Distance from current position to edge equation */
+ float N[3]; /* Triangle normal */
+ float P[3]; /* Random walk position */
+ float tmp[3];
+ float min_dst, max_dst;
+ float cos_a1, cos_a2;
+ float len;
+ int imax;
+ int imin;
+ int imid;
+ int i;
+ ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->hit));
+
+ fX_set_dX(P, rwalk->vtx.P);
+
+ /* Fetch triangle vertices */
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the edge vector */
+ f3_sub(E[0], v1.value, v0.value);
+ f3_sub(E[1], v2.value, v1.value);
+ f3_sub(E[2], v0.value, v2.value);
+
+ /* Compute the triangle normal */
+ f3_cross(N, E[1], E[0]);
+
+ /* Compute the 3D edge equation */
+ f3_normalize(E[0], f3_cross(E[0], E[0], N));
+ f3_normalize(E[1], f3_cross(E[1], E[1], N));
+ f3_normalize(E[2], f3_cross(E[2], E[2], N));
+ E[0][3] = -f3_dot(E[0], v0.value);
+ E[1][3] = -f3_dot(E[1], v1.value);
+ E[2][3] = -f3_dot(E[2], v2.value);
+
+ /* Compute the distance from current position to the edges */
+ dst[0] = f3_dot(E[0], P) + E[0][3];
+ dst[1] = f3_dot(E[1], P) + E[1][3];
+ dst[2] = f3_dot(E[2], P) + E[2][3];
+
+ /* Retrieve the min and max distance from random walk position to triangle
+ * edges */
+ min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
+ max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
+
+ /* Sort the edges with respect to their distance to the random walk position */
+ FOR_EACH(i, 0, 3) {
+ if(dst[i] == min_dst) {
+ imin = i;
+ } else if(dst[i] == max_dst) {
+ imax = i;
+ } else {
+ imid = i;
+ }
+ }
+ (void)imax;
+
+ /* TODO if the current position is near a vertex, one should move toward the
+ * farthest edge along its normal to avoid too small displacement */
+
+ /* Compute the distance `dst' from the current position to the edges to move
+ * to, along the normal of the edge from which the random walk is the nearest
+ *
+ * +. cos(a) = d / dst => dst = d / cos_a
+ * / `*.
+ * / | `*.
+ * / dst| a /`*.
+ * / | / `*.
+ * / | / d `*.
+ * / |/ `*.
+ * +---------o----------------+ */
+ cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
+ cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
+ dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
+ dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
+
+ /* Compute the maximum displacement distance into the triangle along the
+ * normal of the edge from which the random walk is the nearest */
+ len = MMIN(dst[imid], dst[imax]);
+ ASSERT(len != FLT_MAX);
+
+ /* Define the displacement distance as the minimum between 10 percent of
+ * delta and len / 2. */
+ len = MMIN(len*0.5f, (float)(delta*0.1));
+ XD(move_pos)(rwalk->vtx.P, E[imin], len);
+}
+#endif
+
+static res_T
+XD(select_reinjection_dir)
+ (const struct sdis_scene* scn,
+ const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
+ struct XD(rwalk)* rwalk, /* Current random walk state */
+ const float dir0[DIM], /* Challenged direction */
+ const float dir1[DIM], /* Challanged direction */
+ const double delta, /* Max reinjection distance */
+ float reinject_dir[DIM], /* Selected direction */
+ float* reinject_dst, /* Effective reinjection distance */
+ int can_move, /* Define of the random wal pos can be moved or not */
+ int* move_pos, /* Define if the current random walk was moved. May be NULL */
+ struct sXd(hit)* reinject_hit) /* Hit along the reinjection dir */
+{
+ struct sdis_interface* interf;
+ struct sdis_medium* mdm0;
+ struct sdis_medium* mdm1;
+ struct hit_filter_data filter_data;
+ struct sXd(hit) hit;
+ struct sXd(hit) hit0;
+ struct sXd(hit) hit1;
+ double tmp[DIM];
+ double dst;
+ double dst0;
+ double dst1;
+ const double delta_adjusted = delta * RAY_RANGE_MAX_SCALE;
+ const float* dir;
+ const float reinject_threshold = (float)delta * REINJECT_DST_MIN_SCALE;
+ float org[DIM];
+ float range[2];
+ enum sdis_side side;
+ int iattempt = 0;
+ const int MAX_ATTEMPTS = can_move ? 2 : 1;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && rwalk && dir0 && dir1 && delta > 0);
+ ASSERT(reinject_dir && reinject_dst && reinject_hit);
+
+ if(move_pos) *move_pos = 0;
+
+ do {
+ f2(range, 0, FLT_MAX);
+ fX_set_dX(org, rwalk->vtx.P);
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = delta * 0.01;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, &filter_data, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, &filter_data, &hit1));
+
+ /* Retrieve the medium at the reinjection pos along dir0 */
+ if(SXD_HIT_NONE(&hit0)) {
+ XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir0, (float)delta);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0);
+ if(res == RES_BAD_OP) { mdm0 = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit0.prim.prim_id);
+ side = fX(dot)(dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm0 = interface_get_medium(interf, side);
+ }
+
+ /* Retrieve the medium at the reinjection pos along dir1 */
+ if(SXD_HIT_NONE(&hit1)) {
+ XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir1, (float)delta);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1);
+ if(res == RES_BAD_OP) { mdm1 = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit1.prim.prim_id);
+ side = fX(dot)(dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm1 = interface_get_medium(interf, side);
+ }
+
+ dst0 = dst1 = -1;
+ if(mdm0 == mdm) { /* Check reinjection consistency */
+ if(hit0.distance <= delta_adjusted) {
+ dst0 = hit0.distance;
+ } else {
+ dst0 = delta;
+ hit0 = SXD_HIT_NULL;
+ }
+ }
+ if(mdm1 == mdm) {/* Check reinjection consistency */
+ if(hit1.distance <= delta_adjusted) {
+ dst1 = hit1.distance;
+ } else {
+ dst1 = delta;
+ hit1 = SXD_HIT_NULL;
+ }
+ }
+
+ /* No valid reinjection. Maybe the random walk is near a sharp corner and
+ * thus the ray-tracing misses the enclosure geometry. Another possibility
+ * is that the random walk lies roughly on an edge. In this case, sampled
+ * reinjecton dirs can intersect the primitive on the other side of the
+ * edge. Normally, this primitive should be filtered by the "hit_filter"
+ * function but this may be not the case due to a "threshold effect". In
+ * both situations, try to slightly move away from the primitive boundaries
+ * and retry to find a valid reinjection. */
+ if(dst0 == -1 && dst1 == -1) {
+ XD(move_away_primitive_boundaries)(rwalk, delta);
+ if(move_pos) *move_pos = 1;
+ }
+ } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS);
+
+ if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */
+ log_err(scn->dev, "%s: no valid reinjection direction at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ if(dst0 == -1) {
+ /* Invalid dir0 -> move along dir1 */
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 == -1) {
+ /* Invalid dir1 -> move along dir0 */
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else if(dst0 < reinject_threshold && dst1 < reinject_threshold) {
+ /* The displacement along dir0 and dir1 are both below the reinjection
+ * threshold that defines a distance under which the temperature gradients
+ * are ignored. Move along the direction that allows the maximum
+ * displacement. */
+ if(dst0 > dst1) {
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ }
+ } else if(dst0 < reinject_threshold) {
+ /* Ingore dir0 that is bellow the reinject threshold */
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 < reinject_threshold) {
+ /* Ingore dir1 that is bellow the reinject threshold */
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ /* All reinjection directions are valid. Choose the first 1 that was
+ * randomly selected by the sample_reinjection_dir procedure and adjust
+ * the displacement distance. */
+ dir = dir0;
+ if(dst0 < dst1) {
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dst = dst1;
+ hit = SXD_HIT_NULL;
+ }
+
+ /* If the displacement distance is too close of a boundary, move to the
+ * boundary in order to avoid numerical uncertainty. */
+ if(!SXD_HIT_NONE(&hit0)
+ && dst0 != dst
+ && eq_eps(dst0, dst, dst0*0.1)) {
+ dst = dst0;
+ hit = hit0;
+ }
+ }
+
+ /* Setup output variable */
+ fX(set)(reinject_dir, dir);
+ *reinject_dst = (float)dst;
+ *reinject_hit = hit;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(select_reinjection_dir_and_check_validity)
+ (const struct sdis_scene* scn,
+ const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
+ struct XD(rwalk)* rwalk, /* Current random walk state */
+ const float dir0[DIM], /* Challenged direction */
+ const float dir1[DIM], /* Challanged direction */
+ const double delta, /* Max reinjection distance */
+ float out_reinject_dir[DIM], /* Selected direction */
+ float* out_reinject_dst, /* Effective reinjection distance */
+ int can_move, /* Define of the random wal pos can be moved or not */
+ int* move_pos, /* Define if the current random walk was moved. May be NULL */
+ int* is_valid, /* Define if the reinjection defines a valid pos */
+ struct sXd(hit)* out_reinject_hit) /* Hit along the reinjection dir */
+{
+ double pos[DIM];
+ struct sdis_medium* reinject_mdm;
+ struct sXd(hit) reinject_hit;
+ float reinject_dir[DIM];
+ float reinject_dst;
+ res_T res = RES_OK;
+ ASSERT(is_valid && out_reinject_dir && out_reinject_dst && out_reinject_hit);
+
+ /* Select a reinjection direction */
+ res = XD(select_reinjection_dir)(scn, mdm, rwalk, dir0, dir1, delta,
+ reinject_dir, &reinject_dst, can_move, move_pos, &reinject_hit);
+ if(res != RES_OK) goto error;
+
+ if(!SXD_HIT_NONE(&reinject_hit)) {
+ *is_valid = 1;
+ } else {
+ /* Check medium consistency at the reinjection position */
+ XD(move_pos)(dX(set)(pos, rwalk->vtx.P), reinject_dir, reinject_dst);
+ res = scene_get_medium_in_closed_boundaries
+ (scn, pos, &reinject_mdm);
+ if(res == RES_BAD_OP) { reinject_mdm = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+
+ *is_valid = reinject_mdm == mdm;
+ }
+
+ if(*is_valid) {
+ fX(set)(out_reinject_dir, reinject_dir);
+ *out_reinject_dst = reinject_dst;
+ *out_reinject_hit = reinject_hit;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/* Check that the interface fragment is consistent with the current state of
* the random walk */
static INLINE int
@@ -111,8 +463,9 @@ XD(solid_solid_boundary_path)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
- struct sXd(hit) hit0, hit1, hit2, hit3;
+ struct sXd(hit) hit0, hit1;
struct sXd(hit)* hit;
+ struct XD(rwalk) rwalk_saved;
struct sdis_interface* interf = NULL;
struct sdis_medium* solid_front = NULL;
struct sdis_medium* solid_back = NULL;
@@ -120,18 +473,21 @@ XD(solid_solid_boundary_path)
double lambda_front, lambda_back;
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;
double power;
- float range0[2], range1[2];
float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
+ float dir_front[DIM], dir_back[DIM];
float* dir;
- float pos[DIM];
- int dim = DIM;
+ float reinject_dst_front, reinject_dst_back;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt;
+ int move;
+ int reinjection_is_valid;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -151,30 +507,57 @@ XD(solid_solid_boundary_path)
/* 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_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));
+ delta_boundary_back = delta_back *sqrt(DIM);
+
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
+
+ /* Sample a reinjection direction and reflect it around the normal. Then
+ * reflect them on the back side of the interface. */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ XD(reflect)(dir2, dir0, rwalk->hit.normal);
+ fX(minus)(dir1, dir0);
+ fX(minus)(dir3, dir2);
+
+ /* Select the reinjection direction and distance for the front side */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front, rwalk,
+ dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front, 1, &move,
+ &reinjection_is_valid, &hit0);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
- /* 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);
+ /* Select the reinjection direction and distance for the back side */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_back, rwalk,
+ dir1, dir3, delta_boundary_back, dir_back, &reinject_dst_back, 1, &move,
+ &reinjection_is_valid, &hit1);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
+
+ /* If random walk was moved by the select_reinjection_dir on back side, one
+ * has to rerun the select_reinjection_dir on front side at the new pos */
+ if(move) {
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front,
+ rwalk, dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front,
+ 0, NULL, &reinjection_is_valid, &hit0);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
+ }
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjection */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid soid/solid reinjection at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
/* Define the reinjection side. Note that the proba should be :
* Lf/Df' / (Lf/Df' + Lb/Db')
@@ -190,38 +573,15 @@ XD(solid_solid_boundary_path)
proba = (lambda_front/reinject_dst_front)
/ (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
if(r < proba) { /* Reinject in front */
- dir = dir0;
+ dir = dir_front;
hit = &hit0;
mdm = solid_front;
reinject_dst = reinject_dst_front;
- delta_boundary = delta_boundary_front;
} else { /* Reinject in back */
- dir = dir1;
+ dir = dir_back;
hit = &hit1;
mdm = solid_back;
reinject_dst = reinject_dst_back;
- delta_boundary = delta_boundary_back;
- }
-
- /* 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 */
@@ -229,7 +589,7 @@ XD(solid_solid_boundary_path)
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 = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
@@ -238,9 +598,9 @@ XD(solid_solid_boundary_path)
}
}
- /* Reinject */
+ /* Perform reinjection. */
XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
- if(eq_epsf(hit->distance, (float)reinject_dst, 1.e-4f)) {
+ if(hit->distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
rwalk->hit = *hit;
@@ -278,8 +638,8 @@ XD(solid_fluid_boundary_path)
struct sdis_medium* mdm_back = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
- struct sXd(hit) hit0 = SXD_HIT_NULL;
- struct sXd(hit) hit1 = SXD_HIT_NULL;
+ struct XD(rwalk) rwalk_saved;
+ struct sXd(hit) hit = SXD_HIT_NULL;
struct sdis_interface_fragment frag_fluid;
double hc;
double hr;
@@ -291,10 +651,13 @@ XD(solid_fluid_boundary_path)
double delta_boundary;
double r;
double tmp;
- float pos[DIM];
float dir0[DIM], dir1[DIM];
- float range[2];
- int dim = DIM;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt;
+ int reinjection_is_valid = 0;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -325,48 +688,44 @@ XD(solid_fluid_boundary_path)
* delta. */
delta_boundary = sqrt(DIM) * delta;
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
- if(solid == mdm_back) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, rwalk->hit.normal);
- /* 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;
+ if(solid == mdm_back) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
}
- delta = delta_boundary;
- dim = 1;
+ /* Select the solid reinjection direction and distance */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid, rwalk,
+ dir0, dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
+ &reinjection_is_valid, &hit);
+ if(res != RES_OK) goto error;
+
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjecton */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid solid/fluid reinjection at {%g, %g %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
}
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = reinject_dst / sqrt(DIM);
+
/* Fetch the boundary properties */
epsilon = interface_side_get_emissivity(interf, &frag_fluid);
hc = interface_get_convection_coef(interf, frag);
@@ -393,8 +752,8 @@ XD(solid_fluid_boundary_path)
/* 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 = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ const double delta_in_meter = reinject_dst * fp_to_meter;
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
@@ -403,13 +762,13 @@ XD(solid_fluid_boundary_path)
}
}
- /* Reinject */
- XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
- if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ /* Perform solid reinjection */
+ XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
+ if(hit.distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
- rwalk->hit = hit0;
- rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit = hit;
+ rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
} else {
T->func = XD(conductive_path);
rwalk->mdm = solid;
@@ -440,6 +799,7 @@ XD(solid_boundary_with_flux_path)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
+ struct XD(rwalk) rwalk_saved;
struct sdis_interface* interf = NULL;
struct sdis_medium* mdm = NULL;
double lambda;
@@ -448,13 +808,15 @@ XD(solid_boundary_with_flux_path)
double delta_in_meter;
double power;
double tmp;
- struct sXd(hit) hit0;
- struct sXd(hit) hit1;
- float pos[DIM];
+ struct sXd(hit) hit;
float dir0[DIM];
float dir1[DIM];
- float range[2];
- int dim = DIM;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt = 0;
+ int reinjection_is_valid = 0;
res_T res = RES_OK;
ASSERT(frag && phi != SDIS_FLUX_NONE);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -476,48 +838,43 @@ XD(solid_boundary_with_flux_path)
* 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);
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
+ /* 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;
+ if(frag->side == SDIS_BACK) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
}
- delta = delta_boundary;
- dim = 1;
+ /* Select the reinjection direction and distance */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, mdm, rwalk, dir0,
+ dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
+ &reinjection_is_valid, &hit);
+ if(res != RES_OK) goto error;
+
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjecton */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid solid/fluid with flux reinjection "
+ "at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
}
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = reinject_dst / sqrt(DIM);
+
/* Handle the flux */
delta_in_meter = delta*fp_to_meter;
tmp = delta_in_meter / lambda;
@@ -530,8 +887,8 @@ XD(solid_boundary_with_flux_path)
/* 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 = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ delta_in_meter = reinject_dst * fp_to_meter;
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
res = green_path_add_power_term(ctx->green_path, mdm, &rwalk->vtx, tmp);
@@ -539,13 +896,14 @@ XD(solid_boundary_with_flux_path)
}
}
- /* 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)) {
+ /* Reinject. If the reinjection move the point too close of a boundary,
+ * assume that the zone is isotherm and move to the boundary. */
+ XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
+ if(hit.distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
- rwalk->hit = hit0;
- rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit = hit;
+ rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
} else {
T->func = XD(conductive_path);
rwalk->mdm = mdm;
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -24,6 +24,178 @@
#include "sdis_Xd_begin.h"
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* Sample the next direction to walk toward and compute the distance to travel.
+ * Return the sampled direction `dir0', the distance to travel along this
+ * direction, the hit `hit0' along `dir0' wrt to the returned distance, the
+ * direction `dir1' used to adjust the displacement distance, and the hit
+ * `hit1' along `dir1' used to adjust the displacement distance. */
+static float
+XD(sample_next_step)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const float pos[DIM],
+ const float delta_solid,
+ float dir0[DIM], /* Sampled direction */
+ float dir1[DIM], /* Direction used to adjust delta */
+ struct sXd(hit)* hit0, /* Hit along the sampled direction */
+ struct sXd(hit)* hit1) /* Hit used to adjust delta */
+{
+ struct sXd(hit) hits[2];
+ float dirs[2][DIM];
+ float range[2];
+ float delta;
+ ASSERT(scn && rng && pos && delta_solid>0 && dir0 && dir1 && hit0 && hit1);
+
+ *hit0 = SXD_HIT_NULL;
+ *hit1 = SXD_HIT_NULL;
+
+#if DIM == 2
+ /* Sample a main direction around 2PI */
+ ssp_ran_circle_uniform_float(rng, dirs[0], NULL);
+#else
+ /* Sample a main direction around 4PI */
+ ssp_ran_sphere_uniform_float(rng, dirs[0], NULL);
+#endif
+
+ /* Negate the sampled dir */
+ fX(minus)(dirs[1], dirs[0]);
+
+ /* Use the previously sampled direction to estimate the minimum distance from
+ * `pos' to the scene boundary */
+ f2(range, 0.f, delta_solid*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[0], range, NULL, &hits[0]));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[1], range, NULL, &hits[1]));
+ if(SXD_HIT_NONE(&hits[0]) && SXD_HIT_NONE(&hits[1])) {
+ delta = delta_solid;
+ } else {
+ delta = MMIN(hits[0].distance, hits[1].distance);
+ }
+
+ if(!SXD_HIT_NONE(&hits[0])
+ && delta != hits[0].distance
+ && eq_eps(hits[0].distance, delta, delta_solid*0.1)) {
+ /* Set delta to the main hit distance if it is roughly equal to it in order
+ * to avoid numerical issues on moving along the main direction. */
+ delta = hits[0].distance;
+ }
+
+ /* Setup outputs */
+ if(delta <= delta_solid*0.1 && hits[1].distance == delta) {
+ /* Snap the random walk to the boundary if delta is too small */
+ fX(set)(dir0, dirs[1]);
+ *hit0 = hits[1];
+ fX(splat)(dir1, (float)INF);
+ *hit1 = SXD_HIT_NULL;
+ } else {
+ fX(set)(dir0, dirs[0]);
+ *hit0 = hits[0];
+ if(delta == hits[0].distance) {
+ fX(set)(dir1, dirs[0]);
+ *hit1 = hits[0];
+ } else if(delta == hits[1].distance) {
+ fX(set)(dir1, dirs[1]);
+ *hit1 = hits[1];
+ } else {
+ fX(splat)(dir1, 0);
+ *hit1 = SXD_HIT_NULL;
+ }
+ }
+
+ return delta;
+}
+
+/* Sample the next direction to walk toward and compute the distance to travel.
+ * If the targeted position does not lie inside the current medium, reject it
+ * and sample a new next step. */
+static res_T
+XD(sample_next_step_robust)
+ (struct sdis_scene* scn,
+ struct sdis_medium* current_mdm,
+ struct ssp_rng* rng,
+ const double pos[DIM],
+ const float delta_solid,
+ float dir0[DIM], /* Sampled direction */
+ float dir1[DIM], /* Direction used to adjust delta */
+ struct sXd(hit)* hit0, /* Hit along the sampled direction */
+ struct sXd(hit)* hit1, /* Hit used to adjust delta */
+ float* out_delta)
+{
+ struct sdis_medium* mdm;
+ float delta;
+ float org[DIM];
+ const size_t MAX_ATTEMPTS = 100;
+ size_t iattempt = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && current_mdm && rng && pos && delta_solid > 0);
+ ASSERT(dir0 && dir1 && hit0 && hit1 && out_delta);
+
+ fX_set_dX(org, pos);
+ do {
+ double pos_next[DIM];
+
+ /* Compute the next step */
+ delta = XD(sample_next_step)
+ (scn, rng, org, delta_solid, dir0, dir1, hit0, hit1);
+
+ /* Retrieve the medium of the next step */
+ if(hit0->distance > delta) {
+ XD(move_pos)(dX(set)(pos_next, pos), dir0, delta);
+ res = scene_get_medium_in_closed_boundaries(scn, pos_next, &mdm);
+ if(res == RES_BAD_OP) { mdm = NULL; res = RES_OK; }
+ if(res != RES_OK) goto error;
+ } else {
+ struct sdis_interface* interf;
+ enum sdis_side side;
+ interf = scene_get_interface(scn, hit0->prim.prim_id);
+ side = fX(dot)(dir0, hit0->normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm = interface_get_medium(interf, side);
+ }
+
+ /* Check medium consistency */
+ if(current_mdm != mdm) {
+#if 0
+#if DIM == 2
+ log_err(scn->dev,
+ "%s: inconsistent medium during the solid random walk at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev,
+ "%s: inconsistent medium during the solid random walk at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+#endif
+ }
+ } while(current_mdm != mdm && ++iattempt < MAX_ATTEMPTS);
+
+ /* Handle error */
+ if(iattempt >= MAX_ATTEMPTS) {
+#if DIM == 2
+ log_err(scn->dev,
+ "%s: could not find a next valid conductive step at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev,
+ "%s: could not find a next valid conductive step at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ *out_delta = delta;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
res_T
XD(conductive_path)
(struct sdis_scene* scn,
@@ -44,8 +216,8 @@ XD(conductive_path)
(void)ctx, (void)istep;
/* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK);
- if(mdm != rwalk->mdm) {
+ res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm);
+ if(res != RES_OK || 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));
res = RES_BAD_OP_IRRECOVERABLE;
@@ -61,7 +233,6 @@ XD(conductive_path)
}
do { /* Solid random walk */
- struct get_medium_info info = GET_MEDIUM_INFO_NULL;
struct sXd(hit) hit0, hit1;
double lambda; /* Thermal conductivity */
double rho; /* Volumic mass */
@@ -70,7 +241,6 @@ XD(conductive_path)
double power_factor = 0;
double power;
float delta, delta_solid; /* Random walk numerical parameter */
- float range[2];
float dir0[DIM], dir1[DIM];
float org[DIM];
@@ -109,39 +279,20 @@ XD(conductive_path)
goto error;
}
-#if DIM == 2
- /* Sample a direction around 2PI */
- ssp_ran_circle_uniform_float(rng, dir0, NULL);
-#else
- /* Sample a direction around 4PI */
- ssp_ran_sphere_uniform_float(rng, dir0, NULL);
-#endif
-
- /* Trace a ray along the sampled direction and its opposite to check if a
- * surface is hit in [0, delta_solid]. */
fX_set_dX(org, rwalk->vtx.P);
- fX(minus)(dir1, dir0);
- hit0 = hit1 = SXD_HIT_NULL;
- range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
-
- 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) {
+
+ /* Sample the direction to walk toward and compute the distance to travel */
+ res = XD(sample_next_step_robust)(scn, mdm, rng, rwalk->vtx.P, delta_solid,
+ dir0, dir1, &hit0, &hit1, &delta);
+ if(res != RES_OK) goto error;
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ if((S3D_HIT_NONE(&hit0) && S3D_HIT_NONE(&hit1))) { /* Hit nothing */
const double delta_in_meter = delta * fp_to_meter;
power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * power_factor;
- }
- } 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) {
+ } else {
const double delta_s_adjusted = delta_solid * RAY_RANGE_MAX_SCALE;
const double delta_s_in_meter = delta_solid * fp_to_meter;
double h;
@@ -226,10 +377,8 @@ XD(conductive_path)
}
}
- /* Define if the random walk hits something along dir0. Multiply delta by
- * the empirical ray range scale factor to ensure that once moved, the
- * random walk does not lie in the uncertainty zone near the geometry */
- if(hit0.distance > delta * RAY_RANGE_MAX_SCALE) {
+ /* Define if the random walk hits something along dir0 */
+ if(hit0.distance > delta) {
rwalk->hit = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
} else {
@@ -245,45 +394,6 @@ XD(conductive_path)
(ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
if(res != RES_OK) goto error;
- /* Fetch the current medium */
- if(SXD_HIT_NONE(&rwalk->hit)) {
- 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);
- mdm = interface_get_medium(interf, rwalk->hit_side);
- }
-
- /* Check random walk consistency */
- if(mdm != rwalk->mdm) {
- log_err(scn->dev,
- "%s: inconsistent medium during the solid random walk.\n", FUNC_NAME);
-#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
- res = RES_BAD_OP;
- goto error;
- }
-
++istep;
/* Keep going while the solid random walk does not hit an interface */
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -34,6 +34,7 @@ XD(register_heat_vertex_in_fluid)
const double weight)
{
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ struct hit_filter_data filter_data;
const float empirical_dst = 0.1f;
const float range[2] = {0, FLT_MAX};
float org[DIM];
@@ -50,7 +51,9 @@ XD(register_heat_vertex_in_fluid)
fX(set)(dir, rwalk->hit.normal);
if(rwalk->hit_side == SDIS_BACK) fX(minus)(dir, dir);
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &rwalk->hit, &hit));
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = 1.e-6;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &filter_data, &hit));
dst = SXD_HIT_NONE(&hit) ? empirical_dst : hit.distance * 0.5f;
vtx = rwalk->vtx;
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -52,6 +52,7 @@ XD(trace_radiative_path)
/* Launch the radiative random walk */
for(;;) {
const struct sdis_interface* interf = NULL;
+ struct hit_filter_data filter_data;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
struct sdis_medium* chk_mdm = NULL;
double alpha;
@@ -63,12 +64,14 @@ XD(trace_radiative_path)
fX_set_dX(pos, rwalk->vtx.P);
/* Trace the radiative ray */
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = 1.e-6;
#if (SDIS_XD_DIMENSION == 2)
SXD(scene_view_trace_ray_3d
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
#else
SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
#endif
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
rwalk->hit_side = SDIS_SIDE_NULL__;
diff --git a/src/sdis_misc.h b/src/sdis_misc.h
@@ -20,6 +20,14 @@
#include <rsys/float3.h>
#include <star/ssp.h>
+struct accum {
+ double sum; /* Sum of MC weights */
+ double sum2; /* Sum of square MC weights */
+ size_t count; /* #accumulated MC weights */
+};
+#define ACCUM_NULL__ {0,0,0}
+static const struct accum ACCUM_NULL = ACCUM_NULL__;
+
/* 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.001f
@@ -31,6 +39,24 @@
#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
+static INLINE void
+sum_accums
+ (const struct accum accums[],
+ const size_t naccums,
+ struct accum* accum)
+{
+ struct accum acc = ACCUM_NULL;
+ size_t i;
+ ASSERT(accums && naccums && accum);
+
+ FOR_EACH(i, 0, naccums) {
+ acc.sum += accums[i].sum;
+ acc.sum2 += accums[i].sum2;
+ acc.count += accums[i].count;
+ }
+ *accum = acc;
+}
+
/* Reflect the V wrt the normal N. By convention V points outward the surface */
static FINLINE float*
reflect_2d(float res[2], const float V[2], const float N[2])
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -32,6 +32,7 @@ ray_realisation_3d
const double fp_to_meter,
const double Tarad,
const double Tref,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -50,9 +51,14 @@ ray_realisation_3d
ctx.Tarad = Tarad;
ctx.Tref3 = Tref*Tref*Tref;
+ ctx.heat_path = heat_path;
f3_set_d3(dir, direction);
+ /* Register the starting position against the heat path */
+ res = register_heat_vertex(heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE);
+ if(res != RES_OK) goto error;
+
res = trace_radiative_path_3d(scn, dir, fp_to_meter, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -46,8 +46,8 @@ probe_realisation_2d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -61,8 +61,8 @@ probe_realisation_3d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
/*******************************************************************************
@@ -79,8 +79,8 @@ boundary_realisation_2d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -94,8 +94,8 @@ boundary_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -140,6 +140,7 @@ ray_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
#endif /* SDIS_REALISATION_H */
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -45,9 +45,9 @@ XD(compute_temperature)
}* stack = NULL;
size_t istack = 0;
#endif
- /* Maximum accepted #failures before stopping the realisation */
struct sdis_heat_vertex* heat_vtx = NULL;
- const size_t MAX_FAILS = 10;
+ /* Maximum accepted #failures before stopping the realisation */
+ const size_t MAX_FAILS = 1;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
@@ -121,8 +121,8 @@ XD(probe_realisation)
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -213,8 +213,8 @@ XD(boundary_realisation)
const double fp_to_meter,
const double Tarad,
const double Tref,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -373,3 +373,14 @@ scene_get_medium
: scene_get_medium_3d(scn, pos, info, out_medium);
}
+res_T
+scene_get_medium_in_closed_boundaries
+ (const struct sdis_scene* scn,
+ const double pos[],
+ struct sdis_medium** out_medium)
+{
+ return scene_is_2d(scn)
+ ? scene_get_medium_in_closed_boundaries_2d(scn, pos, out_medium)
+ : scene_get_medium_in_closed_boundaries_3d(scn, pos, out_medium);
+}
+
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -22,8 +22,14 @@
#include "sdis_medium_c.h"
#include "sdis_scene_c.h"
+#include <star/ssp.h>
+#include <rsys/float22.h>
+#include <rsys/float33.h>
#include <rsys/rsys.h>
+/* Emperical cos threshold defining if an angle is sharp */
+#define SHARP_ANGLE_COS_THRESOLD -0.70710678 /* ~ cos(3*PI/4) */
+
/*******************************************************************************
* Define the helper functions and the data types used by the scene
* independently of its dimension, i.e. 2D or 3D.
@@ -148,6 +154,7 @@ clear_properties(struct sdis_scene* scn)
/* Vector macros generic to SDIS_SCENE_DIMENSION */
#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+#define fXX_mulfX CONCAT(CONCAT(CONCAT(CONCAT(f, DIM), DIM), _mulf), DIM)
/* Macro making generic its subimitted name to SDIS_SCENE_DIMENSION */
#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
@@ -162,40 +169,261 @@ clear_properties(struct sdis_scene* scn)
* Helper functions
******************************************************************************/
#if DIM == 2
-/* Check that `hit' roughly lies on a vertex. For segments, a simple but
- * approximative way is to test that its position have at least one barycentric
- * coordinate roughly equal to 0 or 1. */
-static FINLINE int
-hit_on_vertex(const struct s2d_hit* hit)
+#define ON_VERTEX_EPSILON 1.e-4f
+/* Check that `hit' roughly lies on a vertex. */
+static INLINE int
+hit_on_vertex
+ (const struct s2d_hit* hit,
+ const float org[3],
+ const float dir[3])
+{
+ struct s2d_attrib v0, v1;
+ float E[2];
+ float hit_pos[2];
+ float segment_len;
+ float hit_len0;
+ float hit_len1;
+ ASSERT(hit && !S2D_HIT_NONE(hit) && org && dir);
+
+ /* Rertieve the segment vertices */
+ S2D(segment_get_vertex_attrib(&hit->prim, 0, S2D_POSITION, &v0));
+ S2D(segment_get_vertex_attrib(&hit->prim, 1, S2D_POSITION, &v1));
+
+ /* Compute the length of the segment */
+ segment_len = f2_len(f2_sub(E, v1.value, v0.value));
+
+ /* Compute the hit position onto the segment */
+ f2_add(hit_pos, org, f2_mulf(hit_pos, dir, hit->distance));
+
+ /* Compute the length from hit position to segment vertices */
+ hit_len0 = f2_len(f2_sub(E, v0.value, hit_pos));
+ hit_len1 = f2_len(f2_sub(E, v1.value, hit_pos));
+
+ if(hit_len0 / segment_len < ON_VERTEX_EPSILON
+ || hit_len1 / segment_len < ON_VERTEX_EPSILON)
+ return 1;
+ return 0;
+}
+
+static int
+hit_shared_vertex
+ (const struct s2d_primitive* seg0,
+ const struct s2d_primitive* seg1,
+ const float pos0[2], /* Tested position onto the segment 0 */
+ const float pos1[2]) /* Tested Position onto the segment 1 */
{
- const float on_vertex_eps = 1.e-4f;
- float v;
- ASSERT(hit && !S2D_HIT_NONE(hit));
- v = 1.f - hit->u;
- return eq_epsf(hit->u, 0.f, on_vertex_eps)
- || eq_epsf(hit->u, 1.f, on_vertex_eps)
- || eq_epsf(v, 0.f, on_vertex_eps)
- || eq_epsf(v, 1.f, on_vertex_eps);
+ struct s2d_attrib seg0_vertices[2]; /* Vertex positions of the segment 0 */
+ struct s2d_attrib seg1_vertices[2]; /* Vertex positions of the segment 1 */
+ float d0[2], d1[2]; /* temporary vector */
+ float seg0_len, seg1_len; /* Length of the segments */
+ float tmp0_len, tmp1_len;
+ float cos_normals;
+ int seg0_vert = -1; /* Id of the shared vertex for the segment 0 */
+ int seg1_vert = -1; /* Id of the shared vertex for the segment 1 */
+ int seg0_ivertex, seg1_ivertex;
+ ASSERT(seg0 && seg1 && pos0 && pos1);
+
+ /* Fetch the vertices of the segment 0 */
+ S2D(segment_get_vertex_attrib(seg0, 0, S2D_POSITION, &seg0_vertices[0]));
+ S2D(segment_get_vertex_attrib(seg0, 1, S2D_POSITION, &seg0_vertices[1]));
+
+ /* Fetch the vertices of the segment 1 */
+ S2D(segment_get_vertex_attrib(seg1, 0, S2D_POSITION, &seg1_vertices[0]));
+ S2D(segment_get_vertex_attrib(seg1, 1, S2D_POSITION, &seg1_vertices[1]));
+
+ /* Look for the vertex shared by the 2 segments */
+ for(seg0_ivertex = 0; seg0_ivertex < 2 && seg0_vert < 0; ++seg0_ivertex) {
+ for(seg1_ivertex = 0; seg1_ivertex < 2 && seg1_vert < 0; ++seg1_ivertex) {
+ const int vertex_eq = f2_eq_eps
+ (seg0_vertices[seg0_ivertex].value,
+ seg1_vertices[seg1_ivertex].value,
+ 1.e-6f);
+ if(vertex_eq) {
+ seg0_vert = seg0_ivertex;
+ seg1_vert = seg1_ivertex;
+ /* We assume that the segments are not degenerated. As a consequence we
+ * can break here since a vertex of the segment 0 can be equal to at most
+ * one vertex of the segment 1 */
+ break;
+ }
+ }}
+
+ /* The segments do not have a common vertex */
+ if(seg0_vert < 0) return 0;
+
+ /* Compute the dirctions from shared vertex to the opposite segment vertex */
+ f2_sub(d0, seg0_vertices[(seg0_vert+1)%2].value, seg0_vertices[seg0_vert].value);
+ f2_sub(d1, seg1_vertices[(seg1_vert+1)%2].value, seg1_vertices[seg1_vert].value);
+
+ /* Compute the cosine between the segments */
+ seg0_len = f2_normalize(d0, d0);
+ seg1_len = f2_normalize(d1, d1);
+ cos_normals = f2_dot(d0, d1);
+
+ /* The angle formed by the 2 segments is sharp. Do not filter the hit */
+ if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+
+ /* Compute the length from pos<0|1> to shared vertex */
+ f2_sub(d0, seg0_vertices[seg0_vert].value, pos0);
+ f2_sub(d1, seg1_vertices[seg1_vert].value, pos1);
+ tmp0_len = f2_len(d0);
+ tmp1_len = f2_len(d1);
+
+ return (eq_epsf(seg0_len, 0, 1.e-6f) || tmp0_len/seg0_len < ON_VERTEX_EPSILON)
+ && (eq_epsf(seg1_len, 0, 1.e-6f) || tmp1_len/seg1_len < ON_VERTEX_EPSILON);
}
#else /* DIM == 3 */
-/* Check that `hit' roughly lies on an edge. For triangular primitives, a
- * simple but approximative way is to test that its position have at least one
- * barycentric coordinate roughly equal to 0 or 1. */
-static FINLINE int
-hit_on_edge(const struct s3d_hit* hit)
+#define ON_EDGE_EPSILON 1.e-4f
+/* Check that `hit' roughly lies on an edge. */
+static INLINE int
+hit_on_edge
+ (const struct s3d_hit* hit,
+ const float org[3],
+ const float dir[3])
{
- const float on_edge_eps = 1.e-4f;
- float w;
- ASSERT(hit && !S3D_HIT_NONE(hit));
- w = 1.f - hit->uv[0] - hit->uv[1];
- return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
- || eq_epsf(w, 0.f, on_edge_eps)
- || eq_epsf(w, 1.f, on_edge_eps);
+ struct s3d_attrib v0, v1, v2;
+ float E0[3], E1[3], N[3];
+ float tri_2area;
+ float hit_2area0;
+ float hit_2area1;
+ float hit_2area2;
+ float hit_pos[3];
+ ASSERT(hit && !S3D_HIT_NONE(hit) && org && dir);
+
+ /* Retrieve the triangle vertices */
+ S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the triangle area * 2 */
+ f3_sub(E0, v1.value, v0.value);
+ f3_sub(E1, v2.value, v0.value);
+ tri_2area = f3_len(f3_cross(N, E0, E1));
+
+ /* Compute the hit position */
+ f3_add(hit_pos, org, f3_mulf(hit_pos, dir, hit->distance));
+
+ /* Compute areas */
+ f3_sub(E0, v0.value, hit_pos);
+ f3_sub(E1, v1.value, hit_pos);
+ hit_2area0 = f3_len(f3_cross(N, E0, E1));
+ f3_sub(E0, v1.value, hit_pos);
+ f3_sub(E1, v2.value, hit_pos);
+ hit_2area1 = f3_len(f3_cross(N, E0, E1));
+ f3_sub(E0, v2.value, hit_pos);
+ f3_sub(E1, v0.value, hit_pos);
+ hit_2area2 = f3_len(f3_cross(N, E0, E1));
+
+ if(hit_2area0 / tri_2area < ON_EDGE_EPSILON
+ || hit_2area1 / tri_2area < ON_EDGE_EPSILON
+ || hit_2area2 / tri_2area < ON_EDGE_EPSILON)
+ return 1;
+
+ return 0;
}
+
+static int
+hit_shared_edge
+ (const struct s3d_primitive* tri0,
+ const struct s3d_primitive* tri1,
+ const float pos0[3], /* Tested position onto the triangle 0 */
+ const float pos1[3]) /* Tested Position onto the triangle 1 */
+{
+ struct s3d_attrib tri0_vertices[3]; /* Vertex positions of the triangle 0 */
+ struct s3d_attrib tri1_vertices[3]; /* Vertex positions of the triangle 1 */
+ float E0[3], E1[3]; /* Temporary variables storing triangle edges */
+ float N0[3], N1[3]; /* Temporary Normals */
+ float tri0_2area, tri1_2area; /* 2*area of the submitted triangles */
+ float tmp0_2area, tmp1_2area;
+ float cos_normals;
+ int tri0_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 0 */
+ int tri1_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 1 */
+ int edge_ivertex = 0; /* Temporary variable */
+ int tri0_ivertex, tri1_ivertex;
+ int iv0, iv1, iv2;
+ ASSERT(tri0 && tri1 && pos0 && pos1);
+
+ /* Fetch the vertices of the triangle 0 */
+ S3D(triangle_get_vertex_attrib(tri0, 0, S3D_POSITION, &tri0_vertices[0]));
+ S3D(triangle_get_vertex_attrib(tri0, 1, S3D_POSITION, &tri0_vertices[1]));
+ S3D(triangle_get_vertex_attrib(tri0, 2, S3D_POSITION, &tri0_vertices[2]));
+
+ /* Fetch the vertices of the triangle 1 */
+ S3D(triangle_get_vertex_attrib(tri1, 0, S3D_POSITION, &tri1_vertices[0]));
+ S3D(triangle_get_vertex_attrib(tri1, 1, S3D_POSITION, &tri1_vertices[1]));
+ S3D(triangle_get_vertex_attrib(tri1, 2, S3D_POSITION, &tri1_vertices[2]));
+
+ /* Look for the vertices shared by the 2 triangles */
+ for(tri0_ivertex=0; tri0_ivertex < 3 && edge_ivertex < 2; ++tri0_ivertex) {
+ for(tri1_ivertex=0; tri1_ivertex < 3 && edge_ivertex < 2; ++tri1_ivertex) {
+ const int vertex_eq = f3_eq_eps
+ (tri0_vertices[tri0_ivertex].value,
+ tri1_vertices[tri1_ivertex].value,
+ 1.e-6f);
+ if(vertex_eq) {
+ tri0_edge[edge_ivertex] = tri0_ivertex;
+ tri1_edge[edge_ivertex] = tri1_ivertex;
+ ++edge_ivertex;
+ /* We assume that the triangles are not degenerated. As a consequence we
+ * can break here since a vertex of the triangle 0 can be equal to at
+ * most one vertex of the triangle 1 */
+ break;
+ }
+ }}
+
+ /* The triangles do not have a common edge */
+ if(edge_ivertex < 2) return 0;
+
+ /* Ensure that the vertices of the shared edge are registered in the right
+ * order regarding the triangle vertices, i.e. (0,1), (1,2) or (2,0) */
+ if((tri0_edge[0]+1)%3 != tri0_edge[1]) SWAP(int, tri0_edge[0], tri0_edge[1]);
+ if((tri1_edge[0]+1)%3 != tri1_edge[1]) SWAP(int, tri1_edge[0], tri1_edge[1]);
+
+ /* Compute the shared edge normal lying in the triangle 0 plane */
+ iv0 = tri0_edge[0];
+ iv1 = tri0_edge[1];
+ iv2 = (tri0_edge[1]+1) % 3;
+ f3_sub(E0, tri0_vertices[iv1].value, tri0_vertices[iv0].value);
+ f3_sub(E1, tri0_vertices[iv2].value, tri0_vertices[iv0].value);
+ f3_cross(N0, E0, E1); /* Triangle 0 normal */
+ tri0_2area = f3_len(N0);
+ f3_cross(N0, N0, E0);
+
+ /* Compute the shared edge normal lying in the triangle 1 plane */
+ iv0 = tri1_edge[0];
+ iv1 = tri1_edge[1];
+ iv2 = (tri1_edge[1]+1) % 3;
+ f3_sub(E0, tri1_vertices[iv1].value, tri1_vertices[iv0].value);
+ f3_sub(E1, tri1_vertices[iv2].value, tri1_vertices[iv0].value);
+ f3_cross(N1, E0, E1);
+ tri1_2area = f3_len(N1);
+ f3_cross(N1, N1, E0);
+
+ /* Compute the cosine between the 2 edge normals */
+ f3_normalize(N0, N0);
+ f3_normalize(N1, N1);
+ cos_normals = f3_dot(N0, N1);
+
+ /* The angle formed by the 2 triangles is sharp */
+ if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+
+ /* Compute the 2 times the area of the (pos0, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri0_vertices[tri0_edge[0]].value, pos0);
+ f3_sub(E1, tri0_vertices[tri0_edge[1]].value, pos0);
+ tmp0_2area = f3_len(f3_cross(N0, E0, E1));
+
+ /* Compute the 2 times the area of the (pos1, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri1_vertices[tri1_edge[0]].value, pos1);
+ f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1);
+ tmp1_2area = f3_len(f3_cross(N1, E0, E1));
+
+ return (eq_epsf(tri0_2area, 0, 1.e-6f) || tmp0_2area/tri0_2area < ON_EDGE_EPSILON)
+ && (eq_epsf(tri1_2area, 0, 1.e-6f) || tmp1_2area/tri1_2area < ON_EDGE_EPSILON);
+}
+#undef ON_EDGE_EPSILON
#endif /* DIM == 2 */
/* Avoid self-intersection for a ray starting from a planar primitive, i.e. a
@@ -206,22 +434,25 @@ XD(hit_filter_function)
const float org[DIM],
const float dir[DIM],
void* ray_data,
- void* filter_data)
+ void* global_data)
{
- const struct sXd(hit)* hit_from = ray_data;
- (void)org, (void)dir, (void)filter_data;
+ const struct hit_filter_data* filter_data = ray_data;
+ const struct sXd(hit)* hit_from = &filter_data->XD(hit);
+ (void)org, (void)dir, (void)global_data;
- if(!hit_from || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */
+ if(!ray_data || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */
if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
- if(eq_epsf(hit->distance, 0, 1.e-6f)) {
+ if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) {
+ float pos[DIM];
+ fX(add)(pos, org, fX(mulf)(pos, dir, hit->distance));
/* If the targeted point is near of the origin, check that it lies on an
- * edge/vertex shared by the 2 primitives. */
+ * edge/vertex shared by the 2 primitives */
#if DIM == 2
- return hit_on_vertex(hit_from) && hit_on_vertex(hit);
+ return hit_shared_vertex(&hit_from->prim, &hit->prim, org, pos);
#else
- return hit_on_edge(hit_from) && hit_on_edge(hit);
+ return hit_shared_edge(&hit_from->prim, &hit->prim, org, pos);
#endif
}
return 0;
@@ -321,7 +552,7 @@ XD(run_analyze)
ASSERT(scn && nprims && indices && interf && nverts && position && out_desc);
res = sencXd(device_create)(scn->dev->logger, scn->dev->allocator,
- scn->dev->nthreads, scn->dev->verbose, &senc);
+ 1/*scn->dev->nthreads*/, scn->dev->verbose, &senc);
if(res != RES_OK) goto error;
res = sencXd(scene_create)(senc,
@@ -801,6 +1032,7 @@ XD(scene_get_medium)
size_t iprim, nprims;
size_t nfailures = 0;
const size_t max_failures = 10;
+ float P[DIM];
/* Range of the parametric coordinate into which positions are challenged */
#if DIM == 2
float st[3];
@@ -821,23 +1053,37 @@ XD(scene_get_medium)
f2(st[2], 5.f/12.f, 5.f/12.f);
#endif
+ fX_set_dX(P, pos);
+
SXD(scene_view_primitives_count(scn->sXd(view), &nprims));
FOR_EACH(iprim, 0, nprims) {
struct sXd(hit) hit;
struct sXd(attrib) attr;
struct sXd(primitive) prim;
+ size_t iprim2;
const float range[2] = {0.f, FLT_MAX};
- float N[DIM], P[DIM], dir[DIM], cos_N_dir;
+ float N[DIM], dir[DIM], cos_N_dir;
size_t istep = 0;
+ /* 1 primitive over 2, take a primitive from the end of the primitive list.
+ * When primitives are sorted in a coherent manner regarding their
+ * position, this strategy avoids to test primitives that are going to be
+ * rejected of the same manner due to possible numerical issues of the
+ * resulting intersection. */
+ if((iprim % 2) == 0) {
+ iprim2 = iprim / 2;
+ } else {
+ iprim2 = nprims - 1 - (iprim / 2);
+ }
+
do {
/* Retrieve a position onto the primitive */
- SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim));
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim2, &prim));
SXD(primitive_get_attrib(&prim, SXD_POSITION, st[istep], &attr));
/* Trace a ray from the random walk vertex toward the retrieved primitive
* position */
- fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos)));
+ fX(normalize)(dir, fX(sub)(dir, attr.value, P));
SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit));
/* Unforeseen error. One has to intersect a primitive ! */
@@ -850,20 +1096,20 @@ XD(scene_get_medium)
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((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit))
+ /* Discard the hit if it is on a vertex/edge, and target a new position
+ * onto the current primitive */
+ } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit, P, dir))
&& ++istep < nsteps);
/* The hits of all targeted positions on the current primitive are on
* vertices. Challenge positions on another primitive. */
- if(istep > nsteps) continue;
+ if(istep >= nsteps) continue;
fX(normalize)(N, hit.normal);
cos_N_dir = fX(dot)(N, dir);
/* Not too close and not roughly orthognonal */
- if(hit.distance > 1.e-6 || absf(cos_N_dir) > 1.e-1f) {
+ if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) {
const struct sdis_interface* interf;
interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
@@ -880,6 +1126,94 @@ XD(scene_get_medium)
}
}
+ if(iprim >= nprims) {
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ if(iprim > 10 && iprim > (size_t)((double)nprims * 0.05)) {
+ log_warn(scn->dev,
+ "%s: performance issue. Up to %lu primitives were tested to define the "
+ "current medium at {%g, %g, %g}.\n",
+ FUNC_NAME, (unsigned long)iprim, SPLIT3(P));
+ }
+
+exit:
+ *out_medium = medium;
+ return res;
+error:
+#if DIM == 2
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+ goto exit;
+}
+
+static INLINE res_T
+XD(scene_get_medium_in_closed_boundaries)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ struct sdis_medium** out_medium)
+{
+ struct sdis_medium* medium = NULL;
+ float P[DIM];
+ float frame[DIM*DIM];
+ float dirs[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
+ int idir;
+ res_T res = RES_OK;
+ ASSERT(scn && pos);
+
+ /* Build a frame that will be used to rotate the main axis by PI/4 around
+ * each axis. This can avoid numerical issues when geometry is discretized
+ * along the main axis */
+#if DIM == 2
+ f22_rotation(frame, (float)PI/4);
+#else
+/* N[0] = N[1] = N[2] = (float)(1.0 / sqrt(3.0));*/
+/* f33_basis(frame, N);*/
+ f33_rotation(frame, (float)PI/4, (float)PI/4, (float)PI/4);
+#endif
+
+ fX_set_dX(P, pos);
+ FOR_EACH(idir, 0, 2*DIM) {
+ struct sXd(hit) hit;
+ float N[DIM];
+ const float range[2] = {0.f, FLT_MAX};
+ float cos_N_dir;
+
+ /* Transform the directions to avoid to be aligned with the axis */
+ fXX_mulfX(dirs[idir], frame, dirs[idir]);
+
+ /* Trace a ray from the random walk vertex toward the retrieved primitive
+ * position */
+ SXD(scene_view_trace_ray(scn->sXd(view), P, dirs[idir], range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(SXD_HIT_NONE(&hit)) continue;
+
+ /* Discard a hits if it lies on an edge/point */
+ if(HIT_ON_BOUNDARY(&hit, P, dirs[idir])) continue;
+
+ fX(normalize)(N, hit.normal);
+ cos_N_dir = fX(dot)(N, dirs[idir]);
+
+ /* Not too close and not roughly orthogonal */
+ if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) {
+ 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);
+ break;
+ }
+ }
+ if(idir >= 2*DIM) {
+ res = XD(scene_get_medium)(scn, pos, NULL, &medium);
+ if(res != RES_OK) goto error;
+ }
+
exit:
*out_medium = medium;
return res;
@@ -893,6 +1227,7 @@ error:
#endif
goto exit;
}
+
#undef SDIS_SCENE_DIMENSION
#undef DIM
#undef sencXd
@@ -909,6 +1244,7 @@ error:
#undef SXD_PRIMITIVE_EQ
#undef fX
#undef fX_set_dX
+#undef fXX_mulfX
#undef XD
#undef HIT_ON_BOUNDARY
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -31,6 +31,12 @@ struct prim_prop {
unsigned back_enclosure; /* Id of the back facing enclosure */
};
+struct hit_filter_data {
+ struct s2d_hit hit_2d;
+ struct s3d_hit hit_3d;
+ double epsilon; /* Threshold defining roughly equal intersections */
+};
+
struct get_medium_info {
/* Targeted position */
float pos_tgt[3];
@@ -93,6 +99,7 @@ enclosure_init(struct mem_allocator* allocator, struct enclosure* enc)
enc->S_over_V = 0;
enc->V = 0;
enc->hc_upper_bound = 0;
+ enc->medium_id = UINT_MAX;
}
static INLINE void
@@ -117,6 +124,7 @@ enclosure_copy(struct enclosure* dst, const struct enclosure* src)
dst->S_over_V = src->S_over_V;
dst->V = src->V;
dst->hc_upper_bound = src->hc_upper_bound;
+ dst->medium_id = src->medium_id;
return darray_uint_copy(&dst->local2global, &src->local2global);
}
@@ -139,6 +147,7 @@ enclosure_copy_and_release(struct enclosure* dst, struct enclosure* src)
dst->S_over_V = src->S_over_V;
dst->V = src->V;
dst->hc_upper_bound = src->hc_upper_bound;
+ dst->medium_id = src->medium_id;
return RES_OK;
}
@@ -223,6 +232,22 @@ scene_get_medium
struct get_medium_info* info, /* May be NULL */
struct sdis_medium** medium);
+/* This function assumes that the tested position lies into finite enclosure.
+ * The medium into which it lies is thus retrieved by tracing a random ray
+ * around the current position. For possible infinite enclosure, one has to use
+ * the `scene_get_medium' function instead that, in counterpart, can be more
+ * time consuming.
+ *
+ * Note that actually, the function internally calls scene_get_medium if no
+ * valid medium is found with the regular procedure. This may be due to
+ * numerical issues or wrong assumptions on the current medium (its boundaries
+ * are opened to infinity). */
+extern LOCAL_SYM res_T
+scene_get_medium_in_closed_boundaries
+ (const struct sdis_scene* scn,
+ const double position[],
+ struct sdis_medium** medium);
+
static INLINE void
scene_get_enclosure_ids
(const struct sdis_scene* scn,
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -17,6 +17,7 @@
#include "sdis_camera.h"
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
#include "sdis_interface_c.h"
#include <star/ssp.h>
@@ -66,28 +67,43 @@ solve_pixel
struct ssp_rng* rng,
struct sdis_medium* mdm,
const struct sdis_camera* cam,
- const double time, /* Observation time */
+ const double time_range[2], /* Observation time */
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
const size_t ipix[2], /* Pixel coordinate in the image plane */
const size_t nrealisations,
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accum)
+ struct sdis_estimator* estimator)
{
- double sum_weights = 0;
- double sum_weights_sqr = 0;
- size_t N = 0; /* #realisations that do not fail */
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
size_t irealisation;
res_T res = RES_OK;
ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(estimator && time_range);
FOR_EACH(irealisation, 0, nrealisations) {
+ struct time t0, t1;
double samp[2]; /* Pixel sample */
double ray_pos[3];
double ray_dir[3];
double w = 0;
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
+ double time;
+ res_T res_simul = RES_OK;
+
+ /* Begin time registration */
+ time_current(&t0);
+
+ time = sample_time(rng, time_range);
+ if(register_paths) {
+ heat_path_init(scn->dev->allocator, &heat_path);
+ pheat_path = &heat_path;
+ }
/* Generate a sample into the pixel to estimate */
samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
@@ -98,23 +114,45 @@ solve_pixel
camera_ray(cam, samp, ray_pos, ray_dir);
/* Launch the realisation */
- res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
- time, fp_to_meter, Tarad, Tref, &w);
- if(res == RES_OK) {
- sum_weights += w;
- sum_weights_sqr += w*w;
- ++N;
- } else if(res == RES_BAD_OP) {
- res = RES_OK;
- } else {
+ res_simul = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, pheat_path, &w);
+
+ /* Handle fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ res = res_simul;
goto error;
}
+
+ if(pheat_path) {
+ pheat_path->status = res_simul == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res != RES_OK) goto error;
+ }
+ }
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ if(res_simul == RES_OK) {
+ /* Update global accumulators */
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp.sum += w; acc_temp.sum2 += w*w; ++acc_temp.count;
+ acc_time.sum += usec; acc_time.sum2 += usec*usec; ++acc_time.count;
+ }
}
- accum->sum_weights = sum_weights;
- accum->sum_weights_sqr = sum_weights_sqr;
- accum->nweights = N;
- accum->nfailures = nrealisations - N;
+ /* Setup the pixel estimator */
+ ASSERT(acc_temp.count == acc_time.count);
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
exit:
return res;
@@ -128,22 +166,23 @@ solve_tile
struct ssp_rng* rng,
struct sdis_medium* mdm,
const struct sdis_camera* cam,
- const double time,
+ const double time_range[2],
const double fp_to_meter,
const double Tarad,
const double Tref,
const size_t origin[2], /* Tile origin in image plane */
const size_t size[2], /* #pixels in X and Y */
const size_t spp, /* #samples per pixel */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accums)
+ struct sdis_estimator_buffer* buf)
{
size_t mcode; /* Morton code of the tile pixel */
size_t npixels;
res_T res = RES_OK;
- ASSERT(scn && rng && mdm && cam && spp && origin && accums && Tref >= 0);
- ASSERT(size &&size[0] && size[1]);
- ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(scn && rng && mdm && cam && spp && origin && Tref >= 0);
+ ASSERT(size &&size[0] && size[1] && buf);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range);
/* Adjust the #pixels to process them wrt a morton order */
npixels = round_up_pow2(MMAX(size[0], size[1]));
@@ -151,19 +190,21 @@ solve_tile
FOR_EACH(mcode, 0, npixels) {
size_t ipix[2];
- struct sdis_accum* accum;
+ struct sdis_estimator* estimator;
ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
if(ipix[0] >= size[0]) continue;
ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
if(ipix[1] >= size[1]) continue;
- accum = accums + ipix[1]*size[0] + ipix[0];
ipix[0] = ipix[0] + origin[0];
ipix[1] = ipix[1] + origin[1];
- res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, ipix,
- spp, pix_sz, accum);
+ /* Fetch the pixel estimator */
+ estimator = estimator_buffer_grab(buf, ipix[0], ipix[1]);
+
+ res = solve_pixel(scn, rng, mdm, cam, time_range, fp_to_meter, Tarad, Tref,
+ ipix, spp, register_paths, pix_sz, estimator);
if(res != RES_OK) goto error;
}
@@ -360,39 +401,48 @@ res_T
sdis_solve_camera
(struct sdis_scene* scn,
const struct sdis_camera* cam,
- const double time,
+ const double time_range[2],
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
const size_t width, /* #pixels in X */
const size_t height, /* #pixels in Y */
const size_t spp, /* #samples per pixel */
- sdis_write_accums_T writer,
- void* writer_data)
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_estimator_buffer** out_buf)
{
#define TILE_SIZE 32 /* definition in X & Y of a tile */
STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
+ struct sdis_estimator_buffer* buf = NULL;
struct sdis_medium* medium = NULL;
- struct darray_accum* tiles = NULL;
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
size_t ntiles_x, ntiles_y, ntiles;
double pix_sz[2]; /* Size of a pixel in the normalized image plane */
int64_t mcode; /* Morton code of a tile */
+ size_t nrealisations;
+ size_t nsuccesses;
+ size_t ix, iy;
size_t i;
ATOMIC res = RES_OK;
- if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
- || !height || !spp || !writer) {
+ if(!scn || !cam || fp_to_meter <= 0 || Tref < 0 || !width || !height || !spp
+ || !out_buf) {
res = RES_BAD_ARG;
goto error;
}
-
if(scene_is_2d(scn)) {
log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
goto error;
}
+ if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0]
+ || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* Retrieve the medium in which the submitted position lies */
res = scene_get_medium(scn, cam->position, NULL, &medium);
@@ -422,15 +472,6 @@ sdis_solve_camera
if(res != RES_OK) goto error;
}
- /* Allocate per thread buffer of accumulations */
- tiles = darray_tile_data_get(&scn->dev->tiles);
- ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
- FOR_EACH(i, 0, scn->dev->nthreads) {
- const size_t naccums = TILE_SIZE * TILE_SIZE;
- res = darray_accum_resize(tiles+i, naccums);
- if(res != RES_OK) goto error;
- }
-
ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
@@ -439,13 +480,16 @@ sdis_solve_camera
pix_sz[0] = 1.0 / (double)width;
pix_sz[1] = 1.0 / (double)height;
+ /* Create the global estimator */
+ res = estimator_buffer_create(scn->dev, width, height, &buf);
+ if(res != RES_OK) goto error;
+
omp_set_num_threads((int)scn->dev->nthreads);
#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};
const int ithread = omp_get_thread_num();
- struct sdis_accum* accums = NULL;
struct ssp_rng* rng = rngs[ithread];
res_T res_local = RES_OK;
@@ -462,25 +506,40 @@ sdis_solve_camera
tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]);
- /* Fetch the accumulations buffer */
- accums = darray_accum_data_get(tiles+ithread);
-
/* Draw the tile */
- res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
- Tref, tile_org, tile_sz, spp, pix_sz, accums);
+ res_local = solve_tile(scn, rng, medium, cam, time_range, fp_to_meter,
+ Tarad, Tref, tile_org, tile_sz, spp, register_paths, pix_sz, buf);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
}
+ }
- /* Write the accumulations */
- res_local = writer(writer_data, tile_org, tile_sz, accums);
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- continue;
+ /* Setup the accumulators of the whole estimator buffer */
+ acc_temp = ACCUM_NULL;
+ acc_time = ACCUM_NULL;
+ nsuccesses = 0;
+ FOR_EACH(iy, 0, height) {
+ FOR_EACH(ix, 0, width) {
+ const struct sdis_estimator* estimator;
+ SDIS(estimator_buffer_at(buf, ix, iy, &estimator));
+ acc_temp.sum += estimator->temperature.sum;
+ acc_temp.sum2 += estimator->temperature.sum2;
+ acc_temp.count += estimator->temperature.count;
+ acc_time.sum += estimator->realisation_time.sum;
+ acc_time.sum2 += estimator->realisation_time.sum2;
+ acc_time.count += estimator->realisation_time.count;
+ nsuccesses += estimator->nrealisations;
}
}
+ nrealisations = width*height*spp;
+ ASSERT(acc_temp.count == acc_time.count);
+ ASSERT(acc_temp.count == nsuccesses);
+ estimator_buffer_setup_realisations_count(buf, nrealisations, nsuccesses);
+ estimator_buffer_setup_temperature(buf, acc_temp.sum, acc_temp.sum2);
+ estimator_buffer_setup_realisation_time(buf, acc_time.sum, acc_time.sum2);
+
exit:
if(rngs) {
FOR_EACH(i, 0, scn->dev->nthreads) {
@@ -489,6 +548,7 @@ exit:
MEM_RM(scn->dev->allocator, rngs);
}
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_buf) *out_buf = buf;
return (res_T)res;
error:
goto exit;
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -21,6 +21,8 @@
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include <rsys/clock_time.h>
+#include <rsys/rsys.h>
#include <star/ssp.h>
#include <omp.h>
@@ -97,7 +99,8 @@ XD(solve_boundary)
struct sdis_green_function** greens = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
- struct accum* accums = NULL;
+ struct accum* acc_temps = NULL;
+ struct accum* acc_times = NULL;
size_t i;
size_t view_nprims;
int64_t irealisation;
@@ -185,9 +188,13 @@ XD(solve_boundary)
if(res != RES_OK) goto error;
}
- /* Create the per thread accumulator */
- accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
- if(!accums) { res = RES_MEM_ERR; goto error; }
+ /* Create the per thread accumulators */
+ acc_temps = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
+ if(!acc_temps) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
/* Create the per thread green function */
if(out_green) {
@@ -208,9 +215,11 @@ XD(solve_boundary)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
- struct accum* accum = &accums[ithread];
+ struct accum* acc_temp = &acc_temps[ithread];
+ struct accum* acc_time = &acc_times[ithread];
struct green_path_handle* pgreen_path = NULL;
struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
struct sdis_heat_path* pheat_path = NULL;
@@ -223,9 +232,13 @@ XD(solve_boundary)
float st[DIM-1];
double time;
res_T res_local = RES_OK;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ /* Begin time registration */
+ time_current(&t0);
+
if(!out_green) {
time = sample_time(rng, time_range);
if(register_paths) {
@@ -266,21 +279,18 @@ XD(solve_boundary)
side = sides[prim.prim_id];
/* Invoke the boundary realisation */
- res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
+ res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w);
- /* Update the MC accumulators */
- if(res_local == RES_OK) {
- accum->sum += w;
- accum->sum2 += w*w;
- ++accum->count;
- } else if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
+ /* Fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
continue;
}
+ /* Register heat path */
if(pheat_path) {
- pheat_path->status = res_local == RES_OK
+ pheat_path->status = res_simul == RES_OK
? SDIS_HEAT_PATH_SUCCEED
: SDIS_HEAT_PATH_FAILED;
@@ -292,18 +302,36 @@ XD(solve_boundary)
if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
}
}
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ /* Update accumulators */
+ if(res_simul == RES_OK) {
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+ }
}
/* Setup the estimated temperature */
if(out_estimator) {
- struct accum acc;
- sum_accums(accums, scn->dev->nthreads, &acc);
- estimator_setup_realisations_count(estimator, nrealisations, acc.count);
- estimator_setup_temperature(estimator, acc.sum, acc.sum2);
+ struct accum acc_temp;
+ struct accum acc_time;
+
+ sum_accums(acc_temps, scn->dev->nthreads, &acc_temp);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_temp.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
}
/* Setup the green function */
if(out_green) {
+ struct accum acc_time;
+
/* Redux the per thread green function into the green of the 1st thread */
green = greens[0]; /* Return the green of the 1st thread */
greens[0] = NULL; /* Make invalid the 1st green for 'on exit' clean up*/
@@ -311,7 +339,8 @@ XD(solve_boundary)
if(res != RES_OK) goto error;
/* Finalize the estimated green */
- res = green_function_finalize(green, rng_proxy);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ res = green_function_finalize(green, rng_proxy, &acc_time);
if(res != RES_OK) goto error;
}
@@ -328,7 +357,8 @@ exit:
}
MEM_RM(scn->dev->allocator, greens);
}
- if(accums) MEM_RM(scn->dev->allocator, accums);
+ if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
if(scene) SXD(scene_ref_put(scene));
if(shape) SXD(shape_ref_put(shape));
if(view) SXD(scene_view_ref_put(view));
@@ -369,6 +399,7 @@ XD(solve_boundary_flux)
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
struct accum* acc_tp = NULL; /* Per thread temperature accumulator */
+ struct accum* acc_ti = NULL; /* Per thread realisation time accumulator */
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
@@ -458,6 +489,7 @@ XD(solve_boundary_flux)
if(!Dst) { res = RES_MEM_ERR; goto error; } \
} (void)0
ALLOC_ACCUMS(acc_tp);
+ ALLOC_ACCUMS(acc_ti);
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
@@ -470,9 +502,11 @@ XD(solve_boundary_flux)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
struct accum* acc_temp = &acc_tp[ithread];
+ struct accum* acc_time = &acc_ti[ithread];
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
@@ -489,9 +523,13 @@ XD(solve_boundary_flux)
double time;
int flux_mask = 0;
res_T res_local = RES_OK;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ /* Begin time registration */
+ time_current(&t0);
+
time = sample_time(rng, time_range);
/* Sample a position onto the boundary */
@@ -545,9 +583,18 @@ XD(solve_boundary_flux)
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_local = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
+
+ res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
solid_side, fp_to_meter, Tarad, Tref, flux_mask, T_brf);
- if(res_local == RES_OK) {
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
+ continue;
+ } else if(res_simul == RES_OK) { /* Update accumulators */
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
const double Tboundary = T_brf[0];
const double Tradiative = T_brf[1];
const double Tfluid = T_brf[2];
@@ -558,6 +605,10 @@ XD(solve_boundary_flux)
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
++acc_temp->count;
+ /* Time */
+ acc_time->sum += usec;
+ acc_time->sum2 += usec*usec;
+ ++acc_time->count;
/* Overwall flux */
acc_flux->sum += w_total;
acc_flux->sum2 += w_total*w_total;
@@ -570,25 +621,25 @@ XD(solve_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
- } else if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
- continue;
}
}
if(res != RES_OK) goto error;
/* Redux the per thread accumulators */
sum_accums(acc_tp, scn->dev->nthreads, &acc_tp[0]);
+ sum_accums(acc_ti, scn->dev->nthreads, &acc_ti[0]);
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
+ ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
estimator_setup_temperature(estimator, acc_tp[0].sum, acc_tp[0].sum2);
+ estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
@@ -599,6 +650,7 @@ exit:
MEM_RM(scn->dev->allocator, rngs);
}
if(acc_tp) MEM_RM(scn->dev->allocator, acc_tp);
+ if(acc_ti) MEM_RM(scn->dev->allocator, acc_ti);
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -20,6 +20,7 @@
#include "sdis_scene_c.h"
#include <rsys/algorithm.h>
+#include <rsys/clock_time.h>
#include <rsys/dynamic_array.h>
#include "sdis_Xd_begin.h"
@@ -215,7 +216,8 @@ XD(solve_medium)
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
struct sdis_estimator* estimator = NULL;
- struct accum* accums = NULL;
+ struct accum* acc_temps = NULL;
+ struct accum* acc_times = NULL;
int64_t irealisation;
int cumul_is_init = 0;
size_t i;
@@ -257,9 +259,13 @@ XD(solve_medium)
if(res != RES_OK) goto error;
}
- /* Create the per thread accumulator */
- accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
- if(!accums) { res = RES_MEM_ERR; goto error; }
+ /* Create the per thread accumulators */
+ acc_temps = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
+ if(!acc_temps) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
/* Compute the enclosure cumulative */
darray_enclosure_cumul_init(scn->dev->allocator, &cumul);
@@ -286,9 +292,11 @@ XD(solve_medium)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
- struct accum* accum = accums + ithread;
+ struct accum* acc_temp = &acc_temps[ithread];
+ struct accum* acc_time = &acc_times[ithread];
struct green_path_handle* pgreen_path = NULL;
struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
const struct enclosure* enc = NULL;
@@ -298,9 +306,12 @@ XD(solve_medium)
double time;
double pos[DIM];
res_T res_local = RES_OK;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ time_current(&t0);
+
if(!out_green) {
/* Sample the time */
time = sample_time(rng, time_range);
@@ -331,19 +342,17 @@ XD(solve_medium)
}
/* Run a probe realisation */
- res_local = XD(probe_realisation)((size_t)irealisation, scn, rng, mdm, pos,
+ res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, mdm, pos,
time, fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &weight);
- if(res_local != RES_OK) {
- if(res_local != RES_BAD_OP) { ATOMIC_SET(&res, res_local); continue; }
- } else {
- accum->sum += weight;
- accum->sum2 += weight*weight;
- ++accum->count;
+
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
+ continue;
}
/* Finalize the registered path */
if(pheat_path) {
- pheat_path->status = res_local == RES_OK
+ pheat_path->status = res_simul == RES_OK
? SDIS_HEAT_PATH_SUCCEED
: SDIS_HEAT_PATH_FAILED;
@@ -355,18 +364,36 @@ XD(solve_medium)
if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
}
}
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ /* Update accumulators */
+ if(res_simul == RES_OK) {
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp->sum += weight; acc_temp->sum2 += weight*weight; ++acc_temp->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+ }
}
if(res != RES_OK) goto error;
/* Setup the estimated temperature */
if(out_estimator) {
- struct accum acc;
- sum_accums(accums, scn->dev->nthreads, &acc);
- estimator_setup_realisations_count(estimator, nrealisations, acc.count);
- estimator_setup_temperature(estimator, acc.sum, acc.sum2);
+ struct accum acc_temp;
+ struct accum acc_time;
+
+ sum_accums(acc_temps, scn->dev->nthreads, &acc_temp);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_temp.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
}
if(out_green) {
+ struct accum acc_time;
+
/* Redux the per thread green function into the green of the 1st thread */
green = greens[0]; /* Return the green of the 1st thread */
greens[0] = NULL; /* Make invalid the 1st green for 'on exit' clean up*/
@@ -374,7 +401,8 @@ XD(solve_medium)
if(res != RES_OK) goto error;
/* Finalize the estimated green */
- res = green_function_finalize(green, rng_proxy);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ res = green_function_finalize(green, rng_proxy, &acc_time);
if(res != RES_OK) goto error;
}
@@ -391,7 +419,8 @@ exit:
}
MEM_RM(scn->dev->allocator, greens);
}
- if(accums) MEM_RM(scn->dev->allocator, accums);
+ if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
if(cumul_is_init) darray_enclosure_cumul_release(&cumul);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -20,6 +20,7 @@
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include <rsys/clock_time.h>
#include <star/ssp.h>
#include <omp.h>
@@ -47,7 +48,8 @@ XD(solve_probe)
struct sdis_green_function** greens = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
- struct accum* accums = NULL;
+ struct accum* acc_temps = NULL;
+ struct accum* acc_times = NULL;
int64_t irealisation = 0;
size_t i;
ATOMIC res = RES_OK;
@@ -88,9 +90,13 @@ XD(solve_probe)
if(res != RES_OK) goto error;
}
- /* Create the per thread accumulator */
- accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
- if(!accums) { res = RES_MEM_ERR; goto error; }
+ /* Create the per thread accumulators */
+ acc_temps = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
+ if(!acc_temps) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
/* Retrieve the medium in which the submitted position lies */
res = scene_get_medium(scn, position, NULL, &medium);
@@ -116,19 +122,25 @@ XD(solve_probe)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
- struct accum* accum = &accums[ithread];
+ struct accum* acc_temp = &acc_temps[ithread];
+ struct accum* acc_time = &acc_times[ithread];
struct green_path_handle* pgreen_path = NULL;
struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
struct sdis_heat_path* pheat_path = NULL;
struct sdis_heat_path heat_path;
double w = NaN;
double time;
- res_T res_local;
+ res_T res_local = RES_OK;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ /* Begin time registration */
+ time_current(&t0);
+
if(!out_green) {
time = sample_time(rng, time_range);
if(register_paths) {
@@ -145,19 +157,17 @@ XD(solve_probe)
pgreen_path = &green_path;
}
- res_local = XD(probe_realisation)((size_t)irealisation, scn, rng, medium,
+ res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, medium,
position, time, fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w);
- if(res_local == RES_OK) {
- accum->sum += w;
- accum->sum2 += w*w;
- ++accum->count;
- } else if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
+
+ /* Handle fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
continue;
}
if(pheat_path) {
- pheat_path->status = res_local == RES_OK
+ pheat_path->status = res_simul == RES_OK
? SDIS_HEAT_PATH_SUCCEED
: SDIS_HEAT_PATH_FAILED;
@@ -169,18 +179,36 @@ XD(solve_probe)
if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
}
}
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ /* Update accumulators */
+ if(res_simul == RES_OK) {
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+ }
}
if(res != RES_OK) goto error;
- /* Setup the estimated temperature */
+ /* Setup the estimated temperature and per realisation time */
if(out_estimator) {
- struct accum acc;
- sum_accums(accums, scn->dev->nthreads, &acc);
- estimator_setup_realisations_count(estimator, nrealisations, acc.count);
- estimator_setup_temperature(estimator, acc.sum, acc.sum2);
+ struct accum acc_temp;
+ struct accum acc_time;
+
+ sum_accums(acc_temps, scn->dev->nthreads, &acc_temp);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_temp.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
}
if(out_green) {
+ struct accum acc_time;
+
/* Redux the per thread green function into the green of the 1st thread */
green = greens[0]; /* Return the green of the 1st thread */
greens[0] = NULL; /* Make invalid the 1st green for 'on exit' clean up*/
@@ -188,7 +216,8 @@ XD(solve_probe)
if(res != RES_OK) goto error;
/* Finalize the estimated green */
- res = green_function_finalize(green, rng_proxy);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ res = green_function_finalize(green, rng_proxy, &acc_time);
if(res != RES_OK) goto error;
}
@@ -205,7 +234,8 @@ exit:
}
MEM_RM(scn->dev->allocator, greens);
}
- if(accums) MEM_RM(scn->dev->allocator, accums);
+ if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_green) *out_green = green;
if(out_estimator) *out_estimator = estimator;
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -20,6 +20,7 @@
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include <rsys/clock_time.h>
#include <star/ssp.h>
#include <omp.h>
@@ -48,7 +49,8 @@ XD(solve_probe_boundary)
struct sdis_green_function** greens = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
- struct accum* accums = NULL;
+ struct accum* acc_temps = NULL;
+ struct accum* acc_times = NULL;
int64_t irealisation = 0;
size_t i;
ATOMIC res = RES_OK;
@@ -131,8 +133,12 @@ XD(solve_probe_boundary)
}
/* Create the per thread accumulator */
- accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
- if(!accums) { res = RES_MEM_ERR; goto error; }
+ acc_temps = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
+ if(!acc_temps) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
/* Create the per thread green function */
if(out_green) {
@@ -154,19 +160,25 @@ XD(solve_probe_boundary)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
- struct accum* accum = &accums[ithread];
+ struct accum* acc_temp = &acc_temps[ithread];
+ struct accum* acc_time = &acc_times[ithread];
struct green_path_handle* pgreen_path = NULL;
struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
struct sdis_heat_path* pheat_path = NULL;
struct sdis_heat_path heat_path;
double w = NaN;
double time;
- res_T res_local;
+ res_T res_local = RES_OK;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ /* Begin time registration */
+ time_current(&t0);
+
if(!out_green) {
time = sample_time(rng, time_range);
if(register_paths) {
@@ -182,19 +194,17 @@ XD(solve_probe_boundary)
pgreen_path = &green_path;
}
- res_local = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
+ res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side,
fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w);
- if(res_local == RES_OK) {
- accum->sum += w;
- accum->sum2 += w*w;
- ++accum->count;
- } else if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
+
+ /* Handle fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
continue;
}
if(pheat_path) {
- pheat_path->status = res_local == RES_OK
+ pheat_path->status = res_simul == RES_OK
? SDIS_HEAT_PATH_SUCCEED
: SDIS_HEAT_PATH_FAILED;
@@ -206,18 +216,36 @@ XD(solve_probe_boundary)
if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
}
}
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ /* Update accumulators */
+ if(res_simul == RES_OK) {
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+ }
}
if(res != RES_OK) goto error;
- /* Setup the estimated temperature */
+ /* Setup the estimated temperature and per realisation time */
if(out_estimator) {
- struct accum acc;
- sum_accums(accums, scn->dev->nthreads, &acc);
- estimator_setup_realisations_count(estimator, nrealisations, acc.count);
- estimator_setup_temperature(estimator, acc.sum, acc.sum2);
+ struct accum acc_temp;
+ struct accum acc_time;
+
+ sum_accums(acc_temps, scn->dev->nthreads, &acc_temp);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_temp.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
}
if(out_green) {
+ struct accum acc_time;
+
/* Redux the per thread green function into the green of the 1st thread */
green = greens[0]; /* Return the green of the 1st thread */
greens[0] = NULL; /* Make invalid the 1st green for 'on exit' clean up*/
@@ -225,7 +253,8 @@ XD(solve_probe_boundary)
if(res != RES_OK) goto error;
/* Finalize the estimated green */
- res = green_function_finalize(green, rng_proxy);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ res = green_function_finalize(green, rng_proxy, &acc_time);
if(res != RES_OK) goto error;
}
@@ -242,7 +271,8 @@ exit:
}
MEM_RM(scn->dev->allocator, greens);
}
- if(accums) MEM_RM(scn->dev->allocator, accums);
+ if(acc_temps) MEM_RM(scn->dev->allocator, acc_temps);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_green) *out_green = green;
if(out_estimator) *out_estimator = estimator;
@@ -279,6 +309,7 @@ XD(solve_probe_boundary_flux)
enum sdis_side solid_side, fluid_side;
struct sdis_interface_fragment frag;
struct accum* acc_tp = NULL; /* Per thread temperature accumulator */
+ struct accum* acc_ti = NULL; /* Per thread realisation time */
struct accum* acc_fl = NULL; /* Per thread flux accumulator */
struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */
struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */
@@ -369,6 +400,7 @@ XD(solve_probe_boundary_flux)
if(!Dst) { res = RES_MEM_ERR; goto error; } \
} (void)0
ALLOC_ACCUMS(acc_tp);
+ ALLOC_ACCUMS(acc_ti);
ALLOC_ACCUMS(acc_fc);
ALLOC_ACCUMS(acc_fl);
ALLOC_ACCUMS(acc_fr);
@@ -387,19 +419,24 @@ XD(solve_probe_boundary_flux)
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
struct accum* acc_temp = &acc_tp[ithread];
+ struct accum* acc_time = &acc_ti[ithread];
struct accum* acc_flux = &acc_fl[ithread];
struct accum* acc_fcon = &acc_fc[ithread];
struct accum* acc_frad = &acc_fr[ithread];
double time, epsilon, hc, hr;
int flux_mask = 0;
double T_brf[3] = { 0, 0, 0 };
- res_T res_local;
+ res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+ /* Begin time registration */
+ time_current(&t0);
+
time = sample_time(rng, time_range);
/* Compute hr and hc */
@@ -412,9 +449,17 @@ XD(solve_probe_boundary_flux)
flux_mask = 0;
if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE;
if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE;
- res_local = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
+ res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time,
solid_side, fp_to_meter, Tarad, Tref, flux_mask, T_brf);
- if(res_local == RES_OK) {
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_simul);
+ continue;
+ } else if(res_simul == RES_OK) { /* Update accumulators */
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
const double Tboundary = T_brf[0];
const double Tradiative = T_brf[1];
const double Tfluid = T_brf[2];
@@ -425,6 +470,10 @@ XD(solve_probe_boundary_flux)
acc_temp->sum += Tboundary;
acc_temp->sum2 += Tboundary*Tboundary;
++acc_temp->count;
+ /* Time */
+ acc_time->sum += usec;
+ acc_time->sum2 += usec*usec;
+ ++acc_time->count;
/* Overwall flux */
acc_flux->sum += w_total;
acc_flux->sum2 += w_total*w_total;
@@ -437,25 +486,25 @@ XD(solve_probe_boundary_flux)
acc_frad->sum += w_rad;
acc_frad->sum2 += w_rad*w_rad;
++acc_frad->count;
- } else if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
- continue;
}
}
if(res != RES_OK) goto error;
/* Redux the per thread accumulators */
sum_accums(acc_tp, scn->dev->nthreads, &acc_tp[0]);
+ sum_accums(acc_ti, scn->dev->nthreads, &acc_ti[0]);
sum_accums(acc_fc, scn->dev->nthreads, &acc_fc[0]);
sum_accums(acc_fr, scn->dev->nthreads, &acc_fr[0]);
sum_accums(acc_fl, scn->dev->nthreads, &acc_fl[0]);
ASSERT(acc_tp[0].count == acc_fl[0].count);
+ ASSERT(acc_tp[0].count == acc_ti[0].count);
ASSERT(acc_tp[0].count == acc_fr[0].count);
ASSERT(acc_tp[0].count == acc_fc[0].count);
/* Setup the estimated values */
estimator_setup_realisations_count(estimator, nrealisations, acc_tp[0].count);
estimator_setup_temperature(estimator, acc_tp[0].sum, acc_tp[0].sum2);
+ estimator_setup_realisation_time(estimator, acc_ti[0].sum, acc_ti[0].sum2);
estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc[0].sum, acc_fc[0].sum2);
estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2);
estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2);
@@ -466,6 +515,7 @@ exit:
MEM_RM(scn->dev->allocator, rngs);
}
if(acc_tp) MEM_RM(scn->dev->allocator, acc_tp);
+ if(acc_ti) MEM_RM(scn->dev->allocator, acc_ti);
if(acc_fc) MEM_RM(scn->dev->allocator, acc_fc);
if(acc_fr) MEM_RM(scn->dev->allocator, acc_fr);
if(acc_fl) MEM_RM(scn->dev->allocator, acc_fl);
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -384,6 +384,7 @@ main(int argc, char** argv)
OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng));
FOR_EACH(isimul, 0, nsimuls) {
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_estimator* estimator;
struct sdis_estimator* estimator2;
struct sdis_green_function* green;
@@ -402,16 +403,18 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
u = (pos[0] + 1) / thickness;
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.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(eq_eps(T.E, ref, 3*T.SE) == 1);
/* Check green function */
OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, Tref, &green));
@@ -426,6 +429,8 @@ main(int argc, char** argv)
OK(sdis_solve_probe(scn, 10, pos, time_range, 1, -1, Tref,
SDIS_HEAT_PATH_ALL, &estimator));
OK(sdis_estimator_ref_put(estimator));
+
+ printf("\n");
}
/* Release memory */
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -390,6 +390,7 @@ main(int argc, char** argv)
OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng));
FOR_EACH(isimul, 0, nsimuls) {
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_estimator* estimator;
struct sdis_estimator* estimator2;
struct sdis_green_function* green;
@@ -407,11 +408,13 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
u = (pos[0] + 1) / thickness;
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
CHK(nfails + nreals == N);
@@ -431,6 +434,8 @@ main(int argc, char** argv)
OK(sdis_solve_probe
(scn, 10, pos, time_range, 1, -1, Tref, SDIS_HEAT_PATH_ALL, &estimator));
OK(sdis_estimator_ref_put(estimator));
+
+ printf("\n");
}
/* Release memory */
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -167,6 +167,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc mc_time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
@@ -266,7 +267,7 @@ main(int argc, char** argv)
/* Test in 3D for various time values. */
nu = (6 * H) / (RHO*CP);
Tinf = (H*(T0 + T1 + T2 + T3 + T4 + T5)) / (6 * H);
- printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ printf(">>> Temperature of the box at (%g %g %g)\n\n", SPLIT3(pos));
FOR_EACH(i, 0, 5) {
double time = i ? (double)i / nu : INF;
double time_range[2];
@@ -279,10 +280,12 @@ main(int argc, char** argv)
/* Solve in 3D */
OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &mc_time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
CHK(nfails + nreals == N);
- printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE);
if(nfails)
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
CHK(eq_eps(T.E, ref, T.SE * 3));
@@ -297,12 +300,13 @@ main(int argc, char** argv)
}
OK(sdis_estimator_ref_put(estimator));
+ printf("\n");
}
/* Test in 2D for various time values. */
nu = (4 * H) / (RHO*CP);
Tinf = (H * (T0 + T1 + T2 + T3)) / (4 * H);
- printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ printf(">>> Temperature of the square at (%g %g)\n\n", SPLIT2(pos));
FOR_EACH(i, 0, 5) {
double time = i ? (double)i / nu : INF;
double time_range[2];
@@ -318,7 +322,9 @@ main(int argc, char** argv)
OK(sdis_estimator_get_failure_count(estimator, &nfails));
CHK(nfails + nreals == N);
OK(sdis_estimator_get_temperature(estimator, &T));
- printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ OK(sdis_estimator_get_realisation_time(estimator, &mc_time));
+ printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE);
if(nfails)
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
CHK(eq_eps(T.E, ref, T.SE * 3));
@@ -333,6 +339,7 @@ main(int argc, char** argv)
}
OK(sdis_estimator_ref_put(estimator));
+ printf("\n");
}
OK(sdis_scene_ref_put(box_scn));
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -177,6 +177,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc mc_time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
@@ -283,7 +284,7 @@ main(int argc, char** argv)
nu = (HC0 + HC1 + HC2 + HC3 + HC4 + HC5) / (RHO * CP);
Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3 + HC4 * T4 + HC5 * T5)
/ (HC0 + HC1 + HC2 + HC3 + HC4 + HC5);
- printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ printf(">>> Temperature of the box at (%g %g %g)\n\n", SPLIT3(pos));
FOR_EACH(i, 0, 5) {
double time = i ? (double)i / nu : INF;
double time_range[2];
@@ -298,7 +299,9 @@ main(int argc, char** argv)
OK(sdis_estimator_get_failure_count(estimator, &nfails));
CHK(nfails + nreals == N);
OK(sdis_estimator_get_temperature(estimator, &T));
- printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ OK(sdis_estimator_get_realisation_time(estimator, &mc_time));
+ printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE);
if(nfails)
printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
CHK(eq_eps(T.E, ref, T.SE * 3));
@@ -313,12 +316,13 @@ main(int argc, char** argv)
}
OK(sdis_estimator_ref_put(estimator));
+ printf("\n");
}
/* Test in 2D for various time values. */
nu = (HC0 + HC1 + HC2 + HC3) / (RHO * CP);
Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3) / (HC0 + HC1 + HC2 + HC3);
- printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ printf(">>> Temperature of the square at (%g %g)\n\n", SPLIT2(pos));
FOR_EACH(i, 0, 5) {
double time = i ? (double)i / nu : INF;
double time_range[2];
@@ -332,7 +336,9 @@ main(int argc, char** argv)
OK(sdis_estimator_get_failure_count(estimator, &nfails));
CHK(nfails + nreals == N);
OK(sdis_estimator_get_temperature(estimator, &T));
- printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ OK(sdis_estimator_get_realisation_time(estimator, &mc_time));
+ printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE);
if(nfails)
printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
CHK(eq_eps(T.E, ref, T.SE * 3));
@@ -347,6 +353,7 @@ main(int argc, char** argv)
}
OK(sdis_estimator_ref_put(estimator));
+ printf("\n");
}
OK(sdis_scene_ref_put(box_scn));
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -135,6 +135,7 @@ solve(struct sdis_scene* scn, const double pos[])
struct sdis_estimator* estimator2;
struct sdis_green_function* green;
struct sdis_mc T;
+ struct sdis_mc time;
size_t nreals;
size_t nfails;
double ref;
@@ -152,6 +153,7 @@ solve(struct sdis_scene* scn, const double pos[])
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_scene_get_dimension(scn, &dim));
@@ -166,6 +168,7 @@ solve(struct sdis_scene* scn, const double pos[])
break;
default: FATAL("Unreachable code.\n"); break;
}
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
printf("Elapsed time = %s\n\n", dump);
diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c
@@ -0,0 +1,391 @@
+/* Copyright (C) 2016-2019 |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 <star/s3dut.h>
+#include <rsys/math.h>
+
+#define Tfluid 350.0 /* Temperature of the fluid in Kelvin */
+#define LAMBDA 10.0 /* Thermal conductivity */
+#define Hcoef 1.0 /* Convection coefficient */
+#define Pw 10000.0 /* Volumetric power */
+#define Nreals 10000 /* #realisations */
+#define UNKNOWN -1
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+compute_aabb
+ (const double* positions,
+ const size_t nvertices,
+ double lower[3],
+ double upper[3])
+{
+ size_t i;
+ CHK(positions && nvertices && lower && upper);
+
+ lower[0] = lower[1] = lower[2] = DBL_MAX;
+ upper[0] = upper[1] = upper[2] =-DBL_MAX;
+
+ FOR_EACH(i, 0, nvertices) {
+ lower[0] = MMIN(lower[0], positions[i*3+0]);
+ lower[1] = MMIN(lower[1], positions[i*3+1]);
+ lower[2] = MMIN(lower[2], positions[i*3+2]);
+ upper[0] = MMAX(upper[0], positions[i*3+0]);
+ upper[1] = MMAX(upper[1], positions[i*3+1]);
+ upper[2] = MMAX(upper[2], positions[i*3+2]);
+ }
+}
+
+static double
+trilinear_temperature(const double pos[3])
+{
+ const double a = 333;
+ const double b = 432;
+ const double c = 579;
+ double x, y, z;
+ CHK(pos[0] >= -10.0 && pos[0] <= 10.0);
+ CHK(pos[1] >= -10.0 && pos[1] <= 10.0);
+ CHK(pos[2] >= -10.0 && pos[2] <= 10.0);
+
+ x = (pos[0] + 10.0) / 20.0;
+ y = (pos[1] + 10.0) / 20.0;
+ z = (pos[2] + 10.0) / 20.0;
+
+ return a*x + b*y + c*z;
+}
+
+static double
+volumetric_temperature(const double pos[3], const double upper[3])
+{
+ const double beta = - 1.0 / 3.0 * Pw / (2*LAMBDA);
+ const double temp =
+ beta * (pos[0]*pos[0] - upper[0]*upper[0])
+ + beta * (pos[1]*pos[1] - upper[1]*upper[1])
+ + beta * (pos[2]*pos[2] - upper[2]*upper[2]);
+ CHK(temp > 0);
+ return temp;
+}
+
+static void
+print_estimation_result
+ (const struct sdis_estimator* estimator, const double Tref)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ CHK(nfails <= Nreals * 0.0005);
+ printf("T = %g ~ %g +/- %g [%g, %g]; #failures = %lu / %lu\n",
+ Tref, T.E, T.SE, T.E - 3*T.SE, T.E + 3*T.SE,
+ (unsigned long)nfails, (unsigned long)Nreals);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ struct s3dut_mesh_data msh;
+ struct sdis_interface* interf;
+};
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ const struct context* ctx = context;
+ CHK(ids && ctx && itri < ctx->msh.nprimitives);
+ ids[0] = ctx->msh.indices[itri*3+0];
+ ids[2] = ctx->msh.indices[itri*3+1];
+ ids[1] = ctx->msh.indices[itri*3+2];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ const struct context* ctx = context;
+ CHK(pos && ctx && ivert < ctx->msh.nvertices);
+ pos[0] = ctx->msh.positions[ivert*3+0];
+ pos[1] = ctx->msh.positions[ivert*3+1];
+ pos[2] = ctx->msh.positions[ivert*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ CHK(bound && ctx && itri < ctx->msh.nprimitives);
+ *bound = ctx->interf;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+enum profile {
+ PROFILE_UNKNOWN,
+ PROFILE_TRILINEAR,
+ PROFILE_VOLUMETRIC_POWER
+};
+struct interf {
+ enum profile profile;
+ double upper[3]; /* Upper bound of the scene */
+ double h;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf;
+ double temperature;
+ CHK(data != NULL && frag != NULL);
+ interf = sdis_data_cget(data);
+ switch(interf->profile) {
+ case PROFILE_UNKNOWN:
+ temperature = UNKNOWN;
+ break;
+ case PROFILE_VOLUMETRIC_POWER:
+ temperature = volumetric_temperature(frag->P, interf->upper);
+ break;
+ case PROFILE_TRILINEAR:
+ temperature = trilinear_temperature(frag->P);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->h;;
+}
+
+/*******************************************************************************
+ * Fluid
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL);
+ (void)data;
+ return Tfluid;
+}
+
+/*******************************************************************************
+ * Solid API
+ ******************************************************************************/
+struct solid {
+ double delta;
+ double cp;
+ double lambda;
+ double rho;
+ double temperature;
+ double power;
+};
+
+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_volumetric_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->power;
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_mesh* msh = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct interf* interf_param = NULL;
+ struct solid* solid_param = NULL;
+ struct context ctx;
+ double probe[3];
+ double lower[3];
+ double upper[3];
+ double size[3];
+ const double time[2] = {DBL_MAX, DBL_MAX};
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* 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_volumetric_power;
+
+ /* Create the solid medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->delta = 0.01;
+ solid_param->lambda = 10;
+ solid_param->cp = 1.0;
+ solid_param->rho = 1.0;
+ solid_param->temperature = -1;
+ solid_param->power = SDIS_VOLUMIC_POWER_NONE;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the interface */
+ OK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = 1;
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, &interf));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid super shape */
+ f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5;
+ f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5;
+ OK(s3dut_create_super_shape(&allocator, &f0, &f1, 1, 128, 64, &msh));
+ OK(s3dut_mesh_get_data(msh, &ctx.msh));
+
+ compute_aabb(ctx.msh.positions, ctx.msh.nvertices, lower, upper);
+
+ /* Create the scene */
+ ctx.interf = interf;
+ OK(sdis_scene_create(dev, ctx.msh.nprimitives, get_indices, get_interface,
+ ctx.msh.nvertices, get_position, &ctx, &scn));
+ /*dump_mesh(stdout, ctx.msh.positions,
+ ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/
+
+ /* Compute the delta of the solid random walk */
+ size[0] = upper[0] - lower[0];
+ size[1] = upper[1] - lower[1];
+ size[2] = upper[2] - lower[2];
+ solid_param->delta = MMIN(MMIN(size[0], size[1]), size[2]) / 20.0;
+
+ interf_param->upper[0] = upper[0];
+ interf_param->upper[1] = upper[1];
+ interf_param->upper[2] = upper[2];
+ interf_param->h = 1;
+
+ probe[0] = 0;
+ probe[1] = 0;
+ probe[2] = 0;
+
+ /* Launch probe estimation with trilinear profile set at interfaces */
+ interf_param->profile = PROFILE_TRILINEAR;
+ OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, trilinear_temperature(probe));
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Launch probe estimation with volumetric power profile set at interfaces */
+ interf_param->profile = PROFILE_VOLUMETRIC_POWER;
+ solid_param->power = Pw;
+ OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, volumetric_temperature(probe, upper));
+ solid_param->power = SDIS_VOLUMIC_POWER_NONE;
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Launch medium integration */
+ interf_param->profile = PROFILE_UNKNOWN;
+ OK(sdis_solve_medium(scn, Nreals, solid, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, Tfluid);
+ /*dump_heat_paths(stdout, estimator);*/
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Release */
+ OK(s3dut_mesh_ref_put(msh));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_scene_ref_put(scn));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -148,15 +148,18 @@ check_estimator
const double ref)
{
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
size_t nreals;
size_t nfails;
CHK(estimator && nrealisations);
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
printf("%g ~ %g +/- %g\n", ref, T.E, T.SE);
- printf("#failures = %lu/%lu\n",
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu/%lu\n\n",
(unsigned long)nfails, (unsigned long)nrealisations);
CHK(nfails + nreals == nrealisations);
CHK(nfails < N/1000);
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -433,39 +433,59 @@ geometry_add_shape
}
static void
-dump_image(const struct sdis_accum_buffer* buf)
+dump_image(const struct sdis_estimator_buffer* buf)
{
- struct sdis_accum_buffer_layout layout = SDIS_ACCUM_BUFFER_LAYOUT_NULL;
struct image img;
double* temps = NULL;
- const struct sdis_accum* accums = NULL;
double Tmax = -DBL_MAX;
double Tmin = DBL_MAX;
double norm;
- size_t i, ix, iy;
+ size_t definition[2];
+ size_t ix, iy;
CHK(buf != NULL);
- OK(sdis_accum_buffer_get_layout(buf, &layout));
+ OK(sdis_estimator_buffer_get_definition(buf, definition));
- temps = mem_alloc(layout.width*layout.height*sizeof(double));
+ temps = mem_alloc(definition[0]*definition[1]*sizeof(double));
CHK(temps != NULL);
- OK(sdis_accum_buffer_map(buf, &accums));
-
/* 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);
+ FOR_EACH(iy, 0, definition[1]) {
+ FOR_EACH(ix, 0, definition[0]) {
+ const struct sdis_estimator* estimator;
+ struct sdis_mc T;
+ struct sdis_mc time;
+ size_t nreals;
+ size_t nfails;
+
+ BA(sdis_estimator_buffer_at(NULL, ix, iy, &estimator));
+ BA(sdis_estimator_buffer_at(buf, definition[0]+1, iy, &estimator));
+ BA(sdis_estimator_buffer_at(buf, ix, definition[1]+1, &estimator));
+ BA(sdis_estimator_buffer_at(buf, ix, iy, NULL));
+ OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ CHK(nreals + nfails == SPP);
+ CHK(T.E > 0);
+ CHK(time.E > 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;
- FOR_EACH(ix, 0, layout.width) {
- row[ix] = row_accums[ix].sum_weights / (double)row_accums[ix].nweights;
+ FOR_EACH(iy, 0, definition[1]) {
+ double* row = temps + iy * definition[0];
+ FOR_EACH(ix, 0, definition[0]) {
+ const struct sdis_estimator* estimator;
+ struct sdis_mc T;
+
+ OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ row[ix] = T.E;
Tmax = MMAX(row[ix], Tmax);
Tmin = MMIN(row[ix], Tmin);
}
@@ -481,10 +501,10 @@ dump_image(const struct sdis_accum_buffer* buf)
OK(image_init(NULL, &img));
OK(image_setup(&img, IMG_WIDTH, IMG_HEIGHT, IMG_WIDTH*3, IMAGE_RGB8, NULL));
- FOR_EACH(iy, 0, layout.height) {
- const double* src_row = temps + iy*layout.width;
+ FOR_EACH(iy, 0, definition[1]) {
+ const double* src_row = temps + iy*definition[0];
char* dst_row = img.pixels + iy*img.pitch;
- FOR_EACH(ix, 0, layout.width) {
+ FOR_EACH(ix, 0, definition[0]) {
unsigned char* pixels = (unsigned char*)
(dst_row + ix * sizeof_image_format(img.format));
const unsigned char T = (unsigned char)
@@ -509,9 +529,11 @@ main(int argc, char** argv)
struct geometry geom = GEOMETRY_NULL;
struct s3dut_mesh* msh = NULL;
struct s3dut_mesh_data msh_data;
- struct sdis_accum_buffer* buf = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_camera* cam = NULL;
struct sdis_device* dev = NULL;
+ struct sdis_estimator_buffer* buf = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid0 = NULL;
struct sdis_medium* fluid1 = NULL;
@@ -522,9 +544,12 @@ main(int argc, char** argv)
struct solid solid_param = SOLID_NULL;
struct interf interface_param = INTERF_NULL;
size_t ntris, npos;
+ size_t nreals, nfails;
+ size_t definition[2];
double pos[3];
double tgt[3];
double up[3];
+ double trange[2] = {INF, INF};
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
@@ -590,6 +615,11 @@ main(int argc, char** argv)
geometry_get_interface, npos, geometry_get_position,
&geom, &scn));
+#if 0
+ dump_mesh(stdout, geom.positions, sa_size(geom.positions)/3, geom.indices,
+ sa_size(geom.indices)/3);
+#endif
+
/* Setup the camera */
d3(pos, 3, 3, 3);
d3(tgt, 0, 0, 0);
@@ -599,28 +629,80 @@ main(int argc, char** argv)
OK(sdis_camera_set_fov(cam, MDEG2RAD(70)));
OK(sdis_camera_look_at(cam, pos, tgt, up));
- /* Create the accum buffer */
- OK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf));
-
#if 0
dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
exit(0);
#endif
+ BA(sdis_solve_camera(NULL, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, NULL, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, NULL, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 0, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, -1, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, 0, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, 0,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ 0, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, NULL));
+
+ trange[0] = -1;
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ trange[0] = 10; trange[1] = 1;
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ trange[0] = trange[1] = INF;
/* Launch the simulation */
- OK(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
- sdis_accum_buffer_write, buf));
+ OK(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+
+ BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals));
+ BA(sdis_estimator_buffer_get_realisation_count(buf, NULL));
+ OK(sdis_estimator_buffer_get_realisation_count(buf, &nreals));
+
+ BA(sdis_estimator_buffer_get_failure_count(NULL, &nfails));
+ BA(sdis_estimator_buffer_get_failure_count(buf, NULL));
+ OK(sdis_estimator_buffer_get_failure_count(buf, &nfails));
+
+ BA(sdis_estimator_buffer_get_temperature(NULL, &T));
+ BA(sdis_estimator_buffer_get_temperature(buf, NULL));
+ OK(sdis_estimator_buffer_get_temperature(buf, &T));
+
+ BA(sdis_estimator_buffer_get_realisation_time(NULL, &time));
+ BA(sdis_estimator_buffer_get_realisation_time(buf, NULL));
+ OK(sdis_estimator_buffer_get_realisation_time(buf, &time));
+
+ CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP);
+
+ fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE);
+ fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE);
+ fprintf(stderr, "#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));
+
+ BA(sdis_estimator_buffer_get_definition(NULL, definition));
+ BA(sdis_estimator_buffer_get_definition(buf, NULL));
+ OK(sdis_estimator_buffer_get_definition(buf, definition));
+ CHK(definition[0] == IMG_WIDTH);
+ CHK(definition[1] == IMG_HEIGHT);
/* Write the image */
dump_image(buf);
/* Release memory */
+ OK(sdis_estimator_buffer_ref_put(buf));
OK(sdis_scene_ref_put(scn));
OK(sdis_camera_ref_put(cam));
OK(sdis_interface_ref_put(interf0));
OK(sdis_interface_ref_put(interf1));
OK(sdis_device_ref_put(dev));
- OK(sdis_accum_buffer_ref_put(buf));
geometry_release(&geom);
check_memory_allocator(&allocator);
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -203,6 +203,7 @@ main(int argc, char** argv)
struct s3dut_mesh* msh0 = NULL;
struct s3dut_mesh* msh1 = NULL;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* solid0 = NULL;
struct sdis_medium* solid1 = NULL;
@@ -362,8 +363,10 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
printf("Shape0 temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, N);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu/%lu\n\n", (unsigned long)nfails, N);
CHK(eq_eps(T.E, Tf0, T.SE));
CHK(nreals + nfails == N);
OK(sdis_estimator_ref_put(estimator));
@@ -372,8 +375,10 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
printf("Shape1 temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, N);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu/%lu\n\n", (unsigned long)nfails, N);
CHK(eq_eps(T.E, Tf1, T.SE));
CHK(nreals + nfails == N);
OK(sdis_estimator_ref_put(estimator));
@@ -398,10 +403,12 @@ main(int argc, char** argv)
BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
ref = Tf0 * v0/v + Tf1 * v1/v;
printf("Shape0 + Shape1 temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, Np);
CHK(eq_eps(T.E, ref, T.SE*3));
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -188,6 +188,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* solid0 = NULL;
struct sdis_medium* solid1 = NULL;
@@ -336,10 +337,12 @@ main(int argc, char** argv)
/* Estimate the temperature of the square */
OK(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
printf("Square temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE);
- printf("#failures = %lu / %lu\n", (unsigned long)nfails, N);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu / %lu\n\n", (unsigned long)nfails, N);
CHK(eq_eps(T.E, Tf0, T.SE));
CHK(nreals + nfails == N);
OK(sdis_estimator_ref_put(estimator));
@@ -347,10 +350,12 @@ main(int argc, char** argv)
/* Estimate the temperature of the disk */
OK(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
printf("Disk temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE);
- printf("#failures = %lu / %lu\n", (unsigned long)nfails, N);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu / %lu\n\n", (unsigned long)nfails, N);
CHK(eq_eps(T.E, Tf1, T.SE));
CHK(nreals + nfails == N);
OK(sdis_estimator_ref_put(estimator));
@@ -369,10 +374,12 @@ main(int argc, char** argv)
BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 0, 0, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
ref = Tf0 * a0/a + Tf1 * a1/a;
printf("Square + Disk temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu / %lu\n", (unsigned long)nfails, Np);
CHK(eq_eps(T.E, ref, 3*T.SE));
CHK(nreals + nfails == Np);
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -249,6 +249,7 @@ main(int argc, char** argv)
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
struct sdis_mc F = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
@@ -379,9 +380,14 @@ main(int argc, char** argv)
BA(sdis_estimator_get_temperature(NULL, &T));
OK(sdis_estimator_get_temperature(estimator, &T));
+ BA(sdis_estimator_get_realisation_time(estimator, NULL));
+ BA(sdis_estimator_get_realisation_time(NULL, &time));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
ref = 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
CHK(nfails + nreals == N);
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -146,6 +146,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -246,11 +247,13 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
/* Print the estimation results */
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("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -143,6 +143,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -244,11 +245,14 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
/* Print the estimation results */
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("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -168,6 +168,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -302,17 +303,19 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
/* Print the estimation results */
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("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
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));
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Check green function */
OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green));
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -165,6 +165,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -295,11 +296,13 @@ main(int argc, char** argv)
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
/* Print the estimation results */
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("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -135,6 +135,7 @@ main(int argc, char** argv)
{
struct mem_allocator allocator;
struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
@@ -204,12 +205,13 @@ main(int argc, char** argv)
OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
OK(sdis_estimator_get_failure_count(estimator, &nfails));
-
OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
ref = 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
CHK(nfails + nreals == N);
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -262,11 +262,15 @@ check_green_function(struct sdis_green_function* green)
CHK(E + SE >= mc.E - mc.SE);
CHK(E - SE <= mc.E + mc.SE);
+ OK(sdis_estimator_get_realisation_time(estimator, &mc));
+ printf("Green per realisation time (in usec) = %g +/- %g\n",
+ mc.E, mc.SE);
+
OK(sdis_estimator_ref_put(estimator));
}
void
-dump_heat_paths(FILE* stream, struct sdis_estimator* estimator)
+dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator)
{
const struct sdis_heat_path* path;
size_t ipath;
@@ -356,6 +360,6 @@ dump_heat_paths(FILE* stream, struct sdis_estimator* estimator)
}
}
fprintf(stream, "LOOKUP_TABLE path_type 2\n");
- fprintf(stream, "0.0 0.0 1.0 1.0\n"); /* 0.0 = Bleu: success */
+ fprintf(stream, "0.0 0.0 1.0 1.0\n"); /* 0.0 = Blue: success */
fprintf(stream, "1.0 0.0 0.0 1.0\n"); /* 1.0 = Red: failure */
}
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -292,7 +292,7 @@ check_green_function
extern LOCAL_SYM void
dump_heat_paths
(FILE* stream,
- struct sdis_estimator* estimator);
+ const struct sdis_estimator* estimator);
#endif /* TEST_SDIS_UTILS_H */
diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c
@@ -0,0 +1,403 @@
+/* Copyright (C) 2016-2019 |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/clock_time.h>
+#include <rsys/math.h>
+
+#define Tf 100.0
+#define Power 10000.0
+#define H 50.0
+#define LAMBDA 100.0
+#define DELTA 0.4/*(1.0/2.0)*/
+#define N 100000
+#define LENGTH 10000.0
+
+/*
+ * 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
+ * surrounding fluid whose temperature is fixed to Tf.
+ *
+ *
+ * _\ Tf
+ * / /
+ * \__/
+ *
+ * ... -----Hboundary----- ...
+ *
+ * Lambda, Power
+ *
+ * ... -----Hboundary----- ...
+ *
+ * _\ Tf
+ * / /
+ * \__/
+ *
+ */
+
+static const double vertices_2d[4/*#vertices*/*2/*#coords per vertex*/] = {
+ LENGTH,-0.5,
+ -LENGTH,-0.5,
+ -LENGTH, 0.5,
+ LENGTH, 0.5
+};
+
+static const double vertices_3d[8/*#vertices*/*3/*#coords per vertex*/] = {
+ -LENGTH,-0.5,-LENGTH,
+ LENGTH,-0.5,-LENGTH,
+ -LENGTH, 0.5,-LENGTH,
+ LENGTH, 0.5,-LENGTH,
+ -LENGTH,-0.5, LENGTH,
+ LENGTH,-0.5, LENGTH,
+ -LENGTH, 0.5, LENGTH,
+ LENGTH, 0.5, LENGTH
+};
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_position_2d(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices_2d[ivert*2+0];
+ pos[1] = vertices_2d[ivert*2+1];
+}
+
+static void
+get_position_3d(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices_3d[ivert*3+0];
+ pos[1] = vertices_3d[ivert*3+1];
+ pos[2] = vertices_3d[ivert*3+2];
+}
+
+static void
+get_interface(const size_t iprim, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iprim];
+}
+
+/*******************************************************************************
+ * 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)
+{
+ char dump[128];
+ struct time t0, t1;
+ 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_2d = NULL;
+ struct sdis_scene* scn_3d = 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[12/*#max primitives*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals, nfails;
+ double pos[3];
+ double time_range[2] = { INF, INF };
+ double Tref;
+ double x;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
+ OK(sdis_data_ref_put(data));
+
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2));
+ OK(sdis_data_ref_put(data));
+
+ /* 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 */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ 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;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_adiabatic));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid fluid1 interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = H;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_solid_fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid fluid2 interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = H;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
+ &interf_solid_fluid2));
+ OK(sdis_data_ref_put(data));
+
+ /* Release the media */
+ OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_medium_ref_put(fluid2));
+ OK(sdis_medium_ref_put(solid));
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_solid_fluid2; /* Bottom */
+ interfaces[1] = interf_adiabatic; /* Left */
+ interfaces[2] = interf_solid_fluid1; /* Top */
+ interfaces[3] = interf_adiabatic; /* Right */
+
+ /* Create the 2D scene */
+ OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, get_interface,
+ square_nvertices, get_position_2d, interfaces, &scn_2d));
+
+ /* Map the interfaces to their box triangles */
+ interfaces[0] = interfaces[1] = interf_adiabatic; /* Front */
+ interfaces[2] = interfaces[3] = interf_adiabatic; /* Left */
+ interfaces[4] = interfaces[5] = interf_adiabatic; /* Back */
+ interfaces[6] = interfaces[7] = interf_adiabatic; /* Right */
+ interfaces[8] = interfaces[9] = interf_solid_fluid1; /* Top */
+ interfaces[10]= interfaces[11]= interf_solid_fluid2; /* Bottom */
+
+ /* Create the 3D scene */
+ OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, get_interface,
+ box_nvertices, get_position_3d, interfaces, &scn_3d));
+
+ /* Release the interfaces */
+ OK(sdis_interface_ref_put(interf_adiabatic));
+ OK(sdis_interface_ref_put(interf_solid_fluid1));
+ OK(sdis_interface_ref_put(interf_solid_fluid2));
+
+ pos[0] = 0;
+ pos[1] = 0;
+ pos[2] = 0;
+
+ x = pos[1];
+ Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA);
+
+ printf(">>> 2D\n");
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn_2d, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Elapsed time = %s\n", dump);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ 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);
+ OK(sdis_estimator_ref_put(estimator));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ printf("\n>>> 3D\n");
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn_3d, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Elapsed time = %s\n", dump);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ 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);
+ OK(sdis_estimator_ref_put(estimator));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ OK(sdis_scene_ref_put(scn_2d));
+ OK(sdis_scene_ref_put(scn_3d));
+ OK(sdis_device_ref_put(dev));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c
@@ -1,363 +0,0 @@
-/* Copyright (C) 2016-2019 |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/clock_time.h>
-#include <rsys/math.h>
-
-#define Tf 100.0
-#define Power 10000.0
-#define H 50.0
-#define LAMBDA 100.0
-#define DELTA (1.0/2.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
- * surrounding fluid whose temperature is fixed to Tf.
- *
- *
- * _\ Tf
- * / /
- * \__/
- *
- * ... -----Hboundary----- ...
- *
- * Lambda, Power
- *
- * ... -----Hboundary----- ...
- *
- * _\ Tf
- * / /
- * \__/
- *
- */
-
-static const double vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
- -1000000.5,-0.5,
- -1000000.5, 0.5,
- 1000000.5, 0.5,
- 1000000.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)
-{
- char dump[128];
- struct time t0, t1;
- 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 time_range[2] = { INF, INF };
- double Tref;
- double x;
- (void)argc, (void)argv;
-
- OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
- OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
-
- /* Create the fluid medium */
- fluid_shader.temperature = fluid_get_temperature;
- fluid_shader.calorific_capacity = dummy_medium_getter;
- fluid_shader.volumic_mass = dummy_medium_getter;
-
- OK(sdis_data_create
- (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = Tf;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
- OK(sdis_data_ref_put(data));
-
- OK(sdis_data_create
- (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = Tf;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2));
- OK(sdis_data_ref_put(data));
-
- /* 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 */
- OK(sdis_data_create
- (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- 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;
- OK(sdis_solid_create(dev, &solid_shader, data, &solid));
- OK(sdis_data_ref_put(data));
-
- /* Setup the interface shader */
- interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.front.temperature = interface_get_temperature;
-
- /* Create the adiabatic interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = 0;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
- &interf_adiabatic));
- OK(sdis_data_ref_put(data));
-
- /* Create the solid fluid1 interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = H;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
- &interf_solid_fluid1));
- OK(sdis_data_ref_put(data));
-
- /* Create the solid fluid2 interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = H;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
- &interf_solid_fluid2));
- OK(sdis_data_ref_put(data));
-
- /* Release the media */
- OK(sdis_medium_ref_put(fluid1));
- OK(sdis_medium_ref_put(fluid2));
- OK(sdis_medium_ref_put(solid));
-
- /* 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 */
- OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
- nvertices, get_position, interfaces, &scn));
-
- /* Release the interfaces */
- OK(sdis_interface_ref_put(interf_adiabatic));
- OK(sdis_interface_ref_put(interf_solid_fluid1));
- OK(sdis_interface_ref_put(interf_solid_fluid2));
-
- pos[0] = 0;
- pos[1] = 0.25;
-
- x = pos[1];
- Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA);
-
- time_current(&t0);
- OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
- printf("Elapsed time = %s\n", dump);
-
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- 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);
- OK(sdis_estimator_ref_put(estimator));
- CHK(nfails + nreals == N);
- CHK(nfails < N/1000);
- CHK(eq_eps(T.E, Tref, T.SE*3));
-
- OK(sdis_scene_ref_put(scn));
- OK(sdis_device_ref_put(dev));
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHK(mem_allocated_size() == 0);
- return 0;
-}
-
