commit 922ac5303c1794806926f5a83b9355e01a63d9ef
parent c994286e89e10ac7cbfbb9b78162a581328e4308
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 5 Apr 2024 10:49:11 +0200
Merge remote-tracking branch 'origin/feature_wos' into develop
Diffstat:
77 files changed, 3388 insertions(+), 997 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -13,3 +13,5 @@ rng_state
tags
*.pc
src/sdis_version.h
+*.vtk
+*.obj
diff --git a/Makefile b/Makefile
@@ -38,6 +38,7 @@ SRC =\
src/sdis_green.c\
src/sdis_heat_path.c\
src/sdis_heat_path_boundary.c\
+ src/sdis_heat_path_conductive.c\
src/sdis_interface.c\
src/sdis_log.c\
src/sdis_medium.c\
@@ -89,6 +90,8 @@ libsdis.o: $(OBJ)
echo "senc3d $(SENC3D_VERSION) not found" >&2; exit 1; fi
@if ! $(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp; then \
echo "star-sp $(SSP_VERSION) not found" >&2; exit 1; fi
+ @if ! $(PKG_CONFIG) --atleast-version $(SWF_VERSION) swf; then \
+ echo "swf $(SWF_VERSION) not found" >&2; exit 1; fi
@echo "config done" > $@
.SUFFIXES: .c .d .o
@@ -195,6 +198,10 @@ TEST_SRC =\
src/test_sdis_source.c\
src/test_sdis_transcient.c\
src/test_sdis_unstationary_atm.c\
+ src/test_sdis_unsteady.c\
+ src/test_sdis_unsteady_1d.c\
+ src/test_sdis_unsteady_analytic_profile.c\
+ src/test_sdis_unsteady_analytic_profile_2d.c\
src/test_sdis_volumic_power.c\
src/test_sdis_volumic_power4.c
TEST_SRC_LONG =\
@@ -268,6 +275,8 @@ test_all: test
clean_test:
@$(SHELL) make.sh clean_test $(TEST_SRC) $(TEST_SRC_MPI) $(TEST_SRC_LONG)
+ rm -f super_shape_2d.obj paths_wos_2d.vtk paths_delta_sphere_2d.vtk
+ rm -f super_shape_3d.obj paths_wos_3d.vtk paths_delta_sphere_3d.vtk
rm -f rng_state
################################################################################
@@ -295,6 +304,9 @@ src/test_sdis_solve_probe3_2d \
src/test_sdis_source.d \
src/test_sdis_transcient.d \
src/test_sdis_unstationary_atm.d \
+src/test_sdis_unsteady.d \
+src/test_sdis_unsteady_1d.d \
+src/test_sdis_unsteady_analytic_profile_2d.d \
src/test_sdis_utils.d \
src/test_sdis_volumic_power.d \
src/test_sdis_volumic_power2.d \
@@ -326,6 +338,9 @@ src/test_sdis_solve_probe3_2d.o \
src/test_sdis_source.o \
src/test_sdis_transcient.o \
src/test_sdis_unstationary_atm.o \
+src/test_sdis_unsteady.o \
+src/test_sdis_unsteady_1d.o \
+src/test_sdis_unsteady_analytic_profile_2d.o \
src/test_sdis_utils.o \
src/test_sdis_volumic_power.o \
src/test_sdis_volumic_power2.o \
@@ -357,6 +372,9 @@ test_sdis_solve_probe3_2d \
test_sdis_source \
test_sdis_transcient \
test_sdis_unstationary_atm \
+test_sdis_unsteady \
+test_sdis_unsteady_1d \
+test_sdis_unsteady_analytic_profile_2d \
test_sdis_volumic_power \
test_sdis_volumic_power2 \
test_sdis_volumic_power2_2d \
@@ -371,18 +389,21 @@ test_sdis_volumic_power4 \
src/test_sdis_draw_external_flux.d \
src/test_sdis_solid_random_walk_robustness.d \
src/test_sdis_solve_probe3.d \
+src/test_sdis_unsteady_analytic_profile.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis_draw_external_flux.o \
src/test_sdis_solid_random_walk_robustness.o \
src/test_sdis_solve_probe3.o \
+src/test_sdis_unsteady_analytic_profile.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
test_sdis_draw_external_flux \
test_sdis_solid_random_walk_robustness \
test_sdis_solve_probe3 \
+test_sdis_unsteady_analytic_profile \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS) $(S3DUT_LIBS)
diff --git a/config.mk b/config.mk
@@ -63,6 +63,10 @@ SSP_VERSION = 0.14
SSP_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags star-sp)
SSP_LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs star-sp)
+SWF_VERSION = 0.0
+SWF_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags swf)
+SWF_LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs swf)
+
# For tests only
S3DUT_VERSION = 0.4
S3DUT_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags s3dut)
@@ -75,6 +79,7 @@ DPDC_CFLAGS =\
$(SENC2D_CFLAGS)\
$(SENC3D_CFLAGS)\
$(SSP_CFLAGS)\
+ $(SWF_CFLAGS)\
$($(DISTRIB_PARALLELISM)_CFLAGS)\
-fopenmp
DPDC_LIBS =\
@@ -84,6 +89,7 @@ DPDC_LIBS =\
$(SENC2D_LIBS)\
$(SENC3D_LIBS)\
$(SSP_LIBS)\
+ $(SWF_LIBS)\
$($(DISTRIB_PARALLELISM)_LIBS)\
-lm\
-fopenmp
diff --git a/src/sdis.h b/src/sdis.h
@@ -50,6 +50,11 @@
#define SDIS_FLUX_NONE DBL_MAX /* <=> No flux */
#define SDIS_PRIMITIVE_NONE SIZE_MAX /* Invalid primitive */
+/* Syntactic sugar used to define whether a temperature is known or not */
+#define SDIS_TEMPERATURE_NONE NaN /* Unknown temperature */
+#define SDIS_TEMPERATURE_IS_KNOWN(Temp) (!IS_NaN(Temp))
+#define SDIS_TEMPERATURE_IS_UNKNOWN(Temp) (IS_NaN(Temp))
+
/* Forward declaration of external opaque data types */
struct logger;
struct mem_allocator;
@@ -91,6 +96,13 @@ enum sdis_scene_dimension {
SDIS_SCENE_3D
};
+enum sdis_diffusion_algorithm {
+ SDIS_DIFFUSION_DELTA_SPHERE,
+ SDIS_DIFFUSION_WOS, /* Walk on Sphere */
+ SDIS_DIFFUSION_ALGORITHMS_COUNT__,
+ SDIS_DIFFUSION_NONE = SDIS_DIFFUSION_ALGORITHMS_COUNT__
+};
+
/* Random walk vertex, i.e. a spatiotemporal position at a given step of the
* random walk. */
struct sdis_rwalk_vertex {
@@ -235,13 +247,14 @@ struct sdis_solid_shader {
* submitted position and time */
sdis_medium_getter_T volumic_power; /* In W.m^-3 */
- /* Initial/limit condition. A temperature < 0 means that the temperature is
- * unknown for the submitted random walk vertex.
+ /* Initial/limit condition. A temperature set to SDIS_TEMPERATURE_NONE
+ * means that the temperature is unknown for the submitted random walk vertex.
* This getter is always called at time >= t0 (see below). */
sdis_medium_getter_T temperature;
- /* The time until the initial condition is maintained for this solid;
- * can neither be negative nor infinity, default is 0. */
+ /* The time until the initial condition is maintained for this solid.
+ * Can be negative or set to +/- infinity to simulate a system that is always
+ * in the initial state or never reaches it, respectively. */
double t0;
};
#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, 0}
@@ -254,12 +267,14 @@ struct sdis_fluid_shader {
sdis_medium_getter_T calorific_capacity; /* In J.K^-1.kg^-1 */
sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */
- /* Initial/limit condition. A temperature < 0 means that the temperature is
- * unknown for the submitted random walk vertex.
+ /* Initial/limit condition. A temperature set to SDIS_TEMPERATURE_NONE
+ * means that the temperature is unknown for the submitted random walk vertex.
* This getter is always called at time >= t0 (see below). */
sdis_medium_getter_T temperature;
- /* The time until the initial condition is maintained for this fluid;
- * can neither be negative nor infinity, default is 0. */
+
+ /* The time until the initial condition is maintained for this fluid.
+ * Can be negative or set to +/- infinity to simulate a system that is always
+ * in the initial state or never reaches it, respectively. */
double t0;
};
#define SDIS_FLUID_SHADER_NULL__ {NULL, NULL, NULL, 0}
@@ -270,8 +285,8 @@ static const struct sdis_fluid_shader SDIS_FLUID_SHADER_NULL =
struct sdis_interface_side_shader {
/* Fixed temperature/flux. May be NULL if the temperature/flux is unknown
* onto the whole interface */
- sdis_interface_getter_T temperature; /* In Kelvin. < 0 <=> Unknown temp */
- sdis_interface_getter_T flux; /* In W.m^-2. SDIS_FLUX_NONE <=> no flux */
+ sdis_interface_getter_T temperature; /* [K]. SDIS_TEMPERATURE_NONE = Unknown */
+ sdis_interface_getter_T flux; /* [W.m^-2]. SDIS_FLUX_NONE = no flux */
/* Control the emissivity of the interface. May be NULL for solid/solid
* interface or if the emissivity is 0 onto the whole interface. */
@@ -439,7 +454,10 @@ struct sdis_ambient_radiative_temperature {
double temperature; /* In Kelvin */
double reference; /* Used to linearise the radiative transfer */
};
-#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ {-1, -1}
+#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ { \
+ SDIS_TEMPERATURE_NONE, \
+ SDIS_TEMPERATURE_NONE \
+}
static const struct sdis_ambient_radiative_temperature
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__;
@@ -475,7 +493,7 @@ struct sdis_scene_create_args {
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
- {0.0, -1.0}, /* Temperature range */ \
+ {SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE}, /* Temperature range */ \
NULL /* source */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
@@ -498,6 +516,8 @@ struct sdis_solve_probe_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
+
/* Signature of the estimated green function. The signature is ignored in an
* ordinary probe estimation. The signature of the green function can be
* queried to verify that it is the expected one with respect to the caller's
@@ -512,6 +532,7 @@ struct sdis_solve_probe_args {
SDIS_HEAT_PATH_NONE, /* Register paths mask */ \
NULL, /* RNG state */ \
SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE, /* Diffusion algorithm */ \
{0} /* Signature */ \
}
static const struct sdis_solve_probe_args SDIS_SOLVE_PROBE_ARGS_DEFAULT =
@@ -552,6 +573,8 @@ struct sdis_solve_probe_boundary_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
+
/* Signature of the estimated green function. The signature is ignored in an
* ordinary probe estimation. The signature of the green function can be
* queried to verify that it is the expected one with respect to the caller's
@@ -568,6 +591,7 @@ struct sdis_solve_probe_boundary_args {
SDIS_HEAT_PATH_NONE, \
NULL, /* RNG state */ \
SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE, /* Diffusion algorithm */ \
{0} /* Signature */ \
}
static const struct sdis_solve_probe_boundary_args
@@ -611,6 +635,8 @@ struct sdis_solve_boundary_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
+
/* Signature of the estimated green function. The signature is ignored in an
* ordinary probe estimation. The signature of the green function can be
* queried to verify that it is the expected one with respect to the caller's
@@ -627,6 +653,7 @@ struct sdis_solve_boundary_args {
SDIS_HEAT_PATH_NONE, \
NULL, /* RNG state */ \
SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE, /* Diffusion algorithm */ \
{0} /* Signature */ \
}
static const struct sdis_solve_boundary_args SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT =
@@ -646,6 +673,8 @@ struct sdis_solve_medium_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
+
/* Signature of the estimated green function. The signature is ignored in an
* ordinary probe estimation. The signature of the green function can be
* queried to verify that it is the expected one with respect to the caller's
@@ -660,6 +689,7 @@ struct sdis_solve_medium_args {
SDIS_HEAT_PATH_NONE, \
NULL, /* RNG state */ \
SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE, /* Diffusion algorithm */ \
{0} /* Signature */ \
}
static const struct sdis_solve_medium_args SDIS_SOLVE_MEDIUM_ARGS_DEFAULT =
@@ -678,6 +708,8 @@ struct sdis_solve_probe_boundary_flux_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \
10000, /* #realisations */ \
@@ -686,7 +718,8 @@ struct sdis_solve_probe_boundary_flux_args {
{DBL_MAX,DBL_MAX}, /* Time range */ \
1, /* Picard order */ \
NULL, /* RNG state */ \
- SSP_RNG_THREEFRY /* RNG type */ \
+ SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE /* Diffusion algorithm */ \
}
static const struct sdis_solve_probe_boundary_flux_args
SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT =
@@ -705,6 +738,8 @@ struct sdis_solve_boundary_flux_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \
10000, /* #realisations */ \
@@ -713,7 +748,8 @@ struct sdis_solve_boundary_flux_args {
{DBL_MAX,DBL_MAX}, /* Time range */ \
1, /* Picard order */ \
NULL, /* RNG state */ \
- SSP_RNG_THREEFRY /* RNG type */ \
+ SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE /* Diffusion algorithm */ \
}
static const struct sdis_solve_boundary_flux_args
SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT =
@@ -734,6 +770,8 @@ struct sdis_solve_camera_args {
struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
enum ssp_rng_type rng_type; /* RNG type to use */
+
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define SDIS_SOLVE_CAMERA_ARGS_DEFAULT__ { \
NULL, /* Camera */ \
@@ -743,7 +781,8 @@ struct sdis_solve_camera_args {
256, /* #realisations per pixel */ \
SDIS_HEAT_PATH_NONE, \
NULL, /* RNG state */ \
- SSP_RNG_THREEFRY /* RNG type */ \
+ SSP_RNG_THREEFRY, /* RNG type */ \
+ SDIS_DIFFUSION_DELTA_SPHERE /* Diffusion algorithm */ \
}
static const struct sdis_solve_camera_args SDIS_SOLVE_CAMERA_ARGS_DEFAULT =
SDIS_SOLVE_CAMERA_ARGS_DEFAULT__;
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -13,63 +13,9 @@
* 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_XD_BEGIN_H
-#define SDIS_XD_BEGIN_H
-
+#include "sdis.h"
#include <rsys/rsys.h>
-/* Forward declaration */
-struct green_path_handle;
-struct sdis_heat_path;
-
-struct rwalk_context {
- struct green_path_handle* green_path;
- struct sdis_heat_path* heat_path;
-
- double Tmin; /* Lower bound temperature */
- double Tmin2; /* Tmin^2 */
- double Tmin3; /* Tmin^3 */
-
- double That; /* Upper bound temperature */
- double That2; /* That^2 */
- double That3; /* That^3 */
-
- /* Maximum branchings i.e. the maximum number of times XD(sample_coupled_path)
- * can be called. It controls the number of ramifications of the heat path and
- * currently is correlated to the Picard order used to estimate the radiative
- * temperature. max_branchings == picard_order-1 */
- size_t max_branchings;
-
- /* Number of heat path branchings */
- size_t nbranchings;
-
- /* Id of the realisation (for debug) */
- size_t irealisation;
-};
-#define RWALK_CONTEXT_NULL__ { \
- NULL, /* Green path */ \
- NULL, /* Heat path */ \
- 0, /* Tmin */ \
- 0, /* Tmin^2 */ \
- 0, /* Tmin^3 */ \
- 0, /* That */ \
- 0, /* That^2 */ \
- 0, /* That^3 */ \
- 0, /* Max #branchings */ \
- SIZE_MAX, /* #branchings */ \
- SIZE_MAX /* realisation id */ \
-}
-static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
-
-static INLINE size_t
-get_picard_order(const struct rwalk_context* ctx)
-{
- ASSERT(ctx);
- return ctx->max_branchings + 1;
-}
-
-#endif /* SDIS_XD_BEGIN_H */
-
#ifdef SDIS_XD_BEGIN_H__
#error "This header is already included without its associated sdis_Xd_end.h file."
#endif
@@ -140,6 +86,8 @@ get_picard_order(const struct rwalk_context* ctx)
#define SDIS_3D_H
#endif
+struct rwalk_context;
+
/* Current state of the random walk */
struct XD(rwalk) {
struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -25,6 +25,7 @@
#include <star/s2d.h>
#include <star/s3d.h>
#include <star/ssp.h>
+#include <star/swf.h>
#include <omp.h>
@@ -251,6 +252,37 @@ error:
}
static INLINE res_T
+setup_starwf(struct sdis_device* dev)
+{
+ struct swf_H_tabulate_args H2d_args = SWF_H2D_TABULATE_ARGS_DEFAULT;
+ struct swf_H_tabulate_args H3d_args = SWF_H3D_TABULATE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(dev);
+
+ H2d_args.allocator = dev->allocator;
+ H3d_args.allocator = dev->allocator;
+
+ res = swf_H2d_tabulate(&H2d_args, &dev->H_2d);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to tabulate H2d function -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = swf_H3d_tabulate(&H3d_args, &dev->H_3d);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to tabulate H3d function -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
setup_mpi(struct sdis_device* dev, const struct sdis_device_create_args* args)
{
ASSERT(dev && args);
@@ -279,6 +311,8 @@ device_release(ref_T* ref)
dev = CONTAINER_OF(ref, struct sdis_device, ref);
if(dev->s2d_dev) S2D(device_ref_put(dev->s2d_dev));
if(dev->s3d_dev) S3D(device_ref_put(dev->s3d_dev));
+ if(dev->H_2d) SWF(tabulation_ref_put(dev->H_2d));
+ if(dev->H_3d) SWF(tabulation_ref_put(dev->H_3d));
if(dev->logger == &dev->logger__) logger_release(&dev->logger__);
ASSERT(flist_name_is_empty(&dev->interfaces_names));
ASSERT(flist_name_is_empty(&dev->media_names));
@@ -342,6 +376,8 @@ sdis_device_create
if(res != RES_OK) goto error;
res = setup_star3d(dev);
if(res != RES_OK) goto error;
+ res = setup_starwf(dev);
+ if(res != RES_OK) goto error;
res = setup_mpi(dev, args);
if(res != RES_OK) goto error;
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -36,6 +36,7 @@
struct mutex;
struct ssp_rng;
struct ssp_rng_proxy;
+struct swf_tabulation;
struct name { FITEM; void* mem; };
#define FITEM_TYPE name
@@ -64,6 +65,9 @@ struct sdis_device {
struct s2d_device* s2d_dev;
struct s3d_device* s3d_dev;
+ struct swf_tabulation* H_2d;
+ struct swf_tabulation* H_3d;
+
ref_T ref;
};
diff --git a/src/sdis_heat_path.c b/src/sdis_heat_path.c
@@ -27,12 +27,6 @@
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_convective_Xd.h"
-/* Generate the conductive path routines */
-#define SDIS_XD_DIMENSION 2
-#include "sdis_heat_path_conductive_Xd.h"
-#define SDIS_XD_DIMENSION 3
-#include "sdis_heat_path_conductive_Xd.h"
-
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h
@@ -23,22 +23,77 @@
#include <rsys/rsys.h>
/* Forward declarations */
+struct green_path_handle;
struct rwalk_2d;
struct rwalk_3d;
-struct rwalk_context;
struct sdis_scene;
struct ssp_rng;
struct temperature_2d;
struct temperature_3d;
+/*******************************************************************************
+ * Context of a random walk, i.e. its data concerning the current system and the
+ * solve parameters.
+ ******************************************************************************/
+ struct rwalk_context {
+ struct green_path_handle* green_path;
+ struct sdis_heat_path* heat_path;
+
+ double Tmin; /* Lower bound temperature */
+ double Tmin2; /* Tmin^2 */
+ double Tmin3; /* Tmin^3 */
+
+ double That; /* Upper bound temperature */
+ double That2; /* That^2 */
+ double That3; /* That^3 */
+
+ /* Maximum branchings i.e. the maximum number of times XD(sample_coupled_path)
+ * can be called. It controls the number of ramifications of the heat path and
+ * currently is correlated to the Picard order used to estimate the radiative
+ * temperature. max_branchings == picard_order-1 */
+ size_t max_branchings;
+
+ /* Number of heat path branchings */
+ size_t nbranchings;
+
+ /* Id of the realisation (for debug) */
+ size_t irealisation;
+
+ /* Algorithm used for the diffusive random walks,
+ * i.e. for sampling conductive paths */
+ enum sdis_diffusion_algorithm diff_algo;
+};
+#define RWALK_CONTEXT_NULL__ { \
+ NULL, /* Green path */ \
+ NULL, /* Heat path */ \
+ 0, /* Tmin */ \
+ 0, /* Tmin^2 */ \
+ 0, /* Tmin^3 */ \
+ 0, /* That */ \
+ 0, /* That^2 */ \
+ 0, /* That^3 */ \
+ 0, /* Max #branchings */ \
+ SIZE_MAX, /* #branchings */ \
+ SIZE_MAX, /* realisation id */ \
+ SDIS_DIFFUSION_NONE /* Diffusion algorithm */ \
+}
+static const struct rwalk_context RWALK_CONTEXT_NULL = RWALK_CONTEXT_NULL__;
+
+static INLINE size_t
+get_picard_order(const struct rwalk_context* ctx)
+{
+ ASSERT(ctx);
+ return ctx->max_branchings + 1;
+}
+
+/*******************************************************************************
+ * Heat path data structure used to record the geometry of sampled paths
+ ******************************************************************************/
/* Generate the dynamic array of heat vertices */
#define DARRAY_NAME heat_vertex
#define DARRAY_DATA struct sdis_heat_vertex
#include <rsys/dynamic_array.h>
-/*******************************************************************************
- * Heat path data structure
- ******************************************************************************/
struct sdis_heat_path {
/* List of the path vertices */
struct darray_heat_vertex vertices;
@@ -293,4 +348,3 @@ boundary_path_3d
struct temperature_3d* temperature);
#endif /* SDIS_HEAT_PATH_H */
-
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -54,7 +54,7 @@ XD(boundary_path)
/* Check if the boundary temperature is known */
tmp = interface_side_get_temperature(interf, &frag);
- if(tmp >= 0) {
+ if(SDIS_TEMPERATURE_IS_KNOWN(tmp)) {
T->value += tmp;
T->done = 1;
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -992,4 +992,40 @@ error:
goto exit;
}
+res_T
+XD(check_Tref)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ const double Tref,
+ const char* func_name)
+{
+ ASSERT(scn && pos && func_name);
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(Tref)) {
+ log_err(scn->dev,
+ "%s: invalid reference temperature at the position `"FORMAT_VECX"'.\n",
+ func_name, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmin)
+ || SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmax)) {
+ log_err(scn->dev,
+ "%s: invalid scene temperature range. "
+ "At least one boundary is unknown.\n",
+ func_name);
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ if(Tref > scn->tmax) {
+ log_err(scn->dev,
+ "%s: invalid maximum temperature `%gK'. The reference temperature `%gK' "
+ "at the position `"FORMAT_VECX"' is greater than this temperature.\n",
+ func_name, scn->tmax, Tref, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ return RES_OK;
+}
+
#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -28,39 +28,13 @@
* Helper functions
******************************************************************************/
static INLINE res_T
-XD(check_Tref)
- (const struct sdis_scene* scn,
- const double pos[DIM],
- const double Tref,
- const char* func_name)
-{
- ASSERT(scn && pos && func_name);
-
- if(Tref < 0) {
- log_err(scn->dev,
- "%s: invalid reference temperature `%gK' at the position `"FORMAT_VECX"'.\n",
- func_name, Tref, SPLITX(pos));
- return RES_BAD_OP_IRRECOVERABLE;
- }
- if(Tref > scn->tmax) {
- log_err(scn->dev,
- "%s: invalid maximum temperature `%gK'. The reference temperature `%gK' "
- "at the position `"FORMAT_VECX"' is greater than this temperature.\n",
- func_name, scn->tmax, Tref, SPLITX(pos));
- return RES_BAD_OP_IRRECOVERABLE;
- }
-
- return RES_OK;
-}
-
-static INLINE res_T
XD(rwalk_get_Tref)
(const struct sdis_scene* scn,
const struct XD(rwalk)* rwalk,
const struct XD(temperature)* T,
double* out_Tref)
{
- double Tref = -1;
+ double Tref = SDIS_TEMPERATURE_NONE;
res_T res = RES_OK;
ASSERT(rwalk && T && out_Tref);
@@ -213,7 +187,13 @@ XD(solid_fluid_boundary_picard1_path)
/* Compute the convective, conductive and the upper bound radiative coef */
h_conv = interface_get_convection_coef(interf, frag);
h_cond = lambda / delta_m;
- h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
+ if(epsilon <= 0) {
+ h_radi_hat = 0; /* No radiative transfert */
+ } else {
+ res = scene_check_temperature_range(scn);
+ if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; }
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
+ }
/* Compute a global upper bound coefficient */
h_hat = h_conv + h_cond + h_radi_hat;
@@ -305,6 +285,12 @@ XD(solid_fluid_boundary_picard1_path)
res = XD(rwalk_get_Tref)(scn, &rwalk_s, &T_s, &Tref_s);
if(res != RES_OK) goto error;
+ /* The reference temperatures must be known, as this is a radiative path.
+ * If this is not the case, an error should be reported before this point.
+ * Hence these assertions to detect unexpected behavior */
+ ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref));
+ ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref_s));
+
h_radi = BOLTZMANN_CONSTANT * epsilon *
( Tref*Tref*Tref
+ Tref*Tref * Tref_s
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -232,7 +232,15 @@ XD(solid_fluid_boundary_picardN_path)
/* Compute the convective, conductive and the upper bound radiative coef */
h_conv = interface_get_convection_coef(interf, frag);
h_cond = lambda / (delta * scn->fp_to_meter);
- h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon;
+ h_radi_hat = epsilon > 0 ? 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon : 0;
+
+ if(epsilon <= 0) {
+ h_radi_hat = 0; /* No radiative transfert */
+ } else {
+ res = scene_check_temperature_range(scn);
+ if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; }
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon;
+ }
/* Compute a global upper bound coefficient */
h_hat = h_conv + h_cond + h_radi_hat;
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -221,6 +221,23 @@ handle_external_net_flux_3d
struct temperature_3d* T);
/*******************************************************************************
+ * Miscellaneous functions
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+check_Tref_2d
+ (const struct sdis_scene* scn,
+ const double pos[2],
+ const double Tref,
+ const char* call_func_name);
+
+extern LOCAL_SYM res_T
+check_Tref_3d
+ (const struct sdis_scene* scn,
+ const double pos[3],
+ const double Tref,
+ const char* call_func_name);
+
+/*******************************************************************************
* Boundary sub-paths
******************************************************************************/
extern LOCAL_SYM res_T
diff --git a/src/sdis_heat_path_conductive.c b/src/sdis_heat_path_conductive.c
@@ -0,0 +1,128 @@
+/* Copyright (C) 2016-2023 |Méso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis_heat_path.h"
+#include "sdis_heat_path_conductive_c.h"
+#include "sdis_log.h"
+#include "sdis_medium_c.h"
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+check_solid_constant_properties
+ (struct sdis_device* dev,
+ const int evaluate_green,
+ const int use_wos_diffusion,
+ const struct solid_props* props_ref,
+ const struct solid_props* props)
+{
+ res_T res = RES_OK;
+ ASSERT(dev && props_ref && props);
+
+ if(props_ref->lambda != props->lambda) {
+ log_err(dev,
+ "%s: invalid thermal conductivity. One assumes a constant conductivity "
+ "for the whole solid.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(props_ref->rho != props->rho) {
+ log_err(dev,
+ "%s: invalid volumic mass. One assumes a constant volumic mass for "
+ "the whole solid.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(props_ref->cp != props->cp) {
+ log_err(dev,
+ "%s: invalid calorific capacity. One assumes a constant calorific "
+ "capacity for the whole solid.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if((evaluate_green || use_wos_diffusion) && props_ref->power != props->power) {
+ log_err(dev,
+ "%s: invalid variable power density. Stardis expects a constant power "
+ "density per solid when using WoS diffusion and/or green function "
+ "evaluation.", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+conductive_path_2d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_2d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_2d* T)
+{
+ res_T res = RES_OK;
+ ASSERT(ctx);
+
+ switch(ctx->diff_algo) {
+ case SDIS_DIFFUSION_DELTA_SPHERE:
+ res = conductive_path_delta_sphere_2d(scn, ctx, rwalk, rng, T);
+ break;
+ case SDIS_DIFFUSION_WOS:
+ res = conductive_path_wos_2d(scn, ctx, rwalk, rng, T);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return res;
+}
+
+res_T
+conductive_path_3d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T)
+{
+ res_T res = RES_OK;
+ ASSERT(ctx);
+
+ switch(ctx->diff_algo) {
+ case SDIS_DIFFUSION_DELTA_SPHERE:
+ res = conductive_path_delta_sphere_3d(scn, ctx, rwalk, rng, T);
+ break;
+ case SDIS_DIFFUSION_WOS:
+ res = conductive_path_wos_3d(scn, ctx, rwalk, rng, T);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return res;
+}
+
+/* Generate the conductive path functions */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_conductive_delta_sphere_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_conductive_delta_sphere_Xd.h"
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_conductive_wos_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_conductive_wos_Xd.h"
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -1,541 +0,0 @@
-/* Copyright (C) 2016-2023 |Méso|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_device_c.h"
-#include "sdis_green.h"
-#include "sdis_heat_path.h"
-#include "sdis_medium_c.h"
-#include "sdis_misc.h"
-#include "sdis_scene_c.h"
-
-#include <star/ssp.h>
-
-#include "sdis_Xd_begin.h"
-
-/*******************************************************************************
- * Non generic helper functions
- ******************************************************************************/
-#ifndef SDIS_HEAT_PATH_CONDUCTIVE_XD_H
-#define SDIS_HEAT_PATH_CONDUCTIVE_XD_H
-
-static res_T
-check_solid_constant_properties
- (struct sdis_device* dev,
- const int evaluate_green,
- const struct solid_props* props_ref,
- const struct solid_props* props)
-{
- res_T res = RES_OK;
- ASSERT(dev && props_ref && props);
-
- if(props_ref->lambda != props->lambda) {
- log_err(dev,
- "%s: invalid thermal conductivity. One assumes a constant conductivity "
- "for the whole solid.\n", FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
-
- if(props_ref->rho != props->rho) {
- log_err(dev,
- "%s: invalid volumic mass. One assumes a constant volumic mass for "
- "the whole solid.\n", FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
-
- if(props_ref->cp != props->cp) {
- log_err(dev,
- "%s: invalid calorific capacity. One assumes a constant calorific "
- "capacity for the whole solid.\n", FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
-
- if(evaluate_green && props_ref->power != props->power) {
- log_err(dev,
- "%s: invalid volumic power. When estimating the green function, a "
- "constant volumic power is assumed for the whole solid.\n",
- FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-#endif /* SDIS_HEAT_PATH_CONDUCTIVE_XD_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, FLT_MIN, 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_warn(scn->dev,
- "%s: inconsistent medium during the solid random walk at {%g, %g}.\n",
- FUNC_NAME, SPLIT2(pos));
-#else
- log_warn(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_warn(scn->dev,
- "%s: could not find a next valid conductive step at {%g, %g}.\n",
- FUNC_NAME, SPLIT2(pos));
-#else
- log_warn(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;
-}
-
-
-/*******************************************************************************
- * Handle the volumic power at a given diffusive step
- ******************************************************************************/
-struct XD(handle_volumic_power_args) {
- /* Forward/backward direction of the sampled diffusive step */
- const float* dir0;
- const float* dir1;
-
- /* Forward/backward intersections along the sampled diffusive step */
- const struct sXd(hit)* hit0;
- const struct sXd(hit)* hit1;
-
- /* Physical properties */
- double power; /* Volumic power */
- double lambda; /* Conductivity */
-
- float delta_solid; /* Challenged length of a diffusive step */
- float delta; /* Current length of the current diffusive step */
-
- size_t picard_order;
-};
-static const struct XD(handle_volumic_power_args)
-XD(HANDLE_VOLUMIC_POWER_ARGS_NULL) = {
- NULL, NULL, NULL, NULL, -1, -1, -1, -1, 0
-};
-
-static INLINE int
-XD(check_handle_volumic_power_args)
- (const struct XD(handle_volumic_power_args)* args)
-{
- ASSERT(args);
- return args
- && args->dir0
- && args->dir1
- && args->hit0
- && args->hit1
- && args->lambda >= 0
- && args->delta_solid > 0
- && args->delta >= 0
- && args->delta_solid >= 0
- && args->picard_order > 0;
-}
-
-static res_T
-XD(handle_volumic_power)
- (const struct sdis_scene* scn,
- const struct XD(handle_volumic_power_args)* args,
- double* out_power_term,
- struct XD(temperature)* T)
-{
- double power_term = 0;
- res_T res = RES_OK;
- ASSERT(scn && out_power_term && T && XD(check_handle_volumic_power_args)(args));
-
- /* No volumic power. Do nothing */
- if(args->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
-
- /* Check that picardN is not enabled when a volumic power is set since in
- * this situation the upper bound of the Monte-Carlo weight required by
- * picardN cannot be known */
- if(args->picard_order > 1) {
- log_err(scn->dev,
- "%s: invalid not null volumic power '%g' kg/m^3. Could not manage a "
- "volumic power when the picard order is not equal to 1; Picard order is "
- "currently set to %lu.\n",
- FUNC_NAME, args->power, (unsigned long)args->picard_order);
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* No forward/backward intersection along the sampled direction */
- if(SXD_HIT_NONE(args->hit0) && SXD_HIT_NONE(args->hit1)) {
- const double delta_in_meter = args->delta * scn->fp_to_meter;
- power_term = delta_in_meter * delta_in_meter / (2.0 * DIM * args->lambda);
- T->value += args->power * power_term;
-
- /* An intersection along this diffusive step is find. Use it to statically
- * correct the power term currently registered */
- } else {
- const double delta_s_adjusted = args->delta_solid * RAY_RANGE_MAX_SCALE;
- const double delta_s_in_meter = args->delta_solid * scn->fp_to_meter;
- double h;
- double h_in_meter;
- double cos_U_N;
- float N[DIM] = {0};
-
- if(args->delta == args->hit0->distance) {
- fX(normalize)(N, args->hit0->normal);
- cos_U_N = fX(dot)(args->dir0, N);
-
- } else {
- ASSERT(args->delta == args->hit1->distance);
- fX(normalize)(N, args->hit1->normal);
- cos_U_N = fX(dot)(args->dir1, N);
- }
-
- h = args->delta * fabs(cos_U_N);
- h_in_meter = h * scn->fp_to_meter;
-
- /* The regular power term */
- power_term = h_in_meter * h_in_meter / (2.0*args->lambda);
-
- /* Add the power corrective term. Be careful to use the adjusted
- * delta_solid to correctly handle the RAY_RANGE_MAX_SCALE factor in the
- * computation of the limit angle. But keep going with the unmodified
- * delta_solid in the corrective term since it was the one that was
- * "wrongly" used in the previous step and that must be corrected. */
- if(h == delta_s_adjusted) {
- power_term +=
- -(delta_s_in_meter * delta_s_in_meter) / (2.0*DIM*args->lambda);
-
- } else if(h < delta_s_adjusted) {
- const double sin_a = h / delta_s_adjusted;
-#if DIM==2
- /* tmp = sin(2a) / (PI - 2*a) */
- const double tmp = sin_a * sqrt(1 - sin_a*sin_a) / acos(sin_a);
- power_term +=
- -(delta_s_in_meter * delta_s_in_meter) / (4.0*args->lambda) * tmp;
-#else
- const double tmp = (sin_a*sin_a*sin_a - sin_a) / (1-sin_a);
- power_term +=
- (delta_s_in_meter * delta_s_in_meter) / (6.0*args->lambda) * tmp;
-#endif
- }
- T->value += args->power * power_term;
- }
-
-exit:
- *out_power_term = power_term;
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Local function
- ******************************************************************************/
-res_T
-XD(conductive_path)
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- double position_start[DIM];
- struct solid_props props_ref = SOLID_PROPS_NULL;
- double green_power_term = 0;
- struct sdis_medium* mdm;
- size_t istep = 0; /* Help for debug */
- res_T res = RES_OK;
- ASSERT(scn && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_SOLID);
- (void)ctx, (void)istep;
-
- /* Check the random walk consistency */
- 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;
- goto error;
- }
- /* Save the submitted position */
- dX(set)(position_start, rwalk->vtx.P);
-
- /* Retrieve the solid properties at the current position. Use them to verify
- * that those that are supposed to be constant by the conductive random walk
- * remain the same. Note that we take care of the same constraints on the
- * solid reinjection since once reinjected, the position of the random walk
- * is that at the beginning of the conductive random walk. Thus, after a
- * reinjection, the next line retrieves the properties of the reinjection
- * position. By comparing them to the properties along the random walk, we
- * thus verify that the properties are constant throughout the random walk
- * with respect to the properties of the reinjected position. */
- solid_get_properties(mdm, &rwalk->vtx, &props_ref);
-
- do { /* Solid random walk */
- struct XD(handle_volumic_power_args) handle_volpow_args =
- XD(HANDLE_VOLUMIC_POWER_ARGS_NULL);
- struct sXd(hit) hit0, hit1;
- struct solid_props props = SOLID_PROPS_NULL;
- double power_term = 0;
- double mu;
- float delta; /* Random walk numerical parameter */
- double delta_m;
- float dir0[DIM], dir1[DIM];
- float org[DIM];
-
- /* Fetch solid properties */
- res = solid_get_properties(mdm, &rwalk->vtx, &props);
- if(res != RES_OK) goto error;
-
- res = check_solid_constant_properties
- (scn->dev, ctx->green_path != NULL, &props_ref, &props);
- if(res != RES_OK) goto error;
-
- /* Check the limit condition
- * REVIEW Rfo: This can be a bug if the random walk comes from a boundary */
- if(props.temperature >= 0) {
- T->value += props.temperature;
- T->done = 1;
-
- if(ctx->green_path) {
- res = green_path_set_limit_vertex
- (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
-
- if(ctx->heat_path) {
- heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
- }
-
- break;
- }
-
- fX_set_dX(org, rwalk->vtx.P);
-
- /* Sample the direction to walk toward and compute the distance to travel */
- res = XD(sample_next_step_robust)(scn, mdm, rng, rwalk->vtx.P,
- (float)props.delta, dir0, dir1, &hit0, &hit1, &delta);
- if(res != RES_OK) goto error;
-
- /* Add the volumic power density to the measured temperature */
- handle_volpow_args.dir0 = dir0;
- handle_volpow_args.dir1 = dir1;
- handle_volpow_args.hit0 = &hit0;
- handle_volpow_args.hit1 = &hit1;
- handle_volpow_args.power = props.power;
- handle_volpow_args.lambda = props.lambda;
- handle_volpow_args.delta_solid = (float)props.delta;
- handle_volpow_args.delta = delta;
- handle_volpow_args.picard_order = get_picard_order(ctx);
- res = XD(handle_volumic_power)(scn, &handle_volpow_args, &power_term, T);
- if(res != RES_OK) goto error;
-
- /* Register the power term for the green function. Delay its registration
- * until the end of the conductive path, i.e. the path is valid */
- if(ctx->green_path && props.power != SDIS_VOLUMIC_POWER_NONE) {
- green_power_term += power_term;
- }
-
- /* Rewind the time */
- delta_m = delta * scn->fp_to_meter;
- mu = (2*DIM*props.lambda)/(props.rho*props.cp*delta_m*delta_m);
- res = XD(time_rewind)(mu, props.t0, rng, rwalk, ctx, T);
- if(res != RES_OK) goto error;
- if(T->done) break; /* Limit condition was reached */
-
- /* 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 {
- rwalk->hit = hit0;
- rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
- }
-
- /* Update the random walk position */
- XD(move_pos)(rwalk->vtx.P, dir0, delta);
-
- /* Register the new vertex against the heat path */
- res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
- SDIS_HEAT_VERTEX_CONDUCTION, (int)ctx->nbranchings);
- if(res != RES_OK) goto error;
-
- ++istep;
-
- /* Keep going while the solid random walk does not hit an interface */
- } while(SXD_HIT_NONE(&rwalk->hit));
-
- /* Register the power term for the green function */
- if(ctx->green_path && props_ref.power != SDIS_VOLUMIC_POWER_NONE) {
- res = green_path_add_power_term
- (ctx->green_path, rwalk->mdm, &rwalk->vtx, green_power_term);
- if(res != RES_OK) goto error;
- }
-
- T->func = XD(boundary_path);
- rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_conductive_c.h b/src/sdis_heat_path_conductive_c.h
@@ -0,0 +1,78 @@
+/* Copyright (C) 2016-2023 |Méso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_HEAT_PATH_CONDUCTIVE_C_H
+#define SDIS_HEAT_PATH_CONDUCTIVE_C_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct rwalk_2d;
+struct rwalk_3d;
+struct rwalk_context;
+struct sdis_device;
+struct sdis_scene;
+struct solid_props;
+struct ssp_rng;
+struct temperature_2d;
+struct temperature_3d;
+
+extern LOCAL_SYM res_T
+check_solid_constant_properties
+ (struct sdis_device* dev,
+ const int evaluate_green,
+ const int use_wos_diffusion,
+ const struct solid_props* props_ref,
+ const struct solid_props* props);
+
+/*******************************************************************************
+ * Conductive paths using the delta sphere algorithm
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+conductive_path_delta_sphere_2d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_2d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_2d* T);
+
+extern LOCAL_SYM res_T
+conductive_path_delta_sphere_3d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T);
+
+/*******************************************************************************
+ * Conductive paths using the walk on sphere algorithm
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+conductive_path_wos_2d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_2d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_2d* T);
+
+extern LOCAL_SYM res_T
+conductive_path_wos_3d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk_3d* rwalk,
+ struct ssp_rng* rng,
+ struct temperature_3d* T);
+
+#endif /* SDIS_HEAT_PATH_CONDUCTIVE_C_H */
diff --git a/src/sdis_heat_path_conductive_delta_sphere_Xd.h b/src/sdis_heat_path_conductive_delta_sphere_Xd.h
@@ -0,0 +1,481 @@
+/* Copyright (C) 2016-2023 |Méso|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_log.h"
+#include "sdis_green.h"
+#include "sdis_interface_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_scene_c.h"
+
+#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, FLT_MIN, 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_warn(scn->dev,
+ "%s: inconsistent medium during the solid random walk at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_warn(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_warn(scn->dev,
+ "%s: could not find a next valid conductive step at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_warn(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;
+}
+
+/*******************************************************************************
+ * Handle the volumic power at a given diffusive step
+ ******************************************************************************/
+struct XD(handle_volumic_power_args) {
+ /* Forward/backward direction of the sampled diffusive step */
+ const float* dir0;
+ const float* dir1;
+
+ /* Forward/backward intersections along the sampled diffusive step */
+ const struct sXd(hit)* hit0;
+ const struct sXd(hit)* hit1;
+
+ /* Physical properties */
+ double power; /* Volumic power */
+ double lambda; /* Conductivity */
+
+ float delta_solid; /* Challenged length of a diffusive step */
+ float delta; /* Current length of the current diffusive step */
+
+ size_t picard_order;
+};
+static const struct XD(handle_volumic_power_args)
+XD(HANDLE_VOLUMIC_POWER_ARGS_NULL) = {
+ NULL, NULL, NULL, NULL, -1, -1, -1, -1, 0
+};
+
+static INLINE int
+XD(check_handle_volumic_power_args)
+ (const struct XD(handle_volumic_power_args)* args)
+{
+ ASSERT(args);
+ return args
+ && args->dir0
+ && args->dir1
+ && args->hit0
+ && args->hit1
+ && args->lambda >= 0
+ && args->delta_solid > 0
+ && args->delta >= 0
+ && args->delta_solid >= 0
+ && args->picard_order > 0;
+}
+
+static res_T
+XD(handle_volumic_power)
+ (const struct sdis_scene* scn,
+ const struct XD(handle_volumic_power_args)* args,
+ double* out_power_term,
+ struct XD(temperature)* T)
+{
+ double power_term = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && out_power_term && T && XD(check_handle_volumic_power_args)(args));
+
+ /* No volumic power. Do nothing */
+ if(args->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
+
+ /* Check that picardN is not enabled when a volumic power is set since in
+ * this situation the upper bound of the Monte-Carlo weight required by
+ * picardN cannot be known */
+ if(args->picard_order > 1) {
+ log_err(scn->dev,
+ "%s: invalid not null volumic power '%g' kg/m^3. Could not manage a "
+ "volumic power when the picard order is not equal to 1; Picard order is "
+ "currently set to %lu.\n",
+ FUNC_NAME, args->power, (unsigned long)args->picard_order);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* No forward/backward intersection along the sampled direction */
+ if(SXD_HIT_NONE(args->hit0) && SXD_HIT_NONE(args->hit1)) {
+ const double delta_in_meter = args->delta * scn->fp_to_meter;
+ power_term = delta_in_meter * delta_in_meter / (2.0 * DIM * args->lambda);
+ T->value += args->power * power_term;
+
+ /* An intersection along this diffusive step is find. Use it to statically
+ * correct the power term currently registered */
+ } else {
+ const double delta_s_adjusted = args->delta_solid * RAY_RANGE_MAX_SCALE;
+ const double delta_s_in_meter = args->delta_solid * scn->fp_to_meter;
+ double h;
+ double h_in_meter;
+ double cos_U_N;
+ float N[DIM] = {0};
+
+ if(args->delta == args->hit0->distance) {
+ fX(normalize)(N, args->hit0->normal);
+ cos_U_N = fX(dot)(args->dir0, N);
+
+ } else {
+ ASSERT(args->delta == args->hit1->distance);
+ fX(normalize)(N, args->hit1->normal);
+ cos_U_N = fX(dot)(args->dir1, N);
+ }
+
+ h = args->delta * fabs(cos_U_N);
+ h_in_meter = h * scn->fp_to_meter;
+
+ /* The regular power term */
+ power_term = h_in_meter * h_in_meter / (2.0*args->lambda);
+
+ /* Add the power corrective term. Be careful to use the adjusted
+ * delta_solid to correctly handle the RAY_RANGE_MAX_SCALE factor in the
+ * computation of the limit angle. But keep going with the unmodified
+ * delta_solid in the corrective term since it was the one that was
+ * "wrongly" used in the previous step and that must be corrected. */
+ if(h == delta_s_adjusted) {
+ power_term +=
+ -(delta_s_in_meter * delta_s_in_meter) / (2.0*DIM*args->lambda);
+
+ } else if(h < delta_s_adjusted) {
+ const double sin_a = h / delta_s_adjusted;
+#if DIM==2
+ /* tmp = sin(2a) / (PI - 2*a) */
+ const double tmp = sin_a * sqrt(1 - sin_a*sin_a) / acos(sin_a);
+ power_term +=
+ -(delta_s_in_meter * delta_s_in_meter) / (4.0*args->lambda) * tmp;
+#else
+ const double tmp = (sin_a*sin_a*sin_a - sin_a) / (1-sin_a);
+ power_term +=
+ (delta_s_in_meter * delta_s_in_meter) / (6.0*args->lambda) * tmp;
+#endif
+ }
+ T->value += args->power * power_term;
+ }
+
+exit:
+ *out_power_term = power_term;
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+XD(conductive_path_delta_sphere)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ double position_start[DIM];
+ struct solid_props props_ref = SOLID_PROPS_NULL;
+ double green_power_term = 0;
+ struct sdis_medium* mdm;
+ size_t istep = 0; /* Help for debug */
+ res_T res = RES_OK;
+ ASSERT(scn && rwalk && rng && T);
+ ASSERT(rwalk->mdm->type == SDIS_SOLID);
+ (void)ctx, (void)istep;
+
+ /* Check the random walk consistency */
+ 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;
+ goto error;
+ }
+ /* Save the submitted position */
+ dX(set)(position_start, rwalk->vtx.P);
+
+ /* Retrieve the solid properties at the current position. Use them to verify
+ * that those that are supposed to be constant by the conductive random walk
+ * remain the same. Note that we take care of the same constraints on the
+ * solid reinjection since once reinjected, the position of the random walk
+ * is that at the beginning of the conductive random walk. Thus, after a
+ * reinjection, the next line retrieves the properties of the reinjection
+ * position. By comparing them to the properties along the random walk, we
+ * thus verify that the properties are constant throughout the random walk
+ * with respect to the properties of the reinjected position. */
+ solid_get_properties(mdm, &rwalk->vtx, &props_ref);
+
+ do { /* Solid random walk */
+ struct XD(handle_volumic_power_args) handle_volpow_args =
+ XD(HANDLE_VOLUMIC_POWER_ARGS_NULL);
+ struct sXd(hit) hit0, hit1;
+ struct solid_props props = SOLID_PROPS_NULL;
+ double power_term = 0;
+ double mu;
+ float delta; /* Random walk numerical parameter */
+ double delta_m;
+ float dir0[DIM], dir1[DIM];
+ float org[DIM];
+
+ /* Fetch solid properties */
+ res = solid_get_properties(mdm, &rwalk->vtx, &props);
+ if(res != RES_OK) goto error;
+
+ res = check_solid_constant_properties
+ (scn->dev, ctx->green_path != NULL, 0/*use WoS?*/, &props_ref, &props);
+ if(res != RES_OK) goto error;
+
+ /* Check the limit condition
+ * REVIEW Rfo: This can be a bug if the random walk comes from a boundary */
+ if(SDIS_TEMPERATURE_IS_KNOWN(props.temperature)) {
+ T->value += props.temperature;
+ T->done = 1;
+
+ if(ctx->green_path) {
+ res = green_path_set_limit_vertex
+ (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+ if(ctx->heat_path) {
+ heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
+ }
+
+ break;
+ }
+
+ fX_set_dX(org, rwalk->vtx.P);
+
+ /* Sample the direction to walk toward and compute the distance to travel */
+ res = XD(sample_next_step_robust)(scn, mdm, rng, rwalk->vtx.P,
+ (float)props.delta, dir0, dir1, &hit0, &hit1, &delta);
+ if(res != RES_OK) goto error;
+
+ /* Add the volumic power density to the measured temperature */
+ handle_volpow_args.dir0 = dir0;
+ handle_volpow_args.dir1 = dir1;
+ handle_volpow_args.hit0 = &hit0;
+ handle_volpow_args.hit1 = &hit1;
+ handle_volpow_args.power = props.power;
+ handle_volpow_args.lambda = props.lambda;
+ handle_volpow_args.delta_solid = (float)props.delta;
+ handle_volpow_args.delta = delta;
+ handle_volpow_args.picard_order = get_picard_order(ctx);
+ res = XD(handle_volumic_power)(scn, &handle_volpow_args, &power_term, T);
+ if(res != RES_OK) goto error;
+
+ /* Register the power term for the green function. Delay its registration
+ * until the end of the conductive path, i.e. the path is valid */
+ if(ctx->green_path && props.power != SDIS_VOLUMIC_POWER_NONE) {
+ green_power_term += power_term;
+ }
+
+ /* Rewind the time */
+ delta_m = delta * scn->fp_to_meter;
+ mu = (2*DIM*props.lambda)/(props.rho*props.cp*delta_m*delta_m);
+ res = XD(time_rewind)(mu, props.t0, rng, rwalk, ctx, T);
+ if(res != RES_OK) goto error;
+ if(T->done) break; /* Limit condition was reached */
+
+ /* 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 {
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ }
+
+ /* Update the random walk position */
+ XD(move_pos)(rwalk->vtx.P, dir0, delta);
+
+ /* Register the new vertex against the heat path */
+ res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
+ SDIS_HEAT_VERTEX_CONDUCTION, (int)ctx->nbranchings);
+ if(res != RES_OK) goto error;
+
+ ++istep;
+
+ /* Keep going while the solid random walk does not hit an interface */
+ } while(SXD_HIT_NONE(&rwalk->hit));
+
+ /* Register the power term for the green function */
+ if(ctx->green_path && props_ref.power != SDIS_VOLUMIC_POWER_NONE) {
+ res = green_path_add_power_term
+ (ctx->green_path, rwalk->mdm, &rwalk->vtx, green_power_term);
+ if(res != RES_OK) goto error;
+ }
+
+ T->func = XD(boundary_path);
+ rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_conductive_wos_Xd.h b/src/sdis_heat_path_conductive_wos_Xd.h
@@ -0,0 +1,669 @@
+/* Copyright (C) 2016-2023 |Méso|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_device_c.h"
+#include "sdis_heat_path_conductive_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
+
+#include <star/swf.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static res_T
+XD(check_medium_consistency)
+ (struct sdis_scene* scn,
+ const struct XD(rwalk)* rwalk)
+{
+ struct sdis_medium* mdm = NULL;
+ res_T res = RES_OK;
+ ASSERT(rwalk);
+
+ res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm);
+ if(res != RES_OK) goto error;
+
+ /* Check medium consistency */
+ if(mdm != rwalk->mdm) {
+ log_err(scn->dev,
+ "%s:%s: invalid solid walk. Unexpected medium (position: "FORMAT_VECX").\n",
+ __FILE__, FUNC_NAME, SPLITX(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(time_travel)
+ (struct sdis_scene* scn,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ const double alpha, /* Diffusivity, i.e. lambda/(rho*cp) */
+ const double t0, /* Initial time [s] */
+ double* distance, /* Displacement [m/fp_to_meter] */
+ struct XD(temperature)* T)
+{
+ double dir[DIM] = {0};
+ double dst = 0; /* Distance [m] */
+ double tau = 0; /* Time [s] */
+ double x = 0;
+ double r = 0;
+ double temperature = 0; /* [k] */
+ double time = 0; /* [s] */
+ res_T res = RES_OK;
+ ASSERT(scn && rwalk && rng && alpha > 0 && distance && T);
+
+ dst = *distance * scn->fp_to_meter;
+ ASSERT(dst >= 0);
+
+ /* No displacement => no time travel */
+ if(distance == 0) goto exit;
+
+ /* Sample x = tau*alpha/distance^2 */
+ r = ssp_rng_canonical(rng);
+ x = swf_tabulation_inverse(XD(scn->dev->H), SWF_QUADRATIC, r);
+
+ /* Retrieve the time to travel */
+ tau = x / alpha * dst * dst;
+ time = MMIN(tau, rwalk->vtx.time - t0);
+
+ /* Increment the elapsed time */
+ rwalk->elapsed_time += time;
+
+ if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */
+
+ /* Let's take a trip back in time */
+ rwalk->vtx.time = MMAX(t0, rwalk->vtx.time - tau);
+
+ /* The path does not reach the initial condition */
+ if(rwalk->vtx.time > t0) goto exit;
+
+ /* The path reaches the initial condition. Sample a distance corresponding to
+ * the travel time to the initial condition.
+ *
+ * TODO while we use the H function to sample the distance, one should use the
+ * U function. For the moment, this function is not available, hence the use
+ * of H. This is not a problem, since we currently assume that the initial
+ * condition is uniform. Position is only used for path geometry */
+ r = ssp_rng_canonical(rng);
+ x = swf_tabulation_inverse(XD(scn->dev->H), SWF_QUADRATIC, r);
+ dst = sqrt(alpha * time / x);
+ *distance = dst; /* Update travel distance */
+
+ /* Uniformly sample a direction and move along it of the distance that
+ * separate the current position ot its initial condition */
+#if DIM == 2
+ ssp_ran_circle_uniform(rng, dir, NULL);
+#else
+ ssp_ran_sphere_uniform(rng, dir, NULL);
+#endif
+ dX(muld)(dir, dir, dst);
+ dX(add)(rwalk->vtx.P, rwalk->vtx.P, dir);
+
+ /* Fetch the initial temperature */
+ temperature = medium_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) {
+ log_err(scn->dev,
+ "%s:%s: the path reaches the initial condition but the "
+ "%s temperature remains unknown -- position=%g, %g, %g\n",
+ __FILE__, FUNC_NAME,
+ medium_type_to_string(sdis_medium_get_type(rwalk->mdm)),
+ SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Update the temperature */
+ T->value += temperature;
+ T->done = 1;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#if DIM == 2
+static INLINE enum sdis_side
+compute_hit_side_2d
+ (const struct s2d_hit* hit,
+ const double pos[2]) /* Position from which intersection occurs */
+{
+ struct s2d_attrib p0, p1; /* Segment positions */
+ double v0[2] = {0}; /* Vector from segment vertex 0 to segment vertex 1 */
+ double v1[2] = {0}; /* Vector from segment vertex 0 to input position */
+ double z = 0;
+
+ /* Check pre-conditions */
+ ASSERT(hit && pos && !S2D_HIT_NONE(hit));
+
+ /* Retrieve the positions of the intersected segment */
+ S2D(segment_get_vertex_attrib(&hit->prim, 0, S2D_POSITION, &p0));
+ S2D(segment_get_vertex_attrib(&hit->prim, 1, S2D_POSITION, &p1));
+
+ v0[0] = p1.value[0] - p0.value[0];
+ v0[1] = p1.value[1] - p0.value[1];
+ v1[0] = pos[0] - p0.value[0];
+ v1[1] = pos[1] - p0.value[1];
+
+ /* Z coordinate of the cross product between v0 and v1. Its sign indicates on
+ * which side of the segment the position lies. */
+ z = d2_cross(v1, v0);
+ return z > 0 ? SDIS_FRONT : SDIS_BACK;
+}
+#endif
+
+#if DIM == 3
+static INLINE enum sdis_side
+compute_hit_side_3d
+ (const struct s3d_hit* hit,
+ const double pos[3]) /* Position from which intersection occurs */
+{
+ struct s3d_attrib v0; /* Position of the 1st triangle vertex */
+ double p[3] = {0}; /* Position of the 1st triangle vertex in double */
+ double N[3] = {0}; /* Normalized triangle normal */
+ double D = 0; /* Last parameter of the plane triangle plane equation */
+ double dst = 0; /* Distance of pos to the plane */
+
+ /* Check pre-conditions */
+ ASSERT(hit && pos && !S3D_HIT_NONE(hit));
+
+ /* Retrieve the positions of the intersected triangle */
+ S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0));
+ d3_set_f3(p, v0.value);
+
+ /* Compute the plane equation of the triangle */
+ d3_set_f3(N, hit->normal);
+ d3_normalize(N, N);
+ D = -d3_dot(N, p);
+
+ /* Calculate the distance of the input position from the plane of the triangle
+ * and use the sign to define which side of the triangle the position is on */
+ dst = d3_dot(N, pos) + D;
+ return dst > 0 ? SDIS_FRONT : SDIS_BACK;
+}
+#endif
+
+/* Verify that the submitted position is in the expected medium */
+static res_T
+XD(check_diffusion_position)
+ (struct sdis_scene* scn,
+ const struct sdis_medium* expected_medium,
+ const double pos[DIM])
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+ enum sdis_side side;
+
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ float wos_pos[DIM] = {0};
+ float wos_radius = 0;
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(scn && expected_medium && pos);
+
+ /* Look for the nearest surface within 1 mm of the position to be checked. By
+ * limiting the search radius we speed up the closest point query. If no
+ * surface is found, we assume that the position is in the intended medium.
+ * We rely on this assumption because this function is used to verify
+ * positions during diffusive random walks. Diffusion algorithms ensure that
+ * positions are in the current medium. This function is only concerned with
+ * numerical problems which, once the new position has been calculated,
+ * position the random walk beyond the medium. In other words, the path jumps
+ * a boundary that lies within the numerical imprecision of the calculation,
+ * i.e. very close to the position to be verified. So, if no surface is found
+ * close to this position, it means that there is no nearby boundary and,
+ * consequently, no numerical problem of this kind could have arisen. */
+ wos_radius = (float)(scn->fp_to_meter * 1.0e-3);
+ fX_set_dX(wos_pos, pos);
+ SXD(scene_view_closest_point(scn->sXd(view), wos_pos, wos_radius, NULL, &hit));
+ if(SXD_HIT_NONE(&hit)) goto exit;
+
+ /* Fetch interface properties and check path consistency */
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ side = XD(compute_hit_side)(&hit, pos);
+ mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
+ if(mdm != expected_medium) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(setup_hit_wos)
+ (struct sdis_scene* scn,
+ const struct sXd(hit)* hit,
+ struct XD(rwalk)* rwalk)
+{
+ /* Geometry */
+ struct sXd(primitive) prim;
+ struct sXd(attrib) attr;
+
+ /* Properties */
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+ enum sdis_side side = SDIS_SIDE_NULL__;
+
+ /* Miscellaneous */
+ double tgt[DIM] = {0}; /* Target point, i.e. hit position */
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(rwalk && hit);
+
+ /* Find intersected position */
+ SXD(scene_view_get_primitive(scn->sXd(view), hit->prim.prim_id, &prim));
+#if DIM == 2
+ SXD(primitive_get_attrib(&prim, SXD_POSITION, hit->u, &attr));
+#else
+ SXD(primitive_get_attrib(&prim, SXD_POSITION, hit->uv, &attr));
+#endif
+
+ /* Calculate on which side the intersection occurs */
+ dX_set_fX(tgt, attr.value);
+ side = XD(compute_hit_side)(hit, rwalk->vtx.P);
+
+ /* Fetch interface properties and check path consistency */
+ interf = scene_get_interface(scn, hit->prim.prim_id);
+ mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
+ if(mdm != rwalk->mdm) {
+ log_err(scn->dev,
+ "%s:%s: the conductive path has reached an invalid interface; "
+ "unexpected medium (position: "FORMAT_VECX"; side: %s).\n",
+ __FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back");
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ /* Random walk update. Do not set the medium to NULL as the intersection is
+ * found regardless of time, so the initial condition could be reached before
+ * the interface. So we can't yet assume that the random walk has left the
+ * current medium */
+ dX(set)(rwalk->vtx.P, tgt);
+ rwalk->hit = *hit;
+ rwalk->hit_side = side;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(setup_hit_rt)
+ (struct sdis_scene* scn,
+ const double pos[DIM],
+ const double dir[DIM],
+ const struct sXd(hit)* hit,
+ struct XD(rwalk)* rwalk)
+{
+ /* Properties */
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+ enum sdis_side side = SDIS_SIDE_NULL__;
+
+ /* Miscellaneous */
+ double tgt[DIM] = {0}; /* Target point, i.e. hit position */
+ double N[DIM] = {0};
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(pos && dir && rwalk && hit);
+ ASSERT(dX(is_normalized)(dir));
+
+ /* Calculate on which side the intersection occurs */
+ dX(muld)(tgt, dir, hit->distance);
+ dX(add)(tgt, tgt, pos);
+ dX_set_fX(N, hit->normal);
+ dX(normalize)(N, N);
+ side = dX(dot)(N, dir) > 0 ? SDIS_BACK : SDIS_FRONT;
+
+ /* Fetch interface properties and check path consistency */
+ interf = scene_get_interface(scn, hit->prim.prim_id);
+ mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
+ if(mdm != rwalk->mdm) {
+ log_err(scn->dev,
+ "%s:%s: the conductive path has reached an invalid interface; "
+ "unexpected medium (position: "FORMAT_VECX"; side: %s).\n",
+ __FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back");
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ /* Random walk update. Do not set the medium to NULL as the intersection is
+ * found regardless of time, so the initial condition could be reached before
+ * the interface. So we can't yet assume that the random walk has left the
+ * current medium */
+ dX(set)(rwalk->vtx.P, tgt);
+ rwalk->hit = *hit;
+ rwalk->hit_side = side;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(sample_next_position)
+ (struct sdis_scene* scn,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ double* distance) /* Displacement distance */
+{
+ /* Intersection */
+ struct sXd(hit) hit = SXD_HIT_NULL;
+
+ /* Walk on sphere */
+ double wos_distance = 0;
+ double wos_epsilon = 0;
+ float wos_pos[DIM] = {0};
+ float wos_radius = 0;
+
+ /* Miscellaneous */
+ res_T res = RES_OK;
+ ASSERT(rwalk && rng && distance);
+
+ /* Find the closest distance from the current position to the geometry */
+ wos_radius = (float)INF;
+ fX_set_dX(wos_pos, rwalk->vtx.P);
+ SXD(scene_view_closest_point(scn->sXd(view), wos_pos, wos_radius, NULL, &hit));
+ CHK(!SXD_HIT_NONE(&hit));
+ wos_distance = hit.distance;
+
+ /* The current position is in the epsilon shell:
+ * move it to the nearest interface position */
+ wos_epsilon = scn->fp_to_meter*1.e-6;
+ if(wos_distance <= wos_epsilon) {
+ res = XD(setup_hit_wos)(scn, &hit, rwalk);
+ if(res != RES_OK) goto error;
+
+ /* Uniformly sample a new position on the surrounding sphere */
+ } else {
+ double pos[DIM] = {0};
+ double dir[DIM] = {0};
+
+#if DIM == 2
+ ssp_ran_circle_uniform(rng, dir, NULL);
+#else
+ ssp_ran_sphere_uniform(rng, dir, NULL);
+#endif
+ dX(muld)(pos, dir, (double)hit.distance);
+ dX(add)(pos, pos, rwalk->vtx.P);
+
+ /* Check that the new position is in the intended medium. Please note that
+ * we do not use the scene_get_medium_in_closed_boundaries function. It uses
+ * the ray-tracing operator, which has its own numerical uncertainty that is
+ * not the same as that of the closest point operator used by this
+ * scattering algorithm. It can therefore return the expected medium,
+ * whereas the nearest point operator would return an inconsistent medium.
+ * The next diffusion step would then detect an error. This is why we use a
+ * new function based on the same geometric operator used in the present
+ * algorithm. */
+ res = XD(check_diffusion_position)(scn, rwalk->mdm, pos);
+
+ /* Diffusion position is valid => move the path to the new position */
+ if(res == RES_OK) {
+ dX(set)(rwalk->vtx.P, pos);
+
+ /* As a result, the new position is detected as being in the wrong medium.
+ * This means that there has been a numerical problem in moving the
+ * position, which has therefore jumped the solid boundary. To solve this
+ * problem, we can move the trajectory on the solid interface along the
+ * direction of displacement. Indeed, we can assume that the position we
+ * want to move to is actually inside the epsilon shell. In this case, the
+ * trajectory will be moved to this interface in the next step anyway. */
+ } else {
+ float rt_pos[DIM] = {0};
+ float rt_dir[DIM] = {0};
+ float rt_range[2] = {0, 0};
+
+ fX_set_dX(rt_pos, rwalk->vtx.P);
+ fX_set_dX(rt_dir, dir);
+ rt_range[0] = 0;
+ rt_range[1] = (float)INF;
+ SXD(scene_view_trace_ray(scn->sXd(view), rt_pos, rt_dir, rt_range, NULL, &hit));
+
+ /* An intersection should be found. If not, we can do nothing and simply
+ * reject the path.
+ *
+ * TODO: we could take the treatment of numerical problems a step further
+ * by sampling other directions and trying to move in them again. But at
+ * present, and until we have proof to the contrary, we assume that the
+ * rejection of a path should not occur, or that it will be so rare that
+ * we don't care to save it. */
+ if(SXD_HIT_NONE(&hit)) {
+ log_err(scn->dev,
+ "%s:%s: unable to find the next diffusion position "
+ "(position: "FORMAT_VECX"; direction: "FORMAT_VECX"; distance: %g\n",
+ __FILE__, FUNC_NAME, SPLITX(pos), SPLITX(dir), wos_distance);
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ res = XD(setup_hit_rt)(scn, rwalk->vtx.P, dir, &hit, rwalk);
+ if(res != RES_OK) goto error;
+ }
+ }
+
+exit:
+ *distance = hit.distance;
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(handle_volumic_power_wos)
+ (struct sdis_scene* scn,
+ const struct solid_props* props,
+ const double distance, /* [m/fp_to_meter] */
+ double* power_term,
+ struct XD(temperature)* T)
+{
+ double dst = distance * scn->fp_to_meter; /* [m] */
+ double term = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && props && distance >= 0 && power_term && T);
+
+ if(props->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
+
+ /* No displacement => no power density */
+ if(distance == 0) goto exit;
+
+ term = dst*dst / (2*DIM*props->lambda);
+ T->value += props->power * term;
+
+exit:
+ *power_term = term;
+ return res;
+}
+
+static res_T
+XD(update_green_path)
+ (struct green_path_handle* green_path,
+ struct XD(rwalk)* rwalk,
+ struct sdis_medium* mdm,
+ const struct solid_props* props,
+ const double power_term,
+ const struct XD(temperature)* T)
+{
+ res_T res = RES_OK;
+ ASSERT(mdm && props && T);
+
+ /* Is the green function estimated? */
+ if(!green_path) goto exit;
+
+ /* Save power term for green function if any */
+ if(props->power != SDIS_VOLUMIC_POWER_NONE) {
+ res = green_path_add_power_term(green_path, mdm, &rwalk->vtx, power_term);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Set the green path limit to the current position if the initial condition
+ * has been reached */
+ if(T->done) {
+ res = green_path_set_limit_vertex
+ (green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+XD(conductive_path_wos)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ /* Properties */
+ struct sdis_medium* mdm = NULL;
+ struct solid_props props_ref = SOLID_PROPS_NULL;
+ struct solid_props props = SOLID_PROPS_NULL;
+ double alpha = 0; /* diffusivity, i.e. lambda/(rho*cp) */
+
+ /* Miscellaneous */
+ size_t ndiffusion_steps = 0; /* For debug */
+ double green_power_term = 0;
+ int green = 0;
+ const int wos = 1;
+ res_T res = RES_OK;
+ (void)ctx; /* Avoid the "unused variable" warning */
+
+ /* Check pre-conditions */
+ ASSERT(scn && ctx && rwalk && rng && T);
+ ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_SOLID);
+
+ /* Is green evaluated evaluated */
+ green = ctx->green_path != NULL;
+
+ /* Keep track of the solid. After conduction, a boundary may have been
+ * reached, so the random walk medium is NULL. However, this medium is still
+ * needed to update the green path. Hence this backup */
+ mdm = rwalk->mdm;
+
+ res = XD(check_medium_consistency)(scn, rwalk);
+ if(res != RES_OK) goto error;
+
+ /* Retrieve the solid properties at the current position. Use them to verify
+ * that those that are supposed to be constant by the conductive random walk
+ * remain the same. Note that we take care of the same constraints on the
+ * solid reinjection since once reinjected, the position of the random walk
+ * is that at the beginning of the conductive random walk. Thus, after a
+ * reinjection, the next line retrieves the properties of the reinjection
+ * position. By comparing them to the properties along the random walk, we
+ * thus verify that the properties are constant throughout the random walk
+ * with respect to the properties of the reinjected position. */
+ solid_get_properties(rwalk->mdm, &rwalk->vtx, &props_ref);
+ props = props_ref;
+
+ /* The algorithm assumes that lambda, rho and cp are constants. The
+ * diffusivity of the material (alpha) can therefore be calculated once */
+ alpha = props_ref.lambda / (props_ref.rho * props_ref.cp);
+
+ /* Sample a diffusive path */
+ for(;;) {
+ double power_term = 0; /* */
+ double dst = 0; /* [m/fp_to_meter] */
+
+ /* Register the new vertex against the heat path */
+ #define REGISTER_HEAT_VERTEX { \
+ res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, \
+ SDIS_HEAT_VERTEX_CONDUCTION, (int)ctx->nbranchings); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+
+ /* The temperature is known */
+ if(SDIS_TEMPERATURE_IS_KNOWN(props.temperature)) {
+ REGISTER_HEAT_VERTEX;
+ T->value += props.temperature;
+ T->done = 1;
+ break;
+ }
+
+ /* Find the next position of the conductive path */
+ res = XD(sample_next_position)(scn, rwalk, rng, &dst);
+ if(res != RES_OK) goto error;
+
+ /* Going back in time */
+ res = XD(time_travel)(scn, rwalk, rng, alpha, props.t0, &dst, T);
+ if(res != RES_OK) goto error;
+
+ /* Add the volumic power density */
+ res = XD(handle_volumic_power_wos)(scn, &props, dst, &power_term, T);
+ if(res != RES_OK) goto error;
+
+ REGISTER_HEAT_VERTEX;
+
+ /* Accumulate the power term */
+ if(green) green_power_term += power_term;
+
+ /* The path reaches the initial condition */
+ if(T->done) {
+ T->func = NULL;
+ break;
+ }
+
+ /* The path reaches a boundary */
+ if(!SXD_HIT_NONE(&rwalk->hit)) {
+ T->func = XD(boundary_path);
+ rwalk->mdm = NULL;
+ break;
+ }
+
+ #undef REGISTER_VERTEX
+
+ /* Retreive and check solid properties at the new position */
+ res = solid_get_properties(rwalk->mdm, &rwalk->vtx, &props);
+ if(res != RES_OK) goto error;
+ res = check_solid_constant_properties(scn->dev, green, wos, &props_ref, &props);
+ if(res != RES_OK) goto error;
+
+ ++ndiffusion_steps; /* For debug */
+ }
+
+ /* Save green function data */
+ res = XD(update_green_path)
+ (ctx->green_path, rwalk, mdm, &props_ref, green_power_term, T);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -119,7 +119,7 @@ XD(handle_known_fluid_temperature)
temperature = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
/* Check if the temperature is known */
- known_temperature = temperature >= 0;
+ known_temperature = SDIS_TEMPERATURE_IS_KNOWN(temperature);
if(!known_temperature) goto exit;
T->value += temperature;
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -78,7 +78,7 @@ XD(trace_radiative_path)
#endif
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
rwalk->hit_side = SDIS_SIDE_NULL__;
- if(scn->trad.temperature >= 0) {
+ if(SDIS_TEMPERATURE_IS_KNOWN(scn->trad.temperature)) {
T->value += scn->trad.temperature;
T->done = 1;
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -120,7 +120,9 @@ interface_side_get_temperature
case SDIS_BACK: shader = &interf->shader.back; break;
default: FATAL("Unreachable code.\n");
}
- return shader->temperature ? shader->temperature(frag, interf->data) : -1;
+ return shader->temperature
+ ? shader->temperature(frag, interf->data)
+ : SDIS_TEMPERATURE_NONE;
}
static INLINE double
@@ -182,7 +184,8 @@ interface_side_get_reference_temperature
default: FATAL("Unreachable code\n"); break;
}
return shader->reference_temperature
- ? shader->reference_temperature(frag, interf->data) : -1;
+ ? shader->reference_temperature(frag, interf->data)
+ : SDIS_TEMPERATURE_NONE;
}
static INLINE int
diff --git a/src/sdis_medium.c b/src/sdis_medium.c
@@ -23,26 +23,30 @@
/*******************************************************************************
* Helper functions
******************************************************************************/
-static int
+static res_T
check_fluid_shader(const struct sdis_fluid_shader* shader)
{
- ASSERT(shader);
- return shader->calorific_capacity
- && shader->volumic_mass
- && shader->temperature
- && 0 <= shader->t0 && shader->t0 < INF;
+ if(!shader
+ || !shader->calorific_capacity
+ || !shader->volumic_mass
+ || !shader->temperature)
+ return RES_BAD_ARG;
+
+ return RES_OK;
}
-static int
+static res_T
check_solid_shader(const struct sdis_solid_shader* shader)
{
- ASSERT(shader);
- return shader->calorific_capacity
- && shader->thermal_conductivity
- && shader->volumic_mass
- && shader->delta
- && shader->temperature
- && 0 <= shader->t0 && shader->t0 < INF;
+ if(!shader
+ || !shader->calorific_capacity
+ || !shader->thermal_conductivity
+ || !shader->volumic_mass
+ || !shader->delta
+ || !shader->temperature)
+ return RES_BAD_ARG;
+
+ return RES_OK;
}
static res_T
@@ -108,14 +112,11 @@ sdis_fluid_create
struct sdis_medium* medium = NULL;
res_T res = RES_OK;
- if(!dev || !shader || !out_medium) {
- res = RES_BAD_ARG;
- goto error;
- }
+ if(!dev || !out_medium) { res = RES_BAD_ARG; goto error; }
- if(!check_fluid_shader(shader)) {
+ res = check_fluid_shader(shader);
+ if(res != RES_OK) {
log_err(dev, "%s: invalid fluid shader.\n", FUNC_NAME);
- res = RES_BAD_ARG;
goto error;
}
@@ -162,14 +163,11 @@ sdis_solid_create
struct sdis_medium* medium = NULL;
res_T res = RES_OK;
- if(!dev || !shader || !out_medium) {
- res = RES_BAD_ARG;
- goto error;
- }
+ if(!dev || !out_medium) { res = RES_BAD_ARG; goto error; }
- if(!check_solid_shader(shader)) {
+ res = check_solid_shader(shader);
+ if(res != RES_OK) {
log_err(dev, "%s: invalid solid shader.\n", FUNC_NAME);
- res = RES_BAD_ARG;
goto error;
}
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -109,17 +109,7 @@ static const struct solid_props SOLID_PROPS_NULL = SOLID_PROPS_NULL__;
******************************************************************************/
DEFINE_MDM_CHK_PROP_FUNC(fluid, calorific_capacity, 0, INF, 0, 1)
DEFINE_MDM_CHK_PROP_FUNC(fluid, volumic_mass, 0, INF, 0, 1)
-DEFINE_MDM_CHK_PROP_FUNC(fluid, temperature, 0, INF, 1, 1)
-
-static INLINE res_T
-fluid_check_t0(struct sdis_device* dev, const double t0)
-{
- if(t0 < 0) {
- log_err(dev, "invalid negative initial time '%g'.\n", t0);
- return RES_BAD_ARG;
- }
- return RES_OK;
-}
+DEFINE_MDM_CHK_PROP_FUNC(fluid, temperature, -INF, INF, 1, 1)
DEFINE_MDM_GET_PROP_FUNC(fluid, calorific_capacity)
DEFINE_MDM_GET_PROP_FUNC(fluid, volumic_mass)
@@ -151,12 +141,6 @@ fluid_check_properties
CHK_PROP(calorific_capacity, props->cp);
#undef CHK_PROP
- /* Do not check the temperature. An invalid temperature means that the
- * temperature is unknown */
-
- res = fluid_check_t0(dev, props->t0);
- if(res != RES_OK) return res;
-
return RES_OK;
}
@@ -181,17 +165,7 @@ DEFINE_MDM_CHK_PROP_FUNC(solid, thermal_conductivity, 0, INF, 0, 1)
DEFINE_MDM_CHK_PROP_FUNC(solid, volumic_mass, 0, INF, 0, 1)
DEFINE_MDM_CHK_PROP_FUNC(solid, delta, 0, INF, 0, 1)
DEFINE_MDM_CHK_PROP_FUNC(solid, volumic_power, -INF, INF, 1, 1)
-DEFINE_MDM_CHK_PROP_FUNC(solid, temperature, 0, INF, 1, 1)
-
-static INLINE res_T
-solid_check_t0(struct sdis_device* dev, const double t0)
-{
- if(t0 < 0) {
- log_err(dev, "invalid negative initial time '%g'.\n", t0);
- return RES_BAD_ARG;
- }
- return RES_OK;
-}
+DEFINE_MDM_CHK_PROP_FUNC(solid, temperature, -INF, INF, 1, 1)
DEFINE_MDM_GET_PROP_FUNC(solid, calorific_capacity)
DEFINE_MDM_GET_PROP_FUNC(solid, thermal_conductivity)
@@ -214,7 +188,6 @@ static INLINE double
solid_get_t0(const struct sdis_medium* mdm)
{
ASSERT(mdm && mdm->type == SDIS_SOLID);
- ASSERT(0 <= mdm->shader.solid.t0 && mdm->shader.solid.t0 < INF);
return mdm->shader.solid.t0;
}
@@ -242,9 +215,6 @@ solid_check_properties
/* Do not check the temperature. An invalid temperature means that the
* temperature is unknown */
- res = solid_check_t0(dev, props->t0);
- if(res != RES_OK) return res;
-
return RES_OK;
}
diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h
@@ -54,7 +54,7 @@ XD(time_rewind)
/* Fetch the initial temperature */
temperature = medium_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(temperature < 0) {
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) {
log_err(rwalk->mdm->dev, "the path reaches the limit condition but the "
"%s temperature remains unknown -- position=%g, %g, %g\n",
medium_type_to_string(sdis_medium_get_type(rwalk->mdm)),
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -20,14 +20,15 @@
* Helper functions
******************************************************************************/
static INLINE int
-check_ray_realisatio_args(const struct ray_realisation_args* args)
+check_ray_realisation_args(const struct ray_realisation_args* args)
{
return args
&& args->rng
&& args->medium
&& args->medium->type == SDIS_FLUID
&& args->time >= 0
- && args->picard_order > 0;
+ && args->picard_order > 0
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
}
/*******************************************************************************
@@ -50,7 +51,7 @@ ray_realisation_3d
struct temperature_3d T = TEMPERATURE_NULL_3d;
float dir[3];
res_T res = RES_OK;
- ASSERT(scn && weight && check_ray_realisatio_args(args));
+ ASSERT(scn && weight && check_ray_realisation_args(args));
d3_set(rwalk.vtx.P, args->position);
rwalk.vtx.time = args->time;
@@ -67,6 +68,7 @@ ray_realisation_3d
ctx.That3 = ctx.That * ctx.That2;
ctx.max_branchings = args->picard_order - 1;
ctx.irealisation = args->irealisation;
+ ctx.diff_algo = args->diff_algo;
f3_set_d3(dir, args->direction);
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -65,9 +65,18 @@ struct probe_realisation_args {
struct green_path_handle* green_path; /* May be NULL */
struct sdis_heat_path* heat_path; /* May be NULL */
size_t irealisation; /* Id of the realisation (for debug) */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define PROBE_REALISATION_ARGS_NULL__ { \
- NULL, NULL, {0,0,0}, -1, 0, NULL, NULL, SIZE_MAX \
+ NULL, /* RNG */ \
+ NULL, /* Medium */ \
+ {0,0,0}, /* Position */ \
+ -1, /* Observation time */ \
+ 0, /* Picard order */ \
+ NULL, /* Green path */ \
+ NULL, /* Heat path */ \
+ SIZE_MAX, /* Realisation ID */ \
+ SDIS_DIFFUSION_NONE /* Diffusion algorithm */ \
}
static const struct probe_realisation_args PROBE_REALISATION_ARGS_NULL =
PROBE_REALISATION_ARGS_NULL__;
@@ -97,9 +106,19 @@ struct boundary_realisation_args {
struct green_path_handle* green_path; /* May be NULL */
struct sdis_heat_path* heat_path; /* May be NULL */
size_t irealisation; /* Id of the realisation (for debug) */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define BOUNDARY_REALISATION_ARGS_NULL__ { \
- NULL, SIZE_MAX, {0,0}, -1, 0, SDIS_SIDE_NULL__, NULL, NULL, SIZE_MAX \
+ NULL, /* RNG */ \
+ SIZE_MAX, /* Primitive ID */ \
+ {0,0}, /* Parametric coordinates */ \
+ -1, /* Observation time */ \
+ 0, /* Picard order */ \
+ SDIS_SIDE_NULL__, /* Interface side */ \
+ NULL, /* Green path */ \
+ NULL, /* Heat path */ \
+ SIZE_MAX, /* Realisation ID */ \
+ SDIS_DIFFUSION_NONE /* Diffusion algorithm */ \
}
static const struct boundary_realisation_args BOUNDARY_REALISATION_ARGS_NULL =
BOUNDARY_REALISATION_ARGS_NULL__;
@@ -128,9 +147,18 @@ struct boundary_flux_realisation_args {
enum sdis_side solid_side; /* Side of the geometric primitive */
int flux_mask; /* Combination of enum flux_flag */
size_t irealisation; /* Id of the realisation (for debug) */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define BOUNDARY_FLUX_REALISATION_ARGS_NULL__ { \
- NULL, SIZE_MAX, {0,0}, -1, 0, SDIS_SIDE_NULL__, 0, SIZE_MAX \
+ NULL, /* RNG */ \
+ SIZE_MAX, /* Primitive ID */ \
+ {0,0}, /* Parametric coordinates */ \
+ -1, /* Observation time */ \
+ 0, /* Picard order */ \
+ SDIS_SIDE_NULL__, /* Interface side */ \
+ 0, /* Flux mask */ \
+ SIZE_MAX, /* Realisation ID */ \
+ SDIS_DIFFUSION_NONE /* Diffusion algorithm */ \
}
static const struct boundary_flux_realisation_args
BOUNDARY_FLUX_REALISATION_ARGS_NULL = BOUNDARY_FLUX_REALISATION_ARGS_NULL__;
@@ -159,9 +187,18 @@ struct ray_realisation_args {
size_t picard_order; /* Picard order to estimate radiative temperature */
struct sdis_heat_path* heat_path; /* May be NULL */
size_t irealisation; /* Id of the realisation (for debug) */
+ enum sdis_diffusion_algorithm diff_algo; /* Diffusion algorithm to be used */
};
#define RAY_REALISATION_ARGS_NULL__ { \
- NULL, NULL, {0,0,0}, {0,0,0}, -1, 0, NULL, SIZE_MAX \
+ NULL, /* RNG */ \
+ NULL, /* Medium */ \
+ {0,0,0}, /* Position */ \
+ {0,0,0}, /* Direction */ \
+ -1, /* Observation time */ \
+ 0, /* Picard order */ \
+ NULL, /* Heat path */ \
+ SIZE_MAX, /* Realisation ID */ \
+ SDIS_DIFFUSION_NONE /* Diffusion algorithm */ \
}
static const struct ray_realisation_args RAY_REALISATION_ARGS_NULL =
RAY_REALISATION_ARGS_NULL__;
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -39,7 +39,8 @@ check_probe_realisation_args(const struct probe_realisation_args* args)
&& args->rng
&& args->medium
&& args->time >= 0
- && args->picard_order > 0;
+ && args->picard_order > 0
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
}
static INLINE int
@@ -53,7 +54,8 @@ check_boundary_realisation_args(const struct boundary_realisation_args* args)
&& args->uv[1] <= 1
&& args->time >= 0
&& args->picard_order > 0
- && (args->side == SDIS_FRONT || args->side == SDIS_BACK);
+ && (args->side == SDIS_FRONT || args->side == SDIS_BACK)
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
}
static INLINE int
@@ -68,7 +70,8 @@ check_boundary_flux_realisation_args
&& args->uv[1] <= 1
&& args->time >= 0
&& args->picard_order > 0
- && (args->solid_side == SDIS_FRONT || args->solid_side == SDIS_BACK);
+ && (args->solid_side == SDIS_FRONT || args->solid_side == SDIS_BACK)
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
}
#endif /* SDIS_REALISATION_XD_H */
@@ -208,7 +211,7 @@ XD(probe_realisation)
/* Check the initial condition. */
rwalk.vtx.time = t0;
tmp = get_initial_temperature(args->medium, &rwalk.vtx);
- if(tmp >= 0) {
+ if(SDIS_TEMPERATURE_IS_KNOWN(tmp)) {
*weight = tmp;
goto exit;
}
@@ -234,11 +237,12 @@ XD(probe_realisation)
ctx.That3 = ctx.That * ctx.That2;
ctx.max_branchings = args->picard_order - 1;
ctx.irealisation = args->irealisation;
+ ctx.diff_algo = args->diff_algo;
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
- ASSERT(T.value >= 0);
+ ASSERT(SDIS_TEMPERATURE_IS_KNOWN(T.value));
*weight = T.value;
exit:
@@ -309,6 +313,7 @@ XD(boundary_realisation)
ctx.That3 = ctx.That * ctx.That2;
ctx.max_branchings = args->picard_order - 1;
ctx.irealisation = args->irealisation;
+ ctx.diff_algo = args->diff_algo;
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) goto error;
@@ -396,6 +401,7 @@ XD(boundary_flux_realisation)
ctx.That3 = That3; \
ctx.max_branchings = args->picard_order - 1; \
ctx.irealisation = args->irealisation; \
+ ctx.diff_algo = args->diff_algo; \
dX(set)(rwalk.vtx.P, P); \
fX(set)(rwalk.hit.normal, N); \
T = XD(TEMPERATURE_NULL); \
@@ -418,7 +424,7 @@ XD(boundary_flux_realisation)
T.func = XD(radiative_path);
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
- ASSERT(T.value >= 0);
+ ASSERT(SDIS_TEMPERATURE_IS_KNOWN(T.value));
result->Tradiative = T.value;
}
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -578,3 +578,41 @@ exit:
error:
goto exit;
}
+
+res_T
+scene_check_temperature_range(const struct sdis_scene* scn)
+{
+ res_T res = RES_OK;
+ ASSERT(scn);
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmin)) {
+ log_err(scn->dev,
+ "%s the defined minimum temperature is unknown "
+ "when it is expected to be known.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmax)) {
+ log_err(scn->dev,
+ "%s the defined maximum temperature is unknown "
+ "when it is expected to be known.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scn->tmin > scn->tmax) {
+ log_err(scn->dev,
+ "%s: defined temperature range degenerated -- [%g, %g] K\n",
+ FUNC_NAME, scn->tmin, scn->tmax);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -290,6 +290,14 @@ extern LOCAL_SYM res_T
scene_check_dimensionality_3d
(const struct sdis_scene* scn);
+/* Check that the temperature range of the scene is well defined, i.e. that the
+ * minimum and maximum temperatures are known and that they define a valid
+ * range. If this is not the case, the function displays an error message and
+ * returns RES_BAD_ARG */
+extern LOCAL_SYM res_T
+scene_check_temperature_range
+ (const struct sdis_scene* scn);
+
static INLINE void
scene_get_enclosure_ids
(const struct sdis_scene* scn,
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -23,6 +23,7 @@
#include "sdis_misc.h"
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include "sdis_heat_path_boundary_c.h" /* check_Tref_<2d|3d> */
#include <rsys/clock_time.h>
#include <star/ssp.h>
@@ -78,6 +79,11 @@ check_solve_boundary_args(const struct sdis_solve_boundary_args* args)
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -405,6 +411,7 @@ XD(solve_boundary)
realis_args.green_path = pgreen_path;
realis_args.heat_path = pheat_path;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
realis_args.uv[0] = uv[0];
#if SDIS_XD_DIMENSION == 3
realis_args.uv[1] = uv[1];
@@ -778,13 +785,20 @@ XD(solve_boundary_flux)
/* Fetch interface parameters */
epsilon = interface_side_get_emissivity(interf, &frag);
- Tref = interface_side_get_reference_temperature(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
- hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ Tref = interface_side_get_reference_temperature(interf, &frag);
+ if(epsilon <= 0) {
+ hr = 0;
+ } else {
+ res_local = XD(check_Tref)(scn, frag.P, Tref, FUNC_NAME);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, &res_local); continue; }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ }
+
frag.side = solid_side;
imposed_flux = interface_side_get_flux(interf, &frag);
imposed_temp = interface_side_get_temperature(interf, &frag);
- if(imposed_temp >= 0) {
+ if(SDIS_TEMPERATURE_IS_KNOWN(imposed_temp)) {
/* Flux computation on T boundaries is not supported yet */
log_err(scn->dev, "%s: Attempt to compute a flux at a Dirichlet boundary "
"(not available yet).\n", FUNC_NAME);
@@ -805,6 +819,7 @@ XD(solve_boundary_flux)
realis_args.solid_side = solid_side;
realis_args.flux_mask = flux_mask;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
realis_args.uv[0] = uv[0];
#if SDIS_XD_DIMENSION == 3
realis_args.uv[1] = uv[1];
@@ -820,10 +835,9 @@ XD(solve_boundary_flux)
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
/* Convective flux from fluid to solid */
- const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ const double w_conv = hc > 0 ? hc * (result.Tfluid - result.Tboundary) : 0;
/* Radiative flux from ambient to solid */
- const double w_rad = (result.Tradiative < 0) ?
- 0 : hr * (result.Tradiative - result.Tboundary);
+ const double w_rad = hr > 0 ? hr * (result.Tradiative - result.Tboundary) : 0;
/* Imposed flux that goes _into_ the solid */
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
const double w_total = w_conv + w_rad + w_imp;
diff --git a/src/sdis_solve_camera.c b/src/sdis_solve_camera.c
@@ -72,6 +72,11 @@ check_solve_camera_args(const struct sdis_solve_camera_args* args)
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -87,6 +92,7 @@ solve_pixel
const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
const size_t picard_order,
+ const enum sdis_diffusion_algorithm diff_algo,
struct sdis_estimator* estimator,
struct pixel* pixel)
{
@@ -135,6 +141,7 @@ solve_pixel
realis_args.picard_order = picard_order;
realis_args.heat_path = pheat_path;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = diff_algo;
d3_set(realis_args.position, ray_pos);
d3_set(realis_args.direction, ray_dir);
res_simul = ray_realisation_3d(scn, &realis_args, &w);
@@ -197,6 +204,7 @@ solve_tile
const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
const size_t picard_order,
+ const enum sdis_diffusion_algorithm diff_algo,
struct sdis_estimator_buffer* buf,
struct tile* tile)
{
@@ -234,7 +242,7 @@ solve_tile
}
res = solve_pixel
(scn, rng, mdm, cam, time_range, ipix_image, spp, register_paths, pix_sz,
- picard_order, estimator, pixel);
+ picard_order, diff_algo, estimator, pixel);
if(res != RES_OK) goto error;
}
@@ -651,7 +659,8 @@ sdis_solve_camera
/* Draw the tile */
res_local = solve_tile
(scn, rng, medium, args->cam, args->time_range, tile_org, tile_sz,
- args->spp, register_paths, pix_sz, args->picard_order, buf, tile);
+ args->spp, register_paths, pix_sz, args->picard_order, args->diff_algo,
+ buf, tile);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -172,6 +172,11 @@ check_solve_medium_args(const struct sdis_solve_medium_args* args)
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -436,6 +441,7 @@ XD(solve_medium)
realis_args.green_path = pgreen_path;
realis_args.heat_path = pheat_path;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
dX(set)(realis_args.position, pos);
res_simul = XD(probe_realisation)(scn, &realis_args, &weight);
if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -64,6 +64,11 @@ check_solve_probe_args(const struct sdis_solve_probe_args* args)
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -143,6 +148,7 @@ XD(solve_one_probe)
realis_args.time = time;
realis_args.picard_order = args->picard_order;
realis_args.irealisation = irealisation;
+ realis_args.diff_algo = args->diff_algo;
dX(set)(realis_args.position, args->position);
res = XD(probe_realisation)(scn, &realis_args, &w);
if(res != RES_OK) goto error;
@@ -327,6 +333,7 @@ XD(solve_probe)
realis_args.green_path = pgreen_path;
realis_args.heat_path = pheat_path;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
dX(set)(realis_args.position, args->position);
res_simul = XD(probe_realisation)(scn, &realis_args, &w);
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -23,6 +23,7 @@
#include "sdis_misc.h"
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include "sdis_heat_path_boundary_c.h" /* check_Tref_<2d|3d> */
#include <rsys/clock_time.h>
#include <star/ssp.h>
@@ -71,6 +72,11 @@ check_solve_probe_boundary_args
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -138,6 +144,11 @@ check_solve_probe_boundary_flux_args
return RES_BAD_ARG;
}
+ /* Check the diffusion algorithm */
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
return RES_OK;
}
@@ -178,6 +189,7 @@ XD(solve_one_probe_boundary)
realis_args.picard_order = args->picard_order;
realis_args.side = args->side;
realis_args.irealisation = irealisation;
+ realis_args.diff_algo = args->diff_algo;
realis_args.uv[0] = args->uv[0];
#if SDIS_XD_DIMENSION == 3
realis_args.uv[1] = args->uv[1];
@@ -363,6 +375,7 @@ XD(solve_probe_boundary)
realis_args.green_path = pgreen_path;
realis_args.heat_path = pheat_path;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
realis_args.uv[0] = args->uv[0];
#if SDIS_XD_DIMENSION == 3
realis_args.uv[1] = args->uv[1];
@@ -840,6 +853,7 @@ XD(solve_probe_boundary_flux)
BOUNDARY_FLUX_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
+ struct sdis_interface_fragment frag_local = frag;
struct ssp_rng* rng = per_thread_rng[ithread];
struct accum* acc_temp = &per_thread_acc_tp[ithread];
struct accum* acc_time = &per_thread_acc_ti[ithread];
@@ -850,9 +864,10 @@ XD(solve_probe_boundary_flux)
double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
- double Tref = -1;
+ double Tref = SDIS_TEMPERATURE_NONE;
size_t n;
int pcent;
+ res_T res_local = RES_OK;
res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
@@ -863,16 +878,23 @@ XD(solve_probe_boundary_flux)
time = sample_time(rng, args->time_range);
/* Compute hr and hc */
- frag.time = time;
- frag.side = fluid_side;
- epsilon = interface_side_get_emissivity(interf, &frag);
- Tref = interface_side_get_reference_temperature(interf, &frag);
- hc = interface_get_convection_coef(interf, &frag);
- hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
- frag.side = solid_side;
- imposed_flux = interface_side_get_flux(interf, &frag);
- imposed_temp = interface_side_get_temperature(interf, &frag);
- if(imposed_temp >= 0) {
+ frag_local.time = time;
+ frag_local.side = fluid_side;
+ epsilon = interface_side_get_emissivity(interf, &frag_local);
+ Tref = interface_side_get_reference_temperature(interf, &frag_local);
+ hc = interface_get_convection_coef(interf, &frag_local);
+ if(epsilon <= 0) {
+ hr = 0;
+ } else {
+ res_local = XD(check_Tref)(scn, frag_local.P, Tref, FUNC_NAME);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, &res_local); continue; }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ }
+
+ frag_local.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag_local);
+ imposed_temp = interface_side_get_temperature(interf, &frag_local);
+ if(SDIS_TEMPERATURE_IS_KNOWN(imposed_temp)) {
/* Flux computation on T boundaries is not supported yet */
log_err(scn->dev,
"%s: Attempt to compute a flux at a Dirichlet boundary "
@@ -894,6 +916,7 @@ XD(solve_probe_boundary_flux)
realis_args.solid_side = solid_side;
realis_args.flux_mask = flux_mask;
realis_args.irealisation = (size_t)irealisation;
+ realis_args.diff_algo = args->diff_algo;
realis_args.uv[0] = args->uv[0];
#if SDIS_XD_DIMENSION == 3
realis_args.uv[1] = args->uv[1];
@@ -909,10 +932,9 @@ XD(solve_probe_boundary_flux)
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
/* Convective flux from fluid to solid */
- const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ const double w_conv = hc > 0 ? hc * (result.Tfluid - result.Tboundary) : 0;
/* Radiative flux from ambient to solid */
- const double w_rad = (result.Tradiative < 0) ?
- 0 : hr * (result.Tradiative - result.Tboundary);
+ const double w_rad = hr > 0 ? hr * (result.Tradiative - result.Tboundary) : 0;
/* Imposed flux that goes _into_ the solid */
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
/* Total flux */
diff --git a/src/test_sdis_compute_power.c b/src/test_sdis_compute_power.c
@@ -19,7 +19,6 @@
#include <rsys/stretchy_array.h>
#include <star/s3dut.h>
-#define UNKOWN_TEMPERATURE -1
#define N 100000ul /* #realisations */
#define POWER0 10
#define POWER1 5
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -19,8 +19,6 @@
#include <rsys/math.h>
#include <star/ssp.h>
-#define UNKNOWN_TEMPERATURE -1
-
/* The scene is composed of a solid cube whose temperature is unknown. The cube
* faces on +/-X are in contact with a fluid and their convection coefficient
* is null while their emissivity is 1. The left and right fluids are enclosed
@@ -134,7 +132,7 @@ static double
temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
CHK(vtx != NULL); (void)data;
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -179,7 +177,7 @@ solid_get_temperature
CHK(data != NULL);
t0 = ((const struct solid*)sdis_data_cget(data))->t0;
if(vtx->time > t0) {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
} else {
return ((const struct solid*)sdis_data_cget(data))->initial_temperature;
}
@@ -395,7 +393,7 @@ main(int argc, char** argv)
OK(sdis_data_ref_put(data));
/* Create the interface that forces to keep in conduction */
- interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.temperature = SDIS_TEMPERATURE_NONE;
interf.convection_coef = -1;
interf.emissivity = -1;
interf.specular_fraction = -1;
@@ -403,7 +401,7 @@ main(int argc, char** argv)
create_interface(dev, solid, solid2, &interf, interfaces+0);
/* Create the interface that emits radiative heat from the solid */
- interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.temperature = SDIS_TEMPERATURE_NONE;
interf.convection_coef = 0;
interf.emissivity = emissivity;
interf.specular_fraction = 1;
@@ -411,7 +409,7 @@ main(int argc, char** argv)
create_interface(dev, solid, fluid, &interf, interfaces+1);
/* Create the interface that forces the radiative heat to bounce */
- interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.temperature = SDIS_TEMPERATURE_NONE;
interf.convection_coef = 0;
interf.emissivity = 0;
interf.specular_fraction = 1;
@@ -483,7 +481,7 @@ main(int argc, char** argv)
struct sdis_estimator* estimator2;
struct sdis_green_function* green;
struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
- double ref = -1;
+ double ref = SDIS_TEMPERATURE_NONE;
size_t nreals = 0;
size_t nfails = 0;
const size_t N = 10000;
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -18,8 +18,6 @@
#include <star/ssp.h>
-#define UNKNOWN_TEMPERATURE -1
-
/* The scene is composed of a solid square whose temperature is unknown. The
* square segments on +/-X are in contact with a fluid and their convection
* coefficient is null while their emissivity is 1. The left and right fluids
@@ -113,7 +111,7 @@ static double
temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
CHK(vtx != NULL); (void)data;
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -163,7 +161,7 @@ struct interfac {
};
static const struct interfac INTERFACE_NULL = {
- 0, {-1, -1, -1, -1}, {-1, -1, -1, -1}
+ 0, {SDIS_TEMPERATURE_NONE, -1, -1, -1}, {-1, -1, -1, -1}
};
static double
@@ -171,7 +169,7 @@ interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
const struct interfac* interf;
- double T = -1;
+ double T = SDIS_TEMPERATURE_NONE;
CHK(data != NULL && frag != NULL);
interf = sdis_data_cget(data);
switch(frag->side) {
@@ -229,7 +227,7 @@ interface_get_reference_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
const struct interfac* interf;
- double T = -1;
+ double T = SDIS_TEMPERATURE_NONE;
CHK(data != NULL && frag != NULL);
interf = sdis_data_cget(data);
switch(frag->side) {
@@ -398,7 +396,7 @@ main(int argc, char** argv)
/* Create the interface that emits radiative heat from the solid */
interf = INTERFACE_NULL;
- interf.back.temperature = UNKNOWN_TEMPERATURE;
+ interf.back.temperature = SDIS_TEMPERATURE_NONE;
interf.back.emissivity = emissivity;
interf.back.specular_fraction = -1; /* Should not be fetched */
interf.back.reference_temperature = Tref;
@@ -406,7 +404,7 @@ main(int argc, char** argv)
/* Create the interface that forces the radiative heat to bounce */
interf = INTERFACE_NULL;
- interf.front.temperature = UNKNOWN_TEMPERATURE;
+ interf.front.temperature = SDIS_TEMPERATURE_NONE;
interf.front.emissivity = 0;
interf.front.specular_fraction = 1;
interf.front.reference_temperature = Tref;
diff --git a/src/test_sdis_contact_resistance.c b/src/test_sdis_contact_resistance.c
@@ -52,7 +52,6 @@
* (0,0,0)///x=X0///
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 10000 /* #realisations */
#define T0 0.0
@@ -80,7 +79,7 @@ fluid_get_temperature
{
(void)data;
CHK(vtx);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -120,7 +119,7 @@ solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
CHK(vtx && data);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
/*******************************************************************************
@@ -314,7 +313,7 @@ main(int argc, char** argv)
/* Create the adiabatic interfaces */
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1));
OK(sdis_interface_create
@@ -342,7 +341,7 @@ main(int argc, char** argv)
interf_shader.thermal_contact_resistance = interface_get_contact_resistance;
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid1, solid2, &interf_shader, data, &interf_R));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_contact_resistance_2.c b/src/test_sdis_contact_resistance_2.c
@@ -45,7 +45,6 @@
* (0,0,0)///x=X0///
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 10000 /* #realisations */
#define T0 0.0
@@ -73,7 +72,7 @@ fluid_get_temperature
{
(void)data;
CHK(vtx);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -113,7 +112,7 @@ solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
CHK(vtx && data);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
/*******************************************************************************
@@ -408,7 +407,7 @@ main(int argc, char** argv)
/* Create the adiabatic interfaces */
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1));
OK(sdis_interface_create
@@ -436,7 +435,7 @@ main(int argc, char** argv)
interf_shader.thermal_contact_resistance = interface_get_contact_resistance;
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid1, solid2, &interf_shader, data, &interf_R));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -48,7 +48,6 @@
* (0,0,0)
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 100000 /* #realisations */
#define Tf_0 280.0
@@ -73,9 +72,9 @@ fluid_get_temperature
{
CHK(vtx != NULL);
if(*((int*)sdis_data_cget(is_stationary))) {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
} else {
- return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+ return vtx->time <= 0 ? Tf_0 : SDIS_TEMPERATURE_NONE;
}
}
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -48,7 +48,6 @@
* (0,0,0)
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 100000 /* #realisations */
#define Tf_0 280.0
@@ -80,9 +79,9 @@ fluid_get_temperature
CHK(vtx != NULL);
CHK(is_stationary != NULL);
if(*((int*)sdis_data_cget(is_stationary))) {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
} else {
- return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+ return vtx->time <= 0 ? Tf_0 : SDIS_TEMPERATURE_NONE;
}
}
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -92,7 +92,7 @@ MEDIUM_PROP(solid, temperature, 310) /* [K] */
MEDIUM_PROP(solid, delta, 1.0/20.0) /* [m] */
MEDIUM_PROP(fluid, calorific_capacity, 2.0) /* [J/K/Kg] */
MEDIUM_PROP(fluid, volumic_mass, 25.0) /* |kg/m^3] */
-MEDIUM_PROP(fluid, temperature, -1/*<=> unknown*/) /* [K] */
+MEDIUM_PROP(fluid, temperature, SDIS_TEMPERATURE_NONE) /* [K] */
#undef MEDIUM_PROP
static struct sdis_medium*
@@ -138,7 +138,7 @@ interface_get_temperature
struct sdis_data* data)
{
(void)frag, (void)data;/* Avoid the "unused variable" warning */
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
diff --git a/src/test_sdis_enclosure_limit_conditions.c b/src/test_sdis_enclosure_limit_conditions.c
@@ -39,7 +39,6 @@
#define CONVECTION_COEF 10
#define DELTA 0.00625
-#define UNKNOWN_TEMPERATURE -1
#define NREALISATIONS 10000
#define Text 360.0
@@ -57,7 +56,7 @@ DEFINE_MEDIUM_GETTER(solid_get_calorific_capacity, 1)
DEFINE_MEDIUM_GETTER(solid_get_thermal_conductivity, 1)
DEFINE_MEDIUM_GETTER(solid_get_volumic_mass, 1)
DEFINE_MEDIUM_GETTER(solid_get_delta, DELTA)
-DEFINE_MEDIUM_GETTER(solid_get_temperature, UNKNOWN_TEMPERATURE)
+DEFINE_MEDIUM_GETTER(solid_get_temperature, SDIS_TEMPERATURE_NONE)
DEFINE_MEDIUM_GETTER(fluid_get_temperature, *((double*)sdis_data_cget(data)))
#undef DEFINE_MEDIUM_GETTER
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -138,7 +138,7 @@ MEDIUM_PROP(medium, volumic_mass, 1700) /* [kj/m^3] */
MEDIUM_PROP(medium, calorific_capacity, 800) /* [J/K/Kg] */
MEDIUM_PROP(solid, thermal_conductivity, 1.15) /* [W/m/K] */
MEDIUM_PROP(solid, delta, 0.1/20.0) /* [m] */
-MEDIUM_PROP(solid, temperature, -1/*<=> unknown*/) /* [K] */
+MEDIUM_PROP(solid, temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
MEDIUM_PROP(fluid, temperature, T_FLUID) /* [K] */
#undef MEDIUM_PROP
@@ -188,7 +188,7 @@ struct interface {
(void)frag, (void)data; /* Avoid the "unused variable" warning */ \
return Val; \
}
-INTERF_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+INTERF_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
INTERF_PROP(reference_temperature, T_REF) /* [K] */
#undef INTERF_PROP
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -43,7 +43,6 @@
* (0,0,0) /////
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 10000
#define PHI 10.0
@@ -96,7 +95,7 @@ solid_get_temperature
(void)data;
CHK(vtx != NULL);
if(vtx->time > 0)
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
else
return T0;
}
@@ -146,7 +145,7 @@ solve
struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
size_t nreals;
size_t nfails;
- double ref = -1;
+ double ref = SDIS_TEMPERATURE_NONE;
const int nsimuls = 4;
int isimul;
enum sdis_scene_dimension dim;
@@ -392,7 +391,7 @@ main(int argc, char** argv)
/* Create the adiabatic interface */
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
interf_props->phi = 0;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_adiabatic));
@@ -409,7 +408,7 @@ main(int argc, char** argv)
/* Create the PHI interface */
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
interf_props->phi = PHI;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_phi));
diff --git a/src/test_sdis_flux2.c b/src/test_sdis_flux2.c
@@ -16,8 +16,6 @@
#include "sdis.h"
#include "test_sdis_utils.h"
-#define UNKNOWN_TEMPERATURE -1
-
/* This test consists in solving the temperature profile in a solid slab
* surrounded by two different convective and radiative temperatures. The
* conductivity of the solid material is known, as well as its thickness. A net
@@ -175,7 +173,7 @@ solid_get_temperature
CHK(vtx && solid);
if(vtx->time > 0) {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
} else {
/* The initial temperature is a linear profile between T1 and T2, where T1
* and T2 are the temperature on the left and right slab boundary,
@@ -489,8 +487,8 @@ main(int argc, char** argv)
interf_props.h = 0;
interf_props.emissivity = 0;
interf_props.phi = SDIS_FLUX_NONE;
- interf_props.temperature = UNKNOWN_TEMPERATURE;
- interf_props.Tref = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
+ interf_props.Tref = SDIS_TEMPERATURE_NONE;
create_interface(dev, solid, dummy, &interf_props, &interfaces[ADIABATIC]);
/* Interfaces with a fixed temperature */
@@ -506,7 +504,7 @@ main(int argc, char** argv)
interf_props.h = 2;
interf_props.emissivity = 1;
interf_props.phi = 10000;
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.Tref = 300;
create_interface
(dev, solid, fluid1, &interf_props, &interfaces[SOLID_FLUID_WITH_FLUX]);
@@ -514,7 +512,7 @@ main(int argc, char** argv)
interf_props.h = 8;
interf_props.emissivity = 1;
interf_props.phi = SDIS_FLUX_NONE;
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.Tref = 300;
create_interface
(dev, solid, fluid2, &interf_props, &interfaces[SOLID_FLUID]);
diff --git a/src/test_sdis_flux_with_h.c b/src/test_sdis_flux_with_h.c
@@ -107,7 +107,7 @@ solid_get_temperature
struct sdis_data* data)
{
(void)data, (void)vtx;
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c
@@ -63,10 +63,12 @@ main(int argc, char** argv)
BA(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
fluid_shader.temperature = DUMMY_FLUID_SHADER.temperature;
- fluid_shader.t0 = -1;
- BA(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+ fluid_shader.t0 = -INF;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+ OK(sdis_medium_ref_put(fluid));
fluid_shader.t0 = INF;
- BA(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+ OK(sdis_medium_ref_put(fluid));
fluid_shader.t0 = DUMMY_FLUID_SHADER.t0;
BA(sdis_fluid_create(dev, &SDIS_FLUID_SHADER_NULL, NULL, &fluid));
@@ -106,10 +108,12 @@ main(int argc, char** argv)
BA(sdis_solid_create(dev, &solid_shader, NULL, &solid));
solid_shader.temperature = DUMMY_SOLID_SHADER.temperature;
- solid_shader.t0 = -1;
- BA(sdis_solid_create(dev, &solid_shader, NULL, &solid));
+ solid_shader.t0 = -INF;
+ OK(sdis_solid_create(dev, &solid_shader, NULL, &solid));
+ OK(sdis_medium_ref_put(solid));
solid_shader.t0 = INF;
- BA(sdis_solid_create(dev, &solid_shader, NULL, &solid));
+ OK(sdis_solid_create(dev, &solid_shader, NULL, &solid));
+ OK(sdis_medium_ref_put(solid));
solid_shader.t0 = DUMMY_SOLID_SHADER.t0;
OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c
@@ -18,7 +18,6 @@
#include <string.h>
-#define UNKNOWN_TEMPERATURE -1
#define N 10000
/* This test consists in solving the stationary temperature profile in a solid
@@ -234,7 +233,7 @@ solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
CHK(vtx && data);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -607,15 +606,15 @@ main(int argc, char** argv)
create_fluid(dev, &fluid);
/* Create the adiabatic interface for the solid */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.h = -1;
interf_props.emissivity = -1;
interf_props.specular_fraction = -1;
- interf_props.Tref = UNKNOWN_TEMPERATURE;
+ interf_props.Tref = SDIS_TEMPERATURE_NONE;
create_interface(dev, solid, dummy, &interf_props, interfaces+ADIABATIC);
/* Create the interface between the solid and the fluid */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.h = 0;
interf_props.emissivity = 1;
interf_props.specular_fraction = 0;
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -383,8 +383,16 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(NULL, &trad));
OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
- CHK(trad.temperature == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature);
- CHK(trad.reference == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference);
+ if(SDIS_TEMPERATURE_IS_KNOWN(trad.temperature)) {
+ CHK(trad.temperature == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature));
+ }
+ if(SDIS_TEMPERATURE_IS_KNOWN(trad.reference)) {
+ CHK(trad.reference == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference));
+ }
trad.temperature = 100;
trad.reference = 110;
@@ -399,8 +407,16 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_get_temperature_range(scn, NULL));
BA(sdis_scene_get_temperature_range(NULL, t_range));
OK(sdis_scene_get_temperature_range(scn, t_range));
- CHK(t_range[0] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]);
- CHK(t_range[1] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[1]);
+ if(SDIS_TEMPERATURE_IS_KNOWN(t_range[0])) {
+ CHK(t_range[0] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]));
+ }
+ if(SDIS_TEMPERATURE_IS_KNOWN(t_range[1])) {
+ CHK(t_range[1] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[1]);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]));
+ }
t_range[0] = 1;
t_range[1] = 100;
@@ -408,8 +424,8 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_set_temperature_range(NULL, t_range));
BA(sdis_scene_set_temperature_range(scn, NULL));
OK(sdis_scene_set_temperature_range(scn, t_range));
- t_range[0] = -1;
- t_range[1] = -1;
+ t_range[0] = SDIS_TEMPERATURE_NONE;
+ t_range[1] = SDIS_TEMPERATURE_NONE;
OK(sdis_scene_get_temperature_range(scn, t_range));
CHK(t_range[0] == 1);
CHK(t_range[1] == 100);
diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c
@@ -24,7 +24,6 @@
#define Hcoef 1.0 /* Convection coefficient */
#define Pw 10000.0 /* Volumetric power */
#define Nreals 10000 /* #realisations */
-#define UNKNOWN -1
/*******************************************************************************
* Helper functions
@@ -160,7 +159,7 @@ interface_get_temperature
interf = sdis_data_cget(data);
switch(interf->profile) {
case PROFILE_UNKNOWN:
- temperature = UNKNOWN;
+ temperature = SDIS_TEMPERATURE_NONE;
break;
case PROFILE_VOLUMETRIC_POWER:
temperature = volumetric_temperature(frag->P, interf->upper);
@@ -305,7 +304,7 @@ main(int argc, char** argv)
solid_param->lambda = 10;
solid_param->cp = 1.0;
solid_param->rho = 1.0;
- solid_param->temperature = -1;
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
solid_param->power = SDIS_VOLUMIC_POWER_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -42,7 +42,6 @@
* (0,0,0) /////
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 10000 /* #realisations */
#define N_dump 10 /* #dumped paths */
@@ -108,7 +107,7 @@ solid_get_temperature
{
(void)data;
CHK(vtx != NULL);
- if(vtx->time > 0) return UNKNOWN_TEMPERATURE;
+ if(vtx->time > 0) return SDIS_TEMPERATURE_NONE;
return Tf;
}
@@ -237,7 +236,7 @@ main(int argc, char** argv)
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
interf_props->hc = 0;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_adiabatic));
OK(sdis_data_ref_put(data));
@@ -255,7 +254,7 @@ main(int argc, char** argv)
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
interf_props->hc = H;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_H));
OK(sdis_data_ref_put(data));
@@ -419,7 +418,7 @@ main(int argc, char** argv)
}
/* The external fluid cannot have an unknown temperature */
- fluid_param->temperature = UNKNOWN_TEMPERATURE;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(SOLVE(box_scn, &probe_args, &estimator));
fluid_param->temperature = Tf;
@@ -445,7 +444,7 @@ main(int argc, char** argv)
}
/* The external fluid cannot have an unknown temperature */
- fluid_param->temperature = UNKNOWN_TEMPERATURE;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(SOLVE(square_scn, &probe_args, &estimator));
fluid_param->temperature = Tf;
@@ -606,7 +605,7 @@ main(int argc, char** argv)
bound_args.register_paths = SDIS_HEAT_PATH_ALL;
/* Check simulation error handling when paths are registered */
- fluid_param->temperature = UNKNOWN_TEMPERATURE;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(SOLVE(box_scn, &bound_args, &estimator));
/* Dump path */
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -62,7 +62,6 @@
* (0,0)///////
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 100000 /* #realisations */
#define Tf 300.0
@@ -133,7 +132,7 @@ solid_get_temperature
{
(void) data;
CHK(vtx != NULL);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
/*******************************************************************************
@@ -289,7 +288,7 @@ main(int argc, char** argv)
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
interf_props->hc = 0;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
interf_props->emissivity = 0;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_adiabatic));
@@ -313,7 +312,7 @@ main(int argc, char** argv)
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
interf_props->hc = H;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
interf_props->emissivity = EPSILON;
interf_props->reference_temperature = Tref;
interf_shader.back.emissivity = interface_get_emissivity;
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -28,7 +28,6 @@
#include <string.h>
-#define UNKOWN_TEMPERATURE -1
#define IMG_WIDTH 157
#define IMG_HEIGHT 53
#define SPP 30 /* #Samples per pixel, i.e. #realisations per pixel */
@@ -152,7 +151,7 @@ struct fluid {
double rho;
double temperature;
};
-static const struct fluid FLUID_NULL = {0, 0, UNKOWN_TEMPERATURE};
+static const struct fluid FLUID_NULL = {0, 0, SDIS_TEMPERATURE_NONE};
static double
fluid_get_calorific_capacity
@@ -188,7 +187,7 @@ struct solid {
double delta;
double temperature;
};
-static const struct solid SOLID_NULL = {0, 0, 0, 0, UNKOWN_TEMPERATURE};
+static const struct solid SOLID_NULL = {0, 0, 0, 0, SDIS_TEMPERATURE_NONE};
static double
solid_get_calorific_capacity
@@ -241,7 +240,7 @@ struct interf {
double reference_temperature;
};
static const struct interf INTERF_NULL = {
- 0, 0, 0, UNKOWN_TEMPERATURE, UNKOWN_TEMPERATURE
+ 0, 0, 0, SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE
};
static double
@@ -595,21 +594,21 @@ main(int argc, char** argv)
solid_param.lambda = 0.1;
solid_param.rho = 1.0;
solid_param.delta = 1.0/20.0;
- solid_param.temperature = UNKOWN_TEMPERATURE;
+ solid_param.temperature = SDIS_TEMPERATURE_NONE;
create_solid(dev, &solid_param, &solid);
/* Create the fluid0/solid interface */
interface_param.hc = 1;
interface_param.epsilon = 0;
interface_param.specular_fraction = 0;
- interface_param.temperature = UNKOWN_TEMPERATURE;
+ interface_param.temperature = SDIS_TEMPERATURE_NONE;
create_interface(dev, solid, fluid0, &interface_param, &interf0);
/* Create the fluid1/solid interface */
interface_param.hc = 0.1;
interface_param.epsilon = 1;
interface_param.specular_fraction = 1;
- interface_param.temperature = UNKOWN_TEMPERATURE;
+ interface_param.temperature = SDIS_TEMPERATURE_NONE;
interface_param.reference_temperature = 300;
create_interface(dev, fluid1, solid, &interface_param, &interf1);
@@ -674,7 +673,7 @@ main(int argc, char** argv)
BA(sdis_solve_camera(scn, &solve_args, &buf));
solve_args.cam = cam;
OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
- trad.temperature = -1;
+ trad.temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
BA(sdis_solve_camera(scn, &solve_args, &buf));
trad.temperature = 300;
@@ -767,7 +766,7 @@ main(int argc, char** argv)
solve_args.rng_type = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type;
pfluid_param = sdis_data_get(sdis_medium_get_data(fluid1));
- pfluid_param->temperature = UNKOWN_TEMPERATURE;
+ pfluid_param->temperature = SDIS_TEMPERATURE_NONE;
/* Check simulation error handling */
BA(sdis_solve_camera(scn, &solve_args, &buf));
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -272,7 +272,7 @@ main(int argc, char** argv)
solid_param->lambda = 0.1;
solid_param->rho = 1.0;
solid_param->delta = 1.0/20.0;
- solid_param->temperature = -1; /* Unknown temperature */
+ solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
OK(sdis_data_ref_put(data));
@@ -284,7 +284,7 @@ main(int argc, char** argv)
solid_param->lambda = 1.0;
solid_param->rho = 1.0;
solid_param->delta = 1.0/20.0;
- solid_param->temperature = -1; /* Unknown temperature */
+ solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
OK(sdis_data_ref_put(data));
@@ -403,7 +403,7 @@ main(int argc, char** argv)
/* Check simulation error handling when paths are registered */
solve_args.nrealisations = 10;
solve_args.register_paths = SDIS_HEAT_PATH_ALL;
- fluid_param->temperature = -1;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_medium(scn, &solve_args, &estimator));
fluid_param->temperature = Tf1;
OK(sdis_solve_medium(scn, &solve_args, &estimator));
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -256,7 +256,7 @@ main(int argc, char** argv)
solid_param->lambda = 0.1;
solid_param->rho = 1.0;
solid_param->delta = 1.0/20.0;
- solid_param->temperature = -1; /* Unknown temperature */
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
OK(sdis_data_ref_put(data));
@@ -268,7 +268,7 @@ main(int argc, char** argv)
solid_param->lambda = 1.0;
solid_param->rho = 1.0;
solid_param->delta = 1.0/20.0;
- solid_param->temperature = -1; /* Unknown temperature */
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -329,7 +329,7 @@ main(int argc, char** argv)
solid_param->lambda = 0.1;
solid_param->rho = 1.0;
solid_param->delta = 1.0/20.0;
- solid_param->temperature = -1; /* Unknown temperature */
+ solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
solid_shader.calorific_capacity = solid_get_calorific_capacity;
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
@@ -397,6 +397,9 @@ main(int argc, char** argv)
solve_args.picard_order = 0;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
solve_args.picard_order = 1;
+ solve_args.diff_algo = SDIS_DIFFUSION_NONE;
+ BA(sdis_solve_probe(scn, &solve_args, &estimator));
+ solve_args.diff_algo = SDIS_DIFFUSION_DELTA_SPHERE;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
BA(sdis_estimator_get_type(estimator, NULL));
@@ -454,7 +457,7 @@ main(int argc, char** argv)
OK(sdis_estimator_ref_put(estimator));
/* The external fluid cannot have an unknown temperature */
- fluid_param->temperature = -1;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
fluid_param->temperature = 300;
@@ -568,7 +571,7 @@ main(int argc, char** argv)
solve_args.register_paths = SDIS_HEAT_PATH_ALL;
/* Check simulation error handling when paths are registered */
- fluid_param->temperature = -1;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
fluid_param->temperature = 300;
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -78,7 +78,7 @@ temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
CHK(vtx != NULL && IS_INF(vtx->time));
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -255,7 +255,7 @@ main(int argc, char** argv)
solve_args.position[2] = 0.5;
solve_args.time_range[0] = INF;
solve_args.time_range[1] = INF;
-
+ solve_args.diff_algo = SDIS_DIFFUSION_WOS;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
ref = 350 * solve_args.position[2] + (1-solve_args.position[2]) * 300;
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -72,7 +72,7 @@ temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
CHK(vtx != NULL);
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -246,6 +246,7 @@ main(int argc, char** argv)
solve_args.position[1] = 0.5;
solve_args.time_range[0] = INF;
solve_args.time_range[1] = INF;
+ solve_args.diff_algo = SDIS_DIFFUSION_WOS;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -97,7 +97,7 @@ temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
CHK(vtx != NULL);
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -94,7 +94,7 @@ temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)data;
CHK(vtx != NULL);
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -116,7 +116,7 @@ solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
(void)vtx, (void)data;
- return -1;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -234,7 +234,7 @@ main(int argc, char** argv)
OK(sdis_green_function_ref_put(green));
/* The external fluid cannot have an unknown temperature */
- fluid_param->temperature = -1;
+ fluid_param->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
diff --git a/src/test_sdis_solve_probe_boundary_list.c b/src/test_sdis_solve_probe_boundary_list.c
@@ -157,7 +157,7 @@ mesh_add_sphere(struct mesh* mesh)
SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
-SOLID_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
SOLID_PROP(delta, 1.0/20.0) /* [m] */
static struct sdis_medium*
diff --git a/src/test_sdis_solve_probe_list.c b/src/test_sdis_solve_probe_list.c
@@ -310,7 +310,7 @@ create_scene
SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
-SOLID_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE) /* [K] */
#undef SOLID_PROP
static double
diff --git a/src/test_sdis_transcient.c b/src/test_sdis_transcient.c
@@ -18,8 +18,6 @@
#include <string.h>
-#define UNKNOWN_TEMPERATURE -1
-
/*
* The scene is composed of a solid cuboid whose temperature is fixed on its 6
* faces. The test consist in checking that the estimated temperature at a
@@ -225,7 +223,7 @@ solid_get_temperature
if(vtx->time <= 0) {
return ((const struct solid*)sdis_data_cget(data))->init_temperature;
} else {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
}
@@ -542,7 +540,7 @@ main(int argc, char** argv)
interfs[3] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[3]);
interfs[4] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[4]);
interfs[5] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[5]);
- interfs[6] = create_interface(dev, solid, solid, &interf_shader, UNKNOWN_TEMPERATURE);
+ interfs[6] = create_interface(dev, solid, solid, &interf_shader, SDIS_TEMPERATURE_NONE);
/* Setup the box scene context */
ctx.indices = box_indices;
diff --git a/src/test_sdis_unstationary_atm.c b/src/test_sdis_unstationary_atm.c
@@ -26,7 +26,7 @@
* The physical configuration is the following: a slab of fluid with known
* thermophysical properties but unknown temperature is located between a
* "ground" and a slab of solid, with also a unknown temperature profile. On
- * the other side of the solid slab, is a "atmosphere" with known temperature,
+ * the other side of the solid slab, is an "atmosphere" with known temperature,
* and known radiative temperature.
*
* Solving the system means: finding the temperature of the ground, of the
@@ -35,7 +35,7 @@
* reference)
*
* The reference for this system comes from a numerical method and is not
- * analytic. Thus the compliance test MC VS reference is not the usual |MC -
+ * analytic. Thus the compliance test MC VS reference is not the usual |MC -
* ref| <= 3*sigma but is |MC -ref| <= (Tmax -Tmin) * 0.01.
*
* 3D 2D
@@ -60,8 +60,6 @@
#define T0_SOLID 300
#define T0_FLUID 300
-#define UNKNOWN_TEMPERATURE -1
-
#define N 10000 /* #realisations */
#define TG 310
@@ -275,7 +273,7 @@ solid_get_temperature
solid = ((struct solid*)sdis_data_cget(data));
if(vtx->time <= solid->t0)
return solid->temperature;
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
struct fluid {
@@ -294,7 +292,7 @@ fluid_get_temperature
fluid = ((struct fluid*)sdis_data_cget(data));
if(vtx->time <= fluid->t0)
return fluid->temperature;
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -756,7 +754,7 @@ main(int argc, char** argv)
solid_props->rho = 1;
solid_props->delta = 1;
solid_props->t0 = INF;
- solid_props->temperature = UNKNOWN_TEMPERATURE;
+ solid_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &dummy_solid));
OK(sdis_data_ref_put(data));
@@ -786,7 +784,7 @@ main(int argc, char** argv)
OK(sdis_data_ref_put(data));
/* Create the adiabatic interfaces */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.h = 0;
interf_props.emissivity = 0;
interf_props.Tref = TREF;
@@ -794,7 +792,7 @@ main(int argc, char** argv)
create_interface(dev, solid, dummy_solid, &interf_props, &interf_adiabatic_2);
/* Create the P interface */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.h = HC;
interf_props.emissivity = 1;
interf_props.Tref = TREF;
@@ -808,7 +806,7 @@ main(int argc, char** argv)
create_interface(dev, fluid, dummy_solid, &interf_props, &interf_TG);
/* Create the TA interface */
- interf_props.temperature = UNKNOWN_TEMPERATURE;
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
interf_props.h = HA;
interf_props.emissivity = 1;
interf_props.Tref = TREF;
diff --git a/src/test_sdis_unsteady.c b/src/test_sdis_unsteady.c
@@ -0,0 +1,275 @@
+/* Copyright (C) 2016-2023 |Méso|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"
+
+/*
+ * The scene is a solid cube whose temperature is fixed on its faces. The test
+ * calculates its temperature at a given position at different observation times
+ * and validates the result against the calculated values by analytically
+ * solving the green functions.
+ *
+ * T3 (0.1,0.1,0.1)
+ * +-------+
+ * /' T4 /|
+ * +-------+ | T1
+ * T0 | +.....|.+
+ * |, T5 |/
+ * +-------+
+ * (0,0,0) T2
+ */
+
+#define T_INIT 280 /* [K] */
+#define T0 310 /* [K] */
+#define T1 320 /* [K] */
+#define T2 330 /* [K] */
+#define T3 310 /* [K] */
+#define T4 320 /* [K] */
+#define T5 300 /* [K] */
+
+#define NREALISATIONS 10000
+#define FP_TO_METER 0.1
+
+struct reference {
+ double pos[3]; /* [m/FP_TO_METER] */
+ double time; /* [s] */
+ double temp; /* [K] */
+};
+
+static const struct reference references[] = {
+ {{0.3, 0.4, 0.6}, 1000.0000, 281.33455593977152},
+ {{0.3, 0.4, 0.6}, 2000.0000, 286.90151817350699},
+ {{0.3, 0.4, 0.6}, 3000.0000, 292.84330866161531},
+ {{0.3, 0.4, 0.6}, 4000.0000, 297.81444160746452},
+ {{0.3, 0.4, 0.6}, 5000.0000, 301.70787295764546},
+ {{0.3, 0.4, 0.6}, 10000.000, 310.78920179442139},
+ {{0.3, 0.4, 0.6}, 20000.000, 313.37629443163121},
+ {{0.3, 0.4, 0.6}, 30000.000, 313.51064004438581},
+ {{0.3, 0.4, 0.6}, 1000000.0, 313.51797642855502}
+};
+static const size_t nreferences = sizeof(references)/sizeof(*references);
+
+/*******************************************************************************
+ * Solid, i.e. medium of the cube
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 2000.0) /* [J/K/kg] */
+SOLID_PROP(thermal_conductivity, 0.5) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 2500.0) /* [kg/m^3] */
+SOLID_PROP(delta, 1.0/60.0)
+#undef SOLID_PROP
+
+static double /* [K] */
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ (void)data;
+ ASSERT(vtx);
+ if(vtx->time <= 0) return T_INIT; /* Initial temperature [K] */
+ return SDIS_TEMPERATURE_NONE;
+}
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the cube. It is defined only
+ * for Stardis compliance: in Stardis, an interface must divide 2 media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is known
+ ******************************************************************************/
+static double /* [K] */
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data;
+ ASSERT(frag);
+
+ if(frag->Ng[0] == 1) return T0;
+ else if(frag->Ng[0] == -1) return T1;
+ else if(frag->Ng[1] == 1) return T2;
+ else if(frag->Ng[1] == -1) return T3;
+ else if(frag->Ng[2] == 1) return T4;
+ else if(frag->Ng[2] == -1) return T5;
+ else FATAL("Unreachable code\n");
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * The scene
+ ******************************************************************************/
+static void
+get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ (void)itri;
+ ASSERT(interf && ctx);
+ *interf = ctx;
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_scene* scn = NULL;
+
+ scn_args.get_indices = box_get_indices;
+ scn_args.get_interface = get_interface;
+ scn_args.get_position = box_get_position;
+ scn_args.nprimitives = box_ntriangles;
+ scn_args.nvertices = box_nvertices;
+ scn_args.context = interf;
+ scn_args.fp_to_meter = FP_TO_METER;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe
+ (struct sdis_scene* scn,
+ const struct reference* ref,
+ const enum sdis_diffusion_algorithm diff_algo)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ CHK(ref);
+
+ args.position[0] = ref->pos[0];
+ args.position[1] = ref->pos[1];
+ args.position[2] = ref->pos[2];
+ args.time_range[0] =
+ args.time_range[1] = ref->time; /* [s] */
+ args.diff_algo = diff_algo;
+ args.nrealisations = NREALISATIONS;
+ OK(sdis_solve_probe(scn, &args, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ printf("T(%g, %g, %g) at %g s = %g ~ %g +/- %g\n",
+ ref->pos[0]*FP_TO_METER, ref->pos[1]*FP_TO_METER, ref->pos[2]*FP_TO_METER,
+ ref->time, ref->temp, T.E, T.SE);
+ CHK(eq_eps(ref->temp, T.E, 3*T.SE));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+static void
+check
+ (struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo)
+{
+ const char* str_algo = NULL;
+ size_t i = 0;
+
+ switch(diff_algo) {
+ case SDIS_DIFFUSION_WOS: str_algo = "Walk on Sphere"; break;
+ case SDIS_DIFFUSION_DELTA_SPHERE: str_algo = "Delta sphere"; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ printf("\n\x1b[1m\x1b[35m%s\x1b[0m\n", str_algo);
+ FOR_EACH(i, 0, nreferences) {
+ check_probe(scn, references + i, diff_algo);
+ }
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL;
+ struct sdis_scene* scene = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scene = create_scene(sdis, interf);
+
+ check(scene, SDIS_DIFFUSION_DELTA_SPHERE);
+ check(scene, SDIS_DIFFUSION_WOS);
+
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scene));
+ OK(sdis_device_ref_put(sdis));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_unsteady_1d.c b/src/test_sdis_unsteady_1d.c
@@ -0,0 +1,268 @@
+/* Copyright (C) 2016-2023 |Méso|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"
+
+/*
+ * The scene is a solid 1d slab simulated by square whose temperature is fixed
+ * on two of its faces, the two others are adiabatic. The test calculates its
+ * temperature at a given position at different observation times and validates
+ * the result against the calculated values by analytically solving the green
+ * functions.
+ *
+ * ///////(0.1,0.1)
+ * +-------+
+ * | | T1
+ * | |
+ * T0 | |
+ * +-------+
+ * (0,0)///////
+ */
+
+#define T_INIT 280 /* [K] */
+#define T0 310 /* [K] */
+#define T1 320 /* [K] */
+
+#define NREALISATIONS 10000
+#define FP_TO_METER 0.1
+
+struct reference {
+ double pos[2]; /* [m/FP_TO_METER] */
+ double time; /* [s] */
+ double temp; /* [K] */
+};
+
+static const struct reference references[] = {
+ {{0.5, 0.5}, 1000.0000, 280.02848664122115},
+ {{0.5, 0.5}, 2000.0000, 280.86935314560424},
+ {{0.5, 0.5}, 3000.0000, 282.88587826961236},
+ {{0.5, 0.5}, 4000.0000, 285.39698306113996},
+ {{0.5, 0.5}, 5000.0000, 287.96909375994932},
+ {{0.5, 0.5}, 10000.000, 298.39293888670881},
+ {{0.5, 0.5}, 20000.000, 308.80965010883347},
+ {{0.5, 0.5}, 30000.000, 312.69280796373141},
+ {{0.5, 0.5}, 1000000.0, 315.00000000000000}
+};
+static const size_t nreferences = sizeof(references)/sizeof(*references);
+
+/*******************************************************************************
+ * Solid, i.e. medium of the cube
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 2000.0) /* [J/K/kg] */
+SOLID_PROP(thermal_conductivity, 0.5) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 2500.0) /* [kg/m^3] */
+SOLID_PROP(delta, 1.0/60.0)
+#undef SOLID_PROP
+
+static double /* [K] */
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ (void)data;
+ ASSERT(vtx);
+ if(vtx->time <= 0) return T_INIT; /* Initial temperature [K] */
+ return SDIS_TEMPERATURE_NONE; /* Unknown temperature */
+}
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the cube. It is defined only
+ * for Stardis compliance: in Stardis, an interface must divide 2 media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface of the system
+ ******************************************************************************/
+static double /* [K] */
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data;
+ ASSERT(frag);
+
+ if(frag->Ng[0] == 1) return T0;
+ else if(frag->Ng[0] == -1) return T1;
+ else if(frag->Ng[1] == 1) return SDIS_TEMPERATURE_NONE;
+ else if(frag->Ng[1] == -1) return SDIS_TEMPERATURE_NONE;
+ else FATAL("Unreachable code\n");
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * The scene
+ ******************************************************************************/
+static void
+get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ (void)itri;
+ ASSERT(interf && ctx);
+ *interf = ctx;
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_scene* scn = NULL;
+
+ scn_args.get_indices = square_get_indices;
+ scn_args.get_interface = get_interface;
+ scn_args.get_position = square_get_position;
+ scn_args.nprimitives = square_nsegments;
+ scn_args.nvertices = square_nvertices;
+ scn_args.context = interf;
+ scn_args.fp_to_meter = FP_TO_METER;
+ OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
+
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe
+ (struct sdis_scene* scn,
+ const struct reference* ref,
+ const enum sdis_diffusion_algorithm diff_algo)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ CHK(ref);
+
+ args.position[0] = ref->pos[0];
+ args.position[1] = ref->pos[1];
+ args.time_range[0] =
+ args.time_range[1] = ref->time; /* [s] */
+ args.diff_algo = diff_algo;
+ args.nrealisations = NREALISATIONS;
+ OK(sdis_solve_probe(scn, &args, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ printf("T(%g, %g) at %g s = %g ~ %g +/- %g\n",
+ ref->pos[0]*FP_TO_METER, ref->pos[1]*FP_TO_METER,
+ ref->time, ref->temp, T.E, T.SE);
+ CHK(eq_eps(ref->temp, T.E, 3*T.SE));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+static void
+check
+ (struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo)
+{
+ const char* str_algo = NULL;
+ size_t i = 0;
+
+ switch(diff_algo) {
+ case SDIS_DIFFUSION_WOS: str_algo = "Walk on Sphere"; break;
+ case SDIS_DIFFUSION_DELTA_SPHERE: str_algo = "Delta sphere"; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ printf("\n\x1b[1m\x1b[35m%s\x1b[0m\n", str_algo);
+ FOR_EACH(i, 0, nreferences) {
+ check_probe(scn, references + i, diff_algo);
+ }
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL;
+ struct sdis_scene* scene = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scene = create_scene(sdis, interf);
+
+ check(scene, SDIS_DIFFUSION_DELTA_SPHERE);
+ check(scene, SDIS_DIFFUSION_WOS);
+
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scene));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_unsteady_analytic_profile.c b/src/test_sdis_unsteady_analytic_profile.c
@@ -0,0 +1,380 @@
+/* Copyright (C) 2016-2023 |Méso|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/mem_allocator.h>
+
+/*
+ * The system is an unsteady-state temperature profile, meaning that at any
+ * point, at any: time, we can analytically calculate the temperature. We
+ * immerse in this temperature field a supershape representing a solid in which
+ * we want to evaluate the temperature by Monte Carlo at a given position and
+ * observation time. On the Monte Carlo side, the temperature of the supershape
+ * is unknown. Only the boundary temperature is fixed at the temperature of the
+ * unsteady trilinear profile mentioned above. Monte Carlo would then have to
+ * find the temperature defined by the unsteady profile.
+ *
+ * T(z) /\ <-- T(x,y,z,t)
+ * | T(y) ___/ \___
+ * |/ \ . T=? /
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ */
+
+#define LAMBDA 0.1 /* [W/(m.K)] */
+#define RHO 25.0 /* [kg/m^3] */
+#define CP 2.0 /* [J/K/kg)] */
+#define DELTA 1.0/20.0 /* [m/fp_to_meter] */
+
+#define NREALISATIONS 100000
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+temperature(const double pos[3], const double time)
+{
+ const double kx = PI/4.0;
+ const double ky = PI/4.0;
+ const double kz = PI/4.0;
+ const double alpha = LAMBDA / (RHO*CP); /* Diffusivity */
+
+ const double A = 0; /* No termal source */
+ const double B1 = 10;
+ const double B2 = 1000;
+
+ double x, y, z, t;
+ double a, b, c;
+ double temp;
+
+ ASSERT(pos);
+
+ x = pos[0];
+ y = pos[1];
+ z = pos[2];
+ t = time;
+
+ a = B1*(x*x*x*z-3*x*y*y*z);
+ b = B2*sin(kx*x)*sin(ky*y)*sin(kz*z)*exp(-alpha*(kx*kx + ky*ky + kz*kz)*t);
+ c = A * x*x*x*x * y*y*y * z*z;
+
+ temp = (a + b - c) / LAMBDA;
+ return temp;
+}
+
+static INLINE void
+dump_s3dut_mesh(FILE* fp, const struct s3dut_mesh* mesh)
+{
+ struct s3dut_mesh_data mesh_data;
+
+ OK(s3dut_mesh_get_data(mesh, &mesh_data));
+ dump_mesh(fp, mesh_data .positions, mesh_data.nvertices,
+ mesh_data.indices, mesh_data.nprimitives);
+}
+
+static void
+dump_paths
+ (FILE* fp,
+ struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo,
+ const double pos[3],
+ const double time,
+ const size_t npaths)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_estimator* estimator = NULL;
+
+ args.nrealisations = npaths;
+ args.position[0] = pos[0];
+ args.position[1] = pos[1];
+ args.position[2] = pos[2];
+ args.time_range[0] = time;
+ args.time_range[1] = time;
+ args.diff_algo = diff_algo;
+ args.register_paths = SDIS_HEAT_PATH_ALL;
+ OK(sdis_solve_probe(scn, &args, &estimator));
+
+ dump_heat_paths(fp, estimator);
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static struct s3dut_mesh*
+create_super_shape(void)
+{
+ struct s3dut_mesh* mesh = NULL;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 1;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &mesh));
+
+ return mesh;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ struct s3dut_mesh_data mesh_data;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {{0}, NULL};
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < context->mesh_data.nprimitives);
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh_data.indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh_data.indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh_data.indices[itri*3+1];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < context->mesh_data.nprimitives);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < context->mesh_data.nvertices);
+ pos[0] = context->mesh_data.positions[ivert*3+0];
+ pos[1] = context->mesh_data.positions[ivert*3+1];
+ pos[2] = context->mesh_data.positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct s3dut_mesh* mesh,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ OK(s3dut_mesh_get_data(mesh, &context.mesh_data));
+ context.interf = interf;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = context.mesh_data.nprimitives;
+ scn_args.nvertices = context.mesh_data.nvertices;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of the super shape
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, CP) /* [J/K/kg] */
+SOLID_PROP(thermal_conductivity, LAMBDA) /* [W/m/K] */
+SOLID_PROP(volumic_mass, RHO) /* [kg/m^3] */
+SOLID_PROP(delta, DELTA) /* [m/fp_to_meter] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
+#undef SOLID_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.t0 = -INF;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is fixed to the unsteady-state profile
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data;
+ return temperature(frag->P, frag->time);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe
+ (struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo,
+ const double pos[3],
+ const double time, /* [s] */
+ const int green)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ args.nrealisations = NREALISATIONS;
+ args.position[0] = pos[0];
+ args.position[1] = pos[1];
+ args.position[2] = pos[2];
+ args.time_range[0] = time;
+ args.time_range[1] = time;
+ args.diff_algo = diff_algo;
+
+ if(!green) {
+ OK(sdis_solve_probe(scn, &args, &estimator));
+ } else {
+ struct sdis_green_function* greenfn = NULL;
+
+ OK(sdis_solve_probe_green_function(scn, &args, &greenfn));
+ OK(sdis_green_function_solve(greenfn, &estimator));
+ OK(sdis_green_function_ref_put(greenfn));
+ }
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = temperature(pos, time);
+ printf("T(%g, %g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), time, ref, T.E, T.SE);
+ CHK(eq_eps(ref, T.E, 3*T.SE));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ FILE* fp = NULL;
+ struct s3dut_mesh* super_shape = NULL;
+ const double pos[3] = {0.2,0.3,0.4}; /* [m/fp_to_meter] */
+ const double time = 5; /* [s] */
+
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ super_shape = create_super_shape();
+
+ /* Save the super shape geometry for debug and visualisation */
+ CHK(fp = fopen("super_shape_3d.obj", "w"));
+ dump_s3dut_mesh(fp, super_shape);
+ CHK(fclose(fp) == 0);
+
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scn = create_scene(sdis, super_shape, interf);
+
+ check_probe(scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 0/*green*/);
+ check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 0/*green*/);
+ check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 1/*green*/);
+
+ /* Write 10 heat paths sampled by the delta sphere algorithm */
+ CHK(fp = fopen("paths_delta_sphere_3d.vtk", "w"));
+ dump_paths(fp, scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 10);
+ CHK(fclose(fp) == 0);
+
+ /* Write 10 heat paths sampled by the WoS algorithm */
+ CHK(fp = fopen("paths_wos_3d.vtk", "w"));
+ dump_paths(fp, scn, SDIS_DIFFUSION_WOS, pos, time, 10);
+ CHK(fclose(fp) == 0);
+
+ OK(s3dut_mesh_ref_put(super_shape));
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scn));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_unsteady_analytic_profile_2d.c b/src/test_sdis_unsteady_analytic_profile_2d.c
@@ -0,0 +1,410 @@
+/* Copyright (C) 2016-2023 |Méso|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/mem_allocator.h>
+#include <rsys/stretchy_array.h>
+
+/*
+ * The system is an unsteady-state temperature profile, meaning that at any
+ * point, at any: time, we can analytically calculate the temperature. We
+ * immerse in this temperature field a supershape representing a solid in which
+ * we want to evaluate the temperature by Monte Carlo at a given position and
+ * observation time. On the Monte Carlo side, the temperature of the supershape
+ * is unknown. Only the boundary temperature is fixed at the temperature of the
+ * unsteady trilinear profile mentioned above. Monte Carlo would then have to
+ * find the temperature defined by the unsteady profile.
+ *
+ * T(y) /\ <-- T(x,y,t)
+ * | ___/ \___
+ * | \ . T=? /
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ *
+ */
+
+#define LAMBDA 0.1 /* [W/(m.K)] */
+#define RHO 25.0 /* [kg/m^3] */
+#define CP 2.0 /* [J/K/kg)] */
+#define DELTA 1.0/20.0 /* [m/fp_to_meter] */
+
+#define NREALISATIONS 100000
+
+struct super_shape {
+ double* positions;
+ size_t* indices;
+};
+#define SUPER_SHAPE_NULL__ {NULL, NULL}
+static const struct super_shape SUPER_SHAPE_NULL = SUPER_SHAPE_NULL__;
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* Unsteady-state temperature profile. Its corresponding power term is :
+ * P(x, y) = A*[12*x^2*y^3 + 5*x^4*y] */
+static double temperature(const double pos[2], const double time)
+{
+ const double kx = PI/4;
+ const double ky = PI/4;
+ const double alpha = LAMBDA / (RHO*CP); /* Diffusivity */
+
+ const double A = 0; /* No termal source */
+ const double B1 = 0;
+ const double B2 = 1000;
+
+ double x, y, t;
+ double a, b, c;
+ double temp;
+ ASSERT(pos);
+
+ x = pos[0];
+ y = pos[1];
+ t = time;
+
+ a = B1*(x*x*x - 3*x*y*y);
+ b = B2*sin(kx*x)*sin(ky*y)*exp(-alpha*(kx*kx + ky*ky)*t);
+ c = A * x*x*x*x + y*y*y;
+
+ temp = (a + b - c) / LAMBDA;
+ return temp;
+}
+
+static void
+dump_paths
+ (FILE* fp,
+ struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo,
+ const double pos[2],
+ const double time,
+ const size_t npaths)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_estimator* estimator = NULL;
+
+ args.nrealisations = npaths;
+ args.position[0] = pos[0];
+ args.position[1] = pos[1];
+ args.time_range[0] = time;
+ args.time_range[1] = time;
+ args.diff_algo = diff_algo;
+ args.register_paths = SDIS_HEAT_PATH_ALL;
+ OK(sdis_solve_probe(scn, &args, &estimator));
+
+ dump_heat_paths(fp, estimator);
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of the super shape
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, CP) /* [J/K/kg] */
+SOLID_PROP(thermal_conductivity, LAMBDA) /* [W/m/K] */
+SOLID_PROP(volumic_mass, RHO) /* [kg/m^3] */
+SOLID_PROP(delta, DELTA) /* [m/fp_to_meter] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
+#undef SOLID_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.t0 = -INF;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is fixed to the unsteady-state profile
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data;
+ return temperature(frag->P, frag->time);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Super shape
+ ******************************************************************************/
+static struct super_shape
+create_super_shape(void)
+{
+ struct super_shape sshape = SUPER_SHAPE_NULL;
+
+ const unsigned nslices = 128;
+ const double a = 1.0;
+ const double b = 1.0;
+ const double n1 = 1.0;
+ const double n2 = 1.0;
+ const double n3 = 1.0;
+ const double m = 6.0;
+ size_t i = 0;
+
+ FOR_EACH(i, 0, nslices) {
+ const double theta = (double)i * (2.0*PI / (double)nslices);
+ const double tmp0 = pow(fabs(1.0/a * cos(m/4.0*theta)), n2);
+ const double tmp1 = pow(fabs(1.0/b * sin(m/4.0*theta)), n3);
+ const double tmp2 = pow(tmp0 + tmp1, 1.0/n1);
+ const double r = 1.0 / tmp2;
+ const double x = cos(theta) * r;
+ const double y = sin(theta) * r;
+ const size_t j = (i + 1) % nslices;
+
+ sa_push(sshape.positions, x);
+ sa_push(sshape.positions, y);
+ sa_push(sshape.indices, i);
+ sa_push(sshape.indices, j);
+ }
+
+ return sshape;
+}
+
+static void
+release_super_shape(struct super_shape* sshape)
+{
+ CHK(sshape);
+ sa_release(sshape->positions);
+ sa_release(sshape->indices);
+}
+
+static INLINE size_t
+super_shape_nsegments(const struct super_shape* sshape)
+{
+ CHK(sshape);
+ return sa_size(sshape->indices) / 2/*#indices per segment*/;
+}
+
+static INLINE size_t
+super_shape_nvertices(const struct super_shape* sshape)
+{
+ CHK(sshape);
+ return sa_size(sshape->positions) / 2/*#coords per vertex*/;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct super_shape* sshape;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL};
+
+static void
+scene_get_indices(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context);
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->sshape->indices[iseg*2+1];
+ ids[1] = (unsigned)context->sshape->indices[iseg*2+0];
+}
+
+static void
+scene_get_interface(const size_t iseg, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ (void)iseg;
+ CHK(interf && context);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context);
+ pos[0] = context->sshape->positions[ivert*2+0];
+ pos[1] = context->sshape->positions[ivert*2+1];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct super_shape* sshape,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.interf = interf;
+ context.sshape = sshape;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = super_shape_nsegments(sshape);
+ scn_args.nvertices = super_shape_nvertices(sshape);
+ scn_args.context = &context;
+ OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe
+ (struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo,
+ const double pos[2],
+ const double time, /* [s] */
+ const int green)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ args.nrealisations = NREALISATIONS;
+ args.position[0] = pos[0];
+ args.position[1] = pos[1];
+ args.time_range[0] = time;
+ args.time_range[1] = time;
+ args.diff_algo = diff_algo;
+
+ if(!green) {
+ OK(sdis_solve_probe(scn, &args, &estimator));
+ } else {
+ struct sdis_green_function* greenfn = NULL;
+
+ OK(sdis_solve_probe_green_function(scn, &args, &greenfn));
+ OK(sdis_green_function_solve(greenfn, &estimator));
+ OK(sdis_green_function_ref_put(greenfn));
+ }
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = temperature(pos, time);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), time, ref, T.E, T.SE);
+ CHK(eq_eps(ref, T.E, 3*T.SE));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ FILE* fp = NULL;
+ struct super_shape sshape = SUPER_SHAPE_NULL;
+ const double pos[3] = {0.2,0.3}; /* [m/fp_to_meter] */
+ const double time = 5; /* [s] */
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ sshape = create_super_shape();
+
+ /* Save the super shape geometry for debug and visualisation */
+ CHK(fp = fopen("super_shape_2d.obj", "w"));
+ dump_segments(fp, sshape.positions, super_shape_nvertices(&sshape),
+ sshape.indices, super_shape_nsegments(&sshape));
+ CHK(fclose(fp) == 0);
+
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scn = create_scene(sdis, &sshape, interf);
+
+ check_probe(scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 0/*green*/);
+ check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 0/*green*/);
+ check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 1/*green*/);
+
+ /* Write 10 heat paths sampled by the delta sphere algorithm */
+ CHK(fp = fopen("paths_delta_sphere_2d.vtk", "w"));
+ dump_paths(fp, scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 10);
+ CHK(fclose(fp) == 0);
+
+ /* Write 10 heat paths sampled by the WoS algorithm */
+ CHK(fp = fopen("paths_wos_2d.vtk", "w"));
+ dump_paths(fp, scn, SDIS_DIFFUSION_WOS, pos, time, 10);
+ CHK(fclose(fp) == 0);
+
+ release_super_shape(&sshape);
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scn));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -428,4 +428,3 @@ check_green_serialization
struct sdis_scene* scn);
#endif /* TEST_SDIS_UTILS_H */
-
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -45,7 +45,6 @@
* (0,0,0) /////
*/
-#define UNKNOWN_TEMPERATURE -1
#define N 10000 /* #realisations */
#define T0 320
@@ -54,6 +53,22 @@
#define DELTA 1.0/55.0
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static const char*
+algo_cstr(const enum sdis_diffusion_algorithm diff_algo)
+{
+ const char* cstr = "none";
+
+ switch(diff_algo) {
+ case SDIS_DIFFUSION_DELTA_SPHERE: cstr = "delta sphere"; break;
+ case SDIS_DIFFUSION_WOS: cstr = "WoS"; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return cstr;
+}
+
+/*******************************************************************************
* Media
******************************************************************************/
struct solid {
@@ -72,7 +87,7 @@ fluid_get_temperature
{
(void)data;
CHK(vtx != NULL);
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
}
static double
@@ -116,7 +131,7 @@ solid_get_temperature
CHK(data != NULL);
t0 = ((const struct solid*)sdis_data_cget(data))->t0;
if(vtx->time > t0) {
- return UNKNOWN_TEMPERATURE;
+ return SDIS_TEMPERATURE_NONE;
} else {
return ((const struct solid*)sdis_data_cget(data))->initial_temperature;
}
@@ -174,7 +189,6 @@ solve
struct sdis_mc time = SDIS_MC_NULL;
size_t nreals;
size_t nfails;
- double ref = -1;
enum sdis_scene_dimension dim;
const int nsimuls = 4;
int isimul;
@@ -182,11 +196,16 @@ solve
OK(sdis_scene_get_dimension(scn, &dim));
FOR_EACH(isimul, 0, nsimuls) {
- int steady = (isimul % 2) == 0;
+ const enum sdis_diffusion_algorithm algo = (isimul / 2) % 2
+ ? SDIS_DIFFUSION_WOS
+ : SDIS_DIFFUSION_DELTA_SPHERE;
+ const int steady = (isimul % 2) == 0;
+ double power = P0 == SDIS_VOLUMIC_POWER_NONE ? 0 : P0;
/* Restore power value */
solid->vpower = P0;
+ solve_args.diff_algo = algo;
solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9);
solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9);
solve_args.position[2] =
@@ -210,42 +229,28 @@ solve
OK(sdis_estimator_get_temperature(estimator, &T));
OK(sdis_estimator_get_realisation_time(estimator, &time));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s "
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ } else {
+ printf(
+ "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
printf("Elapsed time = %s\n", dump);
printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- if(steady) CHK(eq_eps(T.E, ref, T.SE * 3));
/* Check green function */
time_current(&t0);
@@ -258,37 +263,21 @@ solve
OK(sdis_estimator_get_failure_count(estimator2, &nfails));
OK(sdis_estimator_get_temperature(estimator2, &T));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Green Steady temperature at (%g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Green Mean temperature at (%g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Green Steady temperature at (%g, %g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Green Mean temperature at (%g, %g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Green steady temperature - pos: %g, %g, %g; Power: %g algo: %s"
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ } else {
+ printf(
+ "Green mean temperature - pos: %g, %g, %g; t: [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
@@ -306,7 +295,7 @@ solve
printf("\n");
/* Check same green used at a different power level */
- solid->vpower = 3 * P0;
+ solid->vpower = power = 3 * P0;
time_current(&t0);
OK(sdis_solve_probe(scn, &solve_args, &estimator));
@@ -318,43 +307,28 @@ solve
OK(sdis_estimator_get_temperature(estimator, &T));
OK(sdis_estimator_get_realisation_time(estimator, &time));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s "
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ } else {
+ printf(
+ "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
printf("Elapsed time = %s\n", dump);
printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- if(steady) CHK(eq_eps(T.E, ref, T.SE * 3));
time_current(&t0);
OK(sdis_green_function_solve(green, &estimator2));
@@ -466,7 +440,7 @@ main(int argc, char** argv)
/* Create the adiabatic interface */
OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
interf_props = sdis_data_get(data);
- interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create
(dev, solid, fluid, &interf_shader, data, &interf_adiabatic));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c
@@ -19,7 +19,6 @@
#define N 10000 /* #realisations */
#define Pw 10000 /* Volumic power */
-#define NONE -1
#define DELTA 0.01
#define DELTA_PSQUARE 0.01
@@ -294,10 +293,10 @@ main(int argc, char** argv)
{{0,-0.55, 0}, 90.250, 90.040}
};
const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
- {{0, 0.85}, 678.170, -1},
- {{0, 0.65}, 1520.84, -1},
- {{0, 0.45}, 1794.57, -1},
- {{0, 0.25}, 1429.74, -1}
+ {{0, 0.85}, 678.170, 0},
+ {{0, 0.65}, 1520.84, 0},
+ {{0, 0.45}, 1794.57, 0},
+ {{0, 0.25}, 1429.74, 0}
};
(void)argc, (void)argv;
@@ -341,7 +340,7 @@ main(int argc, char** argv)
solid_param->lambda = 1;
solid_param->delta = DELTA;
solid_param->P = SDIS_VOLUMIC_POWER_NONE;
- solid_param->T = -1;
+ solid_param->T = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
OK(sdis_data_ref_put(data));
@@ -354,7 +353,7 @@ main(int argc, char** argv)
solid_param->lambda = 10;
solid_param->delta = DELTA_PSQUARE;
solid_param->P = Pw;
- solid_param->T = -1;
+ solid_param->T = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
OK(sdis_data_ref_put(data));
@@ -387,7 +386,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = 5;
- interf_param->temperature = NONE;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
&interf_solid1_fluid1));
OK(sdis_data_ref_put(data));
@@ -397,7 +396,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = 10;
- interf_param->temperature = NONE;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data,
&interf_solid1_fluid2));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c
@@ -19,11 +19,10 @@
#define N 10000 /* #realisations */
#define Pw 10000 /* Volumic power */
-#define NONE -1
/* H delta T */
-#define Tboundary1 NONE
-#define Tboundary2 NONE
+#define Tboundary1 SDIS_TEMPERATURE_NONE
+#define Tboundary2 SDIS_TEMPERATURE_NONE
#define DELTA 0.01
#define Tref 286.83 /* In Celsius. Computed with Syrthes at the position 0.5 */
@@ -373,7 +372,7 @@ main(int argc, char** argv)
solid_param->lambda = 1;
solid_param->delta = DELTA;
solid_param->P = SDIS_VOLUMIC_POWER_NONE;
- solid_param->T = -1;
+ solid_param->T = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
OK(sdis_data_ref_put(data));
@@ -386,7 +385,7 @@ main(int argc, char** argv)
solid_param->lambda = 10;
solid_param->delta = DELTA;
solid_param->P = Pw;
- solid_param->T = -1;
+ solid_param->T = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
OK(sdis_data_ref_put(data));
diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c
@@ -41,11 +41,10 @@
#define Tb 1207.1122
/* Fixed temperatures */
-#define UNKNOWN_TEMPERATURE -1
-#define Tsolid1_fluid UNKNOWN_TEMPERATURE /*Tp1*/
-#define Tsolid2_fluid UNKNOWN_TEMPERATURE /*Tp2*/
-#define Tsolid_solid1 UNKNOWN_TEMPERATURE /*Ta*/
-#define Tsolid_solid2 UNKNOWN_TEMPERATURE /*Tb*/
+#define Tsolid1_fluid SDIS_TEMPERATURE_NONE /*Tp1*/
+#define Tsolid2_fluid SDIS_TEMPERATURE_NONE /*Tp2*/
+#define Tsolid_solid1 SDIS_TEMPERATURE_NONE /*Ta*/
+#define Tsolid_solid2 SDIS_TEMPERATURE_NONE /*Tb*/
#define PROBE_POS 1.8
@@ -304,7 +303,7 @@ main(int argc, char** argv)
solid_param->lambda = LAMBDA1;
solid_param->delta = DELTA1;
solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
- solid_param->temperature = -1;
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
OK(sdis_data_ref_put(data));
@@ -317,7 +316,7 @@ main(int argc, char** argv)
solid_param->lambda = LAMBDA2;
solid_param->delta = DELTA2;
solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
- solid_param->temperature = -1;
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
OK(sdis_data_ref_put(data));
@@ -330,7 +329,7 @@ main(int argc, char** argv)
solid_param->lambda = LAMBDA;
solid_param->delta = DELTA;
solid_param->volumic_power = Pw;
- solid_param->temperature = -1;
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid));
OK(sdis_data_ref_put(data));
@@ -363,7 +362,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = 0;
- interf_param->temperature = -1;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data,
&interf_solid_adiabatic));
OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data,
diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c
@@ -270,7 +270,7 @@ main(int argc, char** argv)
solid_param->lambda = LAMBDA;
solid_param->delta = DELTA;
solid_param->volumic_power = Power;
- solid_param->temperature = -1;
+ solid_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_solid_create(dev, &solid_shader, data, &solid));
OK(sdis_data_ref_put(data));
@@ -283,7 +283,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = 0;
- interf_param->temperature = -1;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
&interf_adiabatic));
OK(sdis_data_ref_put(data));
@@ -293,7 +293,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = H;
- interf_param->temperature = -1;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
&interf_solid_fluid1));
OK(sdis_data_ref_put(data));
@@ -303,7 +303,7 @@ main(int argc, char** argv)
NULL, &data));
interf_param = sdis_data_get(data);
interf_param->h = H;
- interf_param->temperature = -1;
+ interf_param->temperature = SDIS_TEMPERATURE_NONE;
OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
&interf_solid_fluid2));
OK(sdis_data_ref_put(data));
