stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit 0d6e9f4fc4dfa1432fd85dcb07ef442cdb2b0460
parent e6e423a19f4bdb581259b3ca300213e300afca7c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  2 Jul 2020 09:31:51 +0200

Merge branch 'feature_rng_state' into develop

Diffstat:
Msrc/sdis.h | 388+++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------
Msrc/sdis_device.c | 57+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/sdis_device_c.h | 10++++++++++
Msrc/sdis_estimator.c | 29+++++++++++++++++++++++++++++
Msrc/sdis_estimator_buffer.c | 32+++++++++++++++++++++++++++++++-
Msrc/sdis_estimator_buffer_c.h | 5+++++
Msrc/sdis_estimator_c.h | 11+++++++++--
Msrc/sdis_green.c | 2+-
Msrc/sdis_realisation_Xd.h | 10+++++-----
Msrc/sdis_solve.c | 213++++++++++++++++++++++++++-----------------------------------------------------
Msrc/sdis_solve_boundary_Xd.h | 169++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/sdis_solve_medium_Xd.h | 67+++++++++++++++++++++++++++++++++++++++++++++----------------------
Msrc/sdis_solve_probe_Xd.h | 65++++++++++++++++++++++++++++++++++++++++++++---------------------
Msrc/sdis_solve_probe_boundary_Xd.h | 154++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------
Msrc/test_sdis_conducto_radiative.c | 29+++++++++++++++++------------
Msrc/test_sdis_conducto_radiative_2d.c | 27++++++++++++++++-----------
Msrc/test_sdis_convection.c | 35+++++++++++++++++++++--------------
Msrc/test_sdis_convection_non_uniform.c | 36++++++++++++++++++++++--------------
Msrc/test_sdis_flux.c | 16++++++++++++----
Msrc/test_sdis_solid_random_walk_robustness.c | 34+++++++++++++++++++++-------------
Msrc/test_sdis_solve_boundary.c | 247++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/test_sdis_solve_boundary_flux.c | 138++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------
Msrc/test_sdis_solve_camera.c | 94++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Msrc/test_sdis_solve_medium.c | 90++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
Msrc/test_sdis_solve_medium_2d.c | 35++++++++++++++++++++---------------
Msrc/test_sdis_solve_probe.c | 90++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------
Msrc/test_sdis_solve_probe2.c | 24+++++++++++++-----------
Msrc/test_sdis_solve_probe2_2d.c | 32+++++++++++++++++---------------
Msrc/test_sdis_solve_probe3.c | 24+++++++++++++-----------
Msrc/test_sdis_solve_probe3_2d.c | 21+++++++++++----------
Msrc/test_sdis_solve_probe_2d.c | 22++++++++++++----------
Msrc/test_sdis_volumic_power.c | 16++++++++++++----
Msrc/test_sdis_volumic_power2.c | 14+++++++++-----
Msrc/test_sdis_volumic_power2_2d.c | 12++++++++----
Msrc/test_sdis_volumic_power3_2d.c | 8+++++++-
Msrc/test_sdis_volumic_power4.c | 13++++++++++---
36 files changed, 1429 insertions(+), 840 deletions(-)

diff --git a/src/sdis.h b/src/sdis.h @@ -49,6 +49,7 @@ struct logger; struct mem_allocator; struct senc2d_scene; struct senc3d_scene; +struct ssp_rng; /* Forward declaration of the Stardis opaque data types. These data types are * ref counted. Once created the caller implicitly owns the created data, i.e. @@ -69,49 +70,20 @@ struct sdis_scene; /* Forward declaration of non ref counted types */ struct sdis_heat_path; +/******************************************************************************* + * Miscellaneous data types + ******************************************************************************/ enum sdis_side { SDIS_FRONT, SDIS_BACK, SDIS_SIDE_NULL__ }; -enum sdis_medium_type { - SDIS_FLUID, - SDIS_SOLID, - SDIS_MEDIUM_TYPES_COUNT__ -}; - -enum sdis_estimator_type { - SDIS_ESTIMATOR_TEMPERATURE, - SDIS_ESTIMATOR_FLUX, - SDIS_ESTIMATOR_TYPES_COUNT__ -}; - enum sdis_scene_dimension { SDIS_SCENE_2D, SDIS_SCENE_3D }; -enum sdis_point_type { - SDIS_FRAGMENT, - SDIS_VERTEX, - SDIS_POINT_TYPES_COUNT__, - SDIS_POINT_NONE = SDIS_POINT_TYPES_COUNT__ -}; - -enum sdis_heat_vertex_type { - SDIS_HEAT_VERTEX_CONDUCTION, - SDIS_HEAT_VERTEX_CONVECTION, - SDIS_HEAT_VERTEX_RADIATIVE -}; - -enum sdis_heat_path_flag { - SDIS_HEAT_PATH_SUCCESS = BIT(0), - SDIS_HEAT_PATH_FAILURE = BIT(1), - SDIS_HEAT_PATH_ALL = SDIS_HEAT_PATH_SUCCESS | SDIS_HEAT_PATH_FAILURE, - SDIS_HEAT_PATH_NONE = 0 -}; - /* Random walk vertex, i.e. a spatiotemporal position at a given step of the * random walk. */ struct sdis_rwalk_vertex { @@ -137,6 +109,15 @@ struct sdis_interface_fragment { static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL = SDIS_INTERFACE_FRAGMENT_NULL__; +/******************************************************************************* + * Estimation data types + ******************************************************************************/ +enum sdis_estimator_type { + SDIS_ESTIMATOR_TEMPERATURE, + SDIS_ESTIMATOR_FLUX, + SDIS_ESTIMATOR_TYPES_COUNT__ +}; + /* Monte-Carlo estimation */ struct sdis_mc { double E; /* Expected value */ @@ -146,6 +127,15 @@ struct sdis_mc { #define SDIS_MC_NULL__ {0, 0, 0} static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__; +/******************************************************************************* + * Data type used to describe physical properties + ******************************************************************************/ +enum sdis_medium_type { + SDIS_FLUID, + SDIS_SOLID, + SDIS_MEDIUM_TYPES_COUNT__ +}; + /* Functor type used to retrieve the spatio temporal physical properties of a * medium. */ typedef double @@ -236,6 +226,22 @@ struct sdis_interface_shader { static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL = SDIS_INTERFACE_SHADER_NULL__; +/******************************************************************************* + * Registered heat path data types + ******************************************************************************/ +enum sdis_heat_vertex_type { + SDIS_HEAT_VERTEX_CONDUCTION, + SDIS_HEAT_VERTEX_CONVECTION, + SDIS_HEAT_VERTEX_RADIATIVE +}; + +enum sdis_heat_path_flag { + SDIS_HEAT_PATH_SUCCESS = BIT(0), + SDIS_HEAT_PATH_FAILURE = BIT(1), + SDIS_HEAT_PATH_ALL = SDIS_HEAT_PATH_SUCCESS | SDIS_HEAT_PATH_FAILURE, + SDIS_HEAT_PATH_NONE = 0 +}; + /* Vertex of heat path v*/ struct sdis_heat_vertex { double P[3]; @@ -247,6 +253,28 @@ struct sdis_heat_vertex { static const struct sdis_heat_vertex SDIS_HEAT_VERTEX_NULL = SDIS_HEAT_VERTEX_NULL__; +/* Functor used to process a heat path registered against the estimator */ +typedef res_T +(*sdis_process_heat_path_T) + (const struct sdis_heat_path* path, + void* context); + +/* Functor used to process the vertices of a heat path */ +typedef res_T +(*sdis_process_heat_vertex_T) + (const struct sdis_heat_vertex* vertex, + void* context); + +/******************************************************************************* + * Green function data types + ******************************************************************************/ +enum sdis_point_type { + SDIS_FRAGMENT, + SDIS_VERTEX, + SDIS_POINT_TYPES_COUNT__, + SDIS_POINT_NONE = SDIS_POINT_TYPES_COUNT__ +}; + /* Path used to estimate the green function */ struct sdis_green_path { /* Internal data. Should not be accessed */ @@ -257,16 +285,17 @@ struct sdis_green_path { static const struct sdis_green_path SDIS_GREEN_PATH_NULL = SDIS_GREEN_PATH_NULL__; +/* Spatio temporal point */ struct sdis_point { union { struct { struct sdis_medium* medium; struct sdis_rwalk_vertex vertex; - } mdmvert; + } mdmvert; /* Medium and a vertex into it */ struct { struct sdis_interface* intface; struct sdis_interface_fragment fragment; - } itfrag; + } itfrag; /* Interface and a fragmetn onto it */ } data; enum sdis_point_type type; }; @@ -296,17 +325,181 @@ typedef res_T const double flux_term, void* context); -/* Functor used to process a heat path registered against the estimator */ -typedef res_T -(*sdis_process_heat_path_T) - (const struct sdis_heat_path* path, - void* context); - -/* Functor used to process the vertices of a heat path */ -typedef res_T -(*sdis_process_heat_vertex_T) - (const struct sdis_heat_vertex* vertex, - void* context); +/******************************************************************************* + * Data types of the input simulation parameters + ******************************************************************************/ +struct sdis_solve_probe_args { + size_t nrealisations; /* #realisations */ + double position[3]; /* Probe position */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_PROBE_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + {0,0,0}, /* Position */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1.0, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + SDIS_HEAT_PATH_NONE, /* Register paths mask */ \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_probe_args SDIS_SOLVE_PROBE_ARGS_DEFAULT = + SDIS_SOLVE_PROBE_ARGS_DEFAULT__; + +/* Arguments of a probe simulation */ +struct sdis_solve_probe_boundary_args { + size_t nrealisations; /* #realisations */ + size_t iprim; /* Identifier of the primitive on which the probe lies */ + double uv[2]; /* Parametric coordinates of the probe onto the primitve */ + double time_range[2]; /* Observation time */ + enum sdis_side side; /* Side of iprim on which the probe lies */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + 0, /* Primitive identifier */ \ + {0,0}, /* UV */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + SDIS_SIDE_NULL__, \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + SDIS_HEAT_PATH_NONE, \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_probe_boundary_args +SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT = + SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__; + +struct sdis_solve_boundary_args { + size_t nrealisations; /* #realisations */ + const size_t* primitives; /* List of boundary primitives to handle */ + const enum sdis_side* sides; /* Per primitive side to consider */ + size_t nprimitives; /* #primitives */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + NULL, /* List or primitive ids */ \ + NULL, /* Per primitive side */ \ + 0, /* #primitives */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + SDIS_HEAT_PATH_NONE, \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_boundary_args SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT = + SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT__; + +struct sdis_solve_medium_args { + size_t nrealisations; /* #realisations */ + struct sdis_medium* medium; /* Medium to solve */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_MEDIUM_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + NULL, /* Medium */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + SDIS_HEAT_PATH_NONE, \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_medium_args SDIS_SOLVE_MEDIUM_ARGS_DEFAULT = + SDIS_SOLVE_MEDIUM_ARGS_DEFAULT__; + +struct sdis_solve_probe_boundary_flux_args { + size_t nrealisations; /* #realisations */ + size_t iprim; /* Identifier of the primitive on which the probe lies */ + double uv[2]; /* Parametric coordinates of the probe onto the primitve */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + 0, /* Primitive identifier */ \ + {0,0}, /* UV */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_probe_boundary_flux_args +SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT = + SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT__; + +struct sdis_solve_boundary_flux_args { + size_t nrealisations; /* #realisations */ + const size_t* primitives; /* List of boundary primitives to handle */ + size_t nprimitives; /* #primitives */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ +}; +#define SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT__ { \ + 10000, /* #realisations */ \ + NULL, /* List or primitive ids */ \ + 0, /* #primitives */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + NULL /* RNG state */ \ +} +static const struct sdis_solve_boundary_flux_args +SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT = + SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT__; + +struct sdis_solve_camera_args { + struct sdis_camera* cam; /* Point of view */ + double time_range[2]; /* Observation time */ + double fp_to_meter; /* Scale from floating point units to meters */ + double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; /* In Kelvin */ + size_t image_resolution[2]; /* Image resolution */ + size_t spp; /* #samples per pixel */ + int register_paths; /* Combination of enum sdis_heat_path_flag */ +}; +#define SDIS_SOLVE_CAMERA_ARGS_DEFAULT__ { \ + NULL, /* Camera */ \ + {DBL_MAX,DBL_MAX}, /* Time range */ \ + 1, /* FP to meter */ \ + -1, /* Ambient radiative temperature */ \ + -1, /* Reference temperature */ \ + {512,512}, /* Image resolution */ \ + 256, /* #realisations per pixel */ \ + SDIS_HEAT_PATH_NONE \ +} +static const struct sdis_solve_camera_args SDIS_SOLVE_CAMERA_ARGS_DEFAULT = + SDIS_SOLVE_CAMERA_ARGS_DEFAULT__; BEGIN_DECLS @@ -394,7 +587,7 @@ sdis_camera_look_at const double up[3]); /******************************************************************************* - * An estimator buffer is 2D array of estimators + * An estimator buffer is 2D array of estimators ******************************************************************************/ SDIS_API res_T sdis_estimator_buffer_ref_get @@ -436,6 +629,11 @@ sdis_estimator_buffer_get_realisation_time (const struct sdis_estimator_buffer* buf, struct sdis_mc* time); +SDIS_API res_T +sdis_estimator_buffer_get_rng_state + (const struct sdis_estimator_buffer* buf, + struct ssp_rng** rng_state); + /******************************************************************************* * A medium encapsulates the properties of either a fluid or a solid. ******************************************************************************/ @@ -730,6 +928,14 @@ sdis_estimator_for_each_path sdis_process_heat_path_T func, void* context); +/* Retrieve the RNG state at the end of the simulation. */ +SDIS_API res_T +sdis_estimator_get_rng_state + (const struct sdis_estimator* estimator, + /* The returned value may be NULL as for instance an estimator retrieved + * from an estimator buffer */ + struct ssp_rng** rng_state); + /******************************************************************************* * The green function saves the estimation of the propagator ******************************************************************************/ @@ -838,91 +1044,43 @@ sdis_heat_path_for_each_vertex SDIS_API res_T sdis_solve_probe (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const double position[3], /* Probe position */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_args* args, struct sdis_estimator** estimator); SDIS_API res_T sdis_solve_probe_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_estimator** estimator); SDIS_API res_T sdis_solve_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const enum sdis_side sides[], /* Per primitive side to consider */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_boundary_args* args, struct sdis_estimator** estimator); SDIS_API res_T sdis_solve_probe_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_probe_boundary_flux_args* args, struct sdis_estimator** estimator); SDIS_API res_T sdis_solve_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_boundary_flux_args* args, struct sdis_estimator** estimator); SDIS_API res_T sdis_solve_camera (struct sdis_scene* scn, - const struct sdis_camera* cam, /* Point of view */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const size_t width, /* Image definition in in X */ - const size_t height, /* Image definition in Y */ - const size_t spp, /* #samples per pixel */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_camera_args* args, struct sdis_estimator_buffer** buf); SDIS_API res_T sdis_solve_medium (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - struct sdis_medium* medium, /* Medium to solve */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_medium_args* args, struct sdis_estimator** estimator); /******************************************************************************* @@ -949,45 +1107,25 @@ sdis_solve_medium SDIS_API res_T sdis_solve_probe_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const double position[3], /* Probe position */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_probe_args* args, struct sdis_green_function** green); SDIS_API res_T sdis_solve_probe_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** green); SDIS_API res_T sdis_solve_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const enum sdis_side sides[], /* Per primitive side to consider */ - const size_t nprimitives, /* #primitives */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_boundary_args* args, struct sdis_green_function** green); SDIS_API res_T sdis_solve_medium_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - struct sdis_medium* medium, /* Medium to solve */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double ambient_radiative_temperature, /* In Kelvin */ - const double reference_temperature, /* In Kelvin */ + const struct sdis_solve_medium_args* args, struct sdis_green_function** green); END_DECLS diff --git a/src/sdis_device.c b/src/sdis_device.c @@ -23,6 +23,7 @@ #include <star/s2d.h> #include <star/s3d.h> +#include <star/ssp.h> #include <omp.h> @@ -138,3 +139,59 @@ sdis_device_ref_put(struct sdis_device* dev) return RES_OK; } +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +create_rng_from_rng_proxy + (struct sdis_device* dev, + const struct ssp_rng_proxy* proxy, + struct ssp_rng** out_rng) +{ + struct ssp_rng_type rng_type; + struct ssp_rng* rng = NULL; + FILE* stream = NULL; + res_T res = RES_OK; + ASSERT(dev && proxy && out_rng); + + stream = tmpfile(); + if(!stream) { + log_err(dev, + "Could not open a temporary stream to store the RNG state.\n"); + res = RES_IO_ERR; + goto error; + } + + SSP(rng_proxy_get_type(proxy, &rng_type)); + res = ssp_rng_create(dev->allocator, &rng_type, &rng); + if(res != RES_OK) { + log_err(dev, "Could not create the RNG -- %s\n", res_to_cstr(res)); + goto error; + } + + res = ssp_rng_proxy_write(proxy, stream); + if(res != RES_OK) { + log_err(dev, "Could not serialize the RNG state -- %s\n", + res_to_cstr(res)); + goto error; + } + + rewind(stream); + res = ssp_rng_read(rng, stream); + if(res != RES_OK) { + log_err(dev, "Could not read the serialized RNG state -- %s\n", + res_to_cstr(res)); + goto error; + } + +exit: + if(out_rng) *out_rng = rng; + if(stream) fclose(stream); + return res; +error: + if(rng) { + SSP(rng_ref_put(rng)); + rng = NULL; + } + goto exit; +} diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h @@ -23,6 +23,10 @@ #include <rsys/logger.h> #include <rsys/ref_count.h> +/* Forward declarations */ +struct ssp_rng; +struct ssp_rng_proxy; + struct name { FITEM; }; #define FITEM_TYPE name #include <rsys/free_list.h> @@ -43,5 +47,11 @@ struct sdis_device { ref_T ref; }; +extern LOCAL_SYM res_T +create_rng_from_rng_proxy + (struct sdis_device* dev, + const struct ssp_rng_proxy* proxy, + struct ssp_rng** out_rng); + #endif /* SDIS_DEVICE_C_H */ diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -16,7 +16,11 @@ #include "sdis.h" #include "sdis_device_c.h" #include "sdis_estimator_c.h" +#include "sdis_log.h" +#include <star/ssp.h> + +#include <rsys/cstr.h> #include <rsys/mutex.h> /******************************************************************************* @@ -32,6 +36,7 @@ estimator_release(ref_T* ref) dev = estimator->dev; darray_heat_path_release(&estimator->paths); if(estimator->mutex) mutex_destroy(estimator->mutex); + if(estimator->rng) SSP(rng_ref_put(estimator->rng)); MEM_RM(dev->allocator, estimator); SDIS(device_ref_put(dev)); } @@ -192,6 +197,16 @@ error: goto exit; } +res_T +sdis_estimator_get_rng_state + (const struct sdis_estimator* estimator, + struct ssp_rng** rng_state) +{ + if(!estimator || !rng_state) return RES_BAD_ARG; + *rng_state = estimator->rng; + return RES_OK; +} + /******************************************************************************* * Local functions ******************************************************************************/ @@ -268,3 +283,17 @@ error: goto exit; } +res_T +estimator_save_rng_state + (struct sdis_estimator* estimator, + const struct ssp_rng_proxy* proxy) +{ + ASSERT(estimator && proxy); + /* Release the previous RNG state if any */ + if(estimator->rng) { + SSP(rng_ref_put(estimator->rng)); + estimator->rng = NULL; + } + return create_rng_from_rng_proxy(estimator->dev, proxy, &estimator->rng); +} + diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c @@ -19,6 +19,8 @@ #include "sdis_estimator_buffer_c.h" #include "sdis_log.h" +#include <star/ssp.h> + struct sdis_estimator_buffer { struct sdis_estimator** estimators; /* Row major per pixe lestimators */ size_t width; @@ -29,6 +31,9 @@ struct sdis_estimator_buffer { size_t nrealisations; /* #successes */ size_t nfailures; + /* State of the RNG after the simulation */ + struct ssp_rng* rng; + ref_T ref; struct sdis_device* dev; }; @@ -54,6 +59,7 @@ estimator_buffer_release(ref_T* ref) MEM_RM(dev->allocator, buf->estimators); } + if(buf->rng) SSP(rng_ref_put(buf->rng)); MEM_RM(dev->allocator, buf); SDIS(device_ref_put(dev)); } @@ -136,13 +142,22 @@ sdis_estimator_buffer_get_temperature res_T sdis_estimator_buffer_get_realisation_time -(const struct sdis_estimator_buffer* buf, struct sdis_mc* mc) + (const struct sdis_estimator_buffer* buf, struct sdis_mc* mc) { if(!buf || !mc) return RES_BAD_ARG; SETUP_MC(mc, &buf->realisation_time); return RES_OK; } +res_T +sdis_estimator_buffer_get_rng_state + (const struct sdis_estimator_buffer* buf, struct ssp_rng** rng_state) +{ + if(!buf || !rng_state) return RES_BAD_ARG; + *rng_state = buf->rng; + return RES_OK; +} + #undef SETUP_MC /******************************************************************************* @@ -243,3 +258,18 @@ estimator_buffer_setup_realisation_time buf->realisation_time.count = buf->nrealisations; } +res_T +estimator_buffer_save_rng_state + (struct sdis_estimator_buffer* buf, + const struct ssp_rng_proxy* proxy) +{ + ASSERT(buf && proxy); + + /* Release the previous RNG state if any */ + if(buf->rng) { + SSP(rng_ref_put(buf->rng)); + buf->rng = NULL; + } + return create_rng_from_rng_proxy(buf->dev, proxy, &buf->rng); +} + diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h @@ -54,5 +54,10 @@ estimator_buffer_setup_realisation_time const double sum, const double sum2); +extern LOCAL_SYM res_T +estimator_buffer_save_rng_state + (struct sdis_estimator_buffer* buf, + const struct ssp_rng_proxy* proxy); + #endif /* SDIS_ESTIMATOR_BUFFER_C_H */ diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -23,6 +23,7 @@ #include <rsys/ref_count.h> /* Forward declarations */ +struct ssp_rng; struct sdis_device; struct sdis_estimator; enum sdis_estimator_type; @@ -44,13 +45,14 @@ struct sdis_estimator { struct mutex* mutex; struct darray_heat_path paths; /* Tracked paths */ + /* State of the RNG after the simulation */ + struct ssp_rng* rng; + enum sdis_estimator_type type; ref_T ref; struct sdis_device* dev; }; -struct sdis_estimator_handle; - /******************************************************************************* * Estimator local API ******************************************************************************/ @@ -66,6 +68,11 @@ estimator_add_and_release_heat_path (struct sdis_estimator* estimator, struct sdis_heat_path* path); +extern LOCAL_SYM res_T +estimator_save_rng_state + (struct sdis_estimator* estimator, + const struct ssp_rng_proxy* proxy); + /* Must be invoked before any others "estimator_setup" functions */ static INLINE void estimator_setup_realisations_count diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -179,7 +179,7 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src) #define HTABLE_DATA struct sdis_interface* #include <rsys/hash_table.h> -/* Generate the hash table that maps and id to a medium */ +/* Generate the hash table that maps an id to a medium */ #define HTABLE_NAME medium #define HTABLE_KEY unsigned #define HTABLE_DATA struct sdis_medium* diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -81,10 +81,10 @@ XD(compute_temperature) * boundary. Indeed, one knows the "right" type of the first vertex only * after the boundary_path execution that defines the sub path to resolve * from the submitted boundary position. Note that if the boundary - * temperature is know, the type is let as it. */ - if(heat_vtx && !T->done) { - if(heat_path_get_last_vertex(ctx->heat_path) != heat_vtx) { - /* Path was reinjected into a solid */ + * temperature is known, the type is let as it. */ + if(heat_vtx && !T->done && T->func != XD(boundary_path)) { + heat_vtx = heat_path_get_last_vertex(ctx->heat_path); + if(T->func == XD(conductive_path)) { heat_vtx->type = SDIS_HEAT_VERTEX_CONDUCTION; } else if(T->func == XD(convective_path)) { heat_vtx->type = SDIS_HEAT_VERTEX_CONVECTION; @@ -228,7 +228,7 @@ XD(boundary_realisation) float st[2]; #endif res_T res = RES_OK; - ASSERT(uv && fp_to_meter > 0 && weight && Tref >= 0 && time >= 0); + ASSERT(uv && fp_to_meter > 0 && weight && time >= 0); T.func = XD(boundary_path); rwalk.hit_side = side; diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -79,6 +79,7 @@ solve_pixel { struct accum acc_temp = ACCUM_NULL; struct accum acc_time = ACCUM_NULL; + struct sdis_heat_path* pheat_path = NULL; size_t irealisation; res_T res = RES_OK; ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0); @@ -91,7 +92,6 @@ solve_pixel double ray_pos[3]; double ray_dir[3]; double w = 0; - struct sdis_heat_path* pheat_path = NULL; struct sdis_heat_path heat_path; double time; res_T res_simul = RES_OK; @@ -131,9 +131,11 @@ solve_pixel /* Check if the path must be saved regarding the register_paths mask */ if(!(register_paths & (int)pheat_path->status)) { heat_path_release(pheat_path); + pheat_path = NULL; } else { /* Register the sampled path */ res = estimator_add_and_release_heat_path(estimator, pheat_path); if(res != RES_OK) goto error; + pheat_path = NULL; } } @@ -155,6 +157,7 @@ solve_pixel estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); exit: + if(pheat_path) heat_path_release(pheat_path); return res; error: goto exit; @@ -220,195 +223,119 @@ error: res_T sdis_solve_probe (struct sdis_scene* scn, - const size_t nrealisations, - const double position[3], - const double time_range[2], - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double Tarad, /* Ambient radiative temperature */ - const double Tref, /* Reference temperature */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_2d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_2d(scn, args, NULL, out_estimator); } else { - return solve_probe_3d(scn, nrealisations, position, time_range, - fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_3d(scn, args, NULL, out_estimator); } } res_T sdis_solve_probe_green_function (struct sdis_scene* scn, - const size_t nrealisations, - const double position[3], - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double Tarad, /* Ambient radiative temperature */ - const double Tref, /* Reference temperature */ + const struct sdis_solve_probe_args* args, struct sdis_green_function** out_green) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_2d(scn, nrealisations, position, NULL, - fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, out_green, NULL); + return solve_probe_2d(scn, args, out_green, NULL); } else { - return solve_probe_3d(scn, nrealisations, position, NULL, - fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, out_green, NULL); + return solve_probe_3d(scn, args, out_green, NULL); } } res_T sdis_solve_probe_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_2d(scn, nrealisations, iprim, uv, time_range, - side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_boundary_2d(scn, args, NULL, out_estimator); } else { - return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, time_range, - side, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_probe_boundary_3d(scn, args, NULL, out_estimator); } } res_T sdis_solve_probe_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** green) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_2d(scn, nrealisations, iprim, uv, NULL, - side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_probe_boundary_2d(scn, args, green, NULL); } else { - return solve_probe_boundary_3d(scn, nrealisations, iprim, uv, NULL, - side, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_probe_boundary_3d(scn, args, green, NULL); } } res_T sdis_solve_boundary (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const enum sdis_side sides[], /* Per primitive side to consider */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_boundary_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_boundary_2d(scn, nrealisations, primitives, sides, nprimitives, - time_range, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_boundary_2d(scn, args, NULL, out_estimator); } else { - return solve_boundary_3d(scn, nrealisations, primitives, sides, nprimitives, - time_range, fp_to_meter, Tarad, Tref, register_paths, NULL, out_estimator); + return solve_boundary_3d(scn, args, NULL, out_estimator); } } res_T sdis_solve_boundary_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const enum sdis_side sides[], /* Per primitive side to consider */ - const size_t nprimitives, /* #primitives */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_boundary_args* args, struct sdis_green_function** green) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_boundary_2d(scn, nrealisations, primitives, sides, - nprimitives, NULL, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, - NULL); + return solve_boundary_2d(scn, args, green, NULL); } else { - return solve_boundary_3d(scn, nrealisations, primitives, sides, - nprimitives, NULL, fp_to_meter, Tarad, Tref, SDIS_HEAT_PATH_NONE, green, - NULL); + return solve_boundary_3d(scn, args, green, NULL); } } res_T sdis_solve_probe_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_probe_boundary_flux_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_probe_boundary_flux_2d(scn, nrealisations, iprim, uv, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_probe_boundary_flux_2d(scn, args, out_estimator); } else { - return solve_probe_boundary_flux_3d(scn, nrealisations, iprim, uv, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_probe_boundary_flux_3d(scn, args, out_estimator); } } res_T sdis_solve_boundary_flux (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_boundary_flux_args* args, struct sdis_estimator** out_estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_boundary_flux_2d(scn, nrealisations, primitives, nprimitives, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_boundary_flux_2d(scn, args, out_estimator); } else { - return solve_boundary_flux_3d(scn, nrealisations, primitives, nprimitives, - time_range, fp_to_meter, Tarad, Tref, out_estimator); + return solve_boundary_flux_3d(scn, args, out_estimator); } } res_T sdis_solve_camera (struct sdis_scene* scn, - const struct sdis_camera* cam, - const double time_range[2], - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const size_t width, /* #pixels in X */ - const size_t height, /* #pixels in Y */ - const size_t spp, /* #samples per pixel */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_camera_args* args, struct sdis_estimator_buffer** out_buf) { #define TILE_SIZE 32 /* definition in X & Y of a tile */ @@ -431,28 +358,36 @@ sdis_solve_camera ATOMIC nsolved_tiles = 0; ATOMIC res = RES_OK; - if(!scn || !cam || fp_to_meter <= 0 || Tref < 0 || !width || !height || !spp - || !out_buf) { + if(!scn + || !args + || !out_buf + || !args->cam + || args->fp_to_meter <= 0 + || !args->image_resolution[0] + || !args->image_resolution[1] + || !args->spp + || args->ambient_radiative_temperature < 0 + || args->reference_temperature < 0 + || args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } + if(scene_is_2d(scn)) { log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME); goto error; } - if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { - res = RES_BAD_ARG; - goto error; - } /* Retrieve the medium in which the submitted position lies */ - res = scene_get_medium(scn, cam->position, NULL, &medium); + res = scene_get_medium(scn, args->cam->position, NULL, &medium); if(res != RES_OK) goto error; if(medium->type != SDIS_FLUID) { log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n", - FUNC_NAME, SPLIT3(cam->position)); + FUNC_NAME, SPLIT3(args->cam->position)); res = RES_BAD_ARG; goto error; } @@ -474,16 +409,17 @@ sdis_solve_camera if(res != RES_OK) goto error; } - ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; - ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; + ntiles_x = (args->image_resolution[0] + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; + ntiles_y = (args->image_resolution[1] + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE; ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y)); ntiles *= ntiles; - pix_sz[0] = 1.0 / (double)width; - pix_sz[1] = 1.0 / (double)height; + pix_sz[0] = 1.0 / (double)args->image_resolution[0]; + pix_sz[1] = 1.0 / (double)args->image_resolution[1]; /* Create the global estimator */ - res = estimator_buffer_create(scn->dev, width, height, &buf); + res = estimator_buffer_create + (scn->dev, args->image_resolution[0], args->image_resolution[1], &buf); if(res != RES_OK) goto error; omp_set_num_threads((int)scn->dev->nthreads); @@ -507,12 +443,14 @@ sdis_solve_camera /* Setup the tile coordinates in the image plane */ tile_org[0] *= TILE_SIZE; tile_org[1] *= TILE_SIZE; - tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]); - tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]); + tile_sz[0] = MMIN(TILE_SIZE, args->image_resolution[0] - tile_org[0]); + tile_sz[1] = MMIN(TILE_SIZE, args->image_resolution[1] - tile_org[1]); /* Draw the tile */ - res_local = solve_tile(scn, rng, medium, cam, time_range, fp_to_meter, - Tarad, Tref, tile_org, tile_sz, spp, register_paths, pix_sz, buf); + res_local = solve_tile(scn, rng, medium, args->cam, args->time_range, + args->fp_to_meter, args->ambient_radiative_temperature, + args->reference_temperature, tile_org, tile_sz, args->spp, + args->register_paths, pix_sz, buf); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; @@ -536,8 +474,8 @@ sdis_solve_camera acc_temp = ACCUM_NULL; acc_time = ACCUM_NULL; nsuccesses = 0; - FOR_EACH(iy, 0, height) { - FOR_EACH(ix, 0, width) { + FOR_EACH(iy, 0, args->image_resolution[1]) { + FOR_EACH(ix, 0, args->image_resolution[0]) { const struct sdis_estimator* estimator; SDIS(estimator_buffer_at(buf, ix, iy, &estimator)); acc_temp.sum += estimator->temperature.sum; @@ -550,12 +488,14 @@ sdis_solve_camera } } - nrealisations = width*height*spp; + nrealisations = args->image_resolution[0]*args->image_resolution[1]*args->spp; ASSERT(acc_temp.count == acc_time.count); ASSERT(acc_temp.count == nsuccesses); estimator_buffer_setup_realisations_count(buf, nrealisations, nsuccesses); estimator_buffer_setup_temperature(buf, acc_temp.sum, acc_temp.sum2); estimator_buffer_setup_realisation_time(buf, acc_time.sum, acc_time.sum2); + res = estimator_buffer_save_rng_state(buf, rng_proxy); + if(res != RES_OK) goto error; exit: if(rngs) { @@ -568,48 +508,35 @@ exit: if(out_buf) *out_buf = buf; return (res_T)res; error: + if(buf) SDIS(estimator_buffer_ref_put(buf)); goto exit; } res_T sdis_solve_medium (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - struct sdis_medium* medium, /* Medium to solve */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_medium_args* args, struct sdis_estimator** estimator) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_medium_2d(scn, nrealisations, medium, time_range, fp_to_meter, - Tarad, Tref, register_paths, NULL, estimator); + return solve_medium_2d(scn, args, NULL, estimator); } else { - return solve_medium_3d(scn, nrealisations, medium, time_range, fp_to_meter, - Tarad, Tref, register_paths, NULL, estimator); + return solve_medium_3d(scn, args, NULL, estimator); } } res_T sdis_solve_medium_green_function (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - struct sdis_medium* medium, /* Medium to solve */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_medium_args* args, struct sdis_green_function** green) { if(!scn) return RES_BAD_ARG; if(scene_is_2d(scn)) { - return solve_medium_2d(scn, nrealisations, medium, NULL, fp_to_meter, Tarad, - Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_medium_2d(scn, args, green, NULL); } else { - return solve_medium_3d(scn, nrealisations, medium, NULL, fp_to_meter, Tarad, - Tref, SDIS_HEAT_PATH_NONE, green, NULL); + return solve_medium_3d(scn, args, green, NULL); } } diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -76,15 +76,7 @@ XD(boundary_get_position)(const unsigned ivert, float pos[DIM], void* context) static res_T XD(solve_boundary) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const enum sdis_side sides[], /* Per primitive side to consider */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_boundary_args* args, struct sdis_green_function** out_green, struct sdis_estimator** out_estimator) { @@ -100,15 +92,18 @@ XD(solve_boundary) struct ssp_rng** rngs = NULL; struct accum* acc_temps = NULL; struct accum* acc_times = NULL; + size_t nrealisations = 0; + int64_t irealisation = 0; size_t i; size_t view_nprims; - int64_t irealisation; int progress = 0; + int register_paths = SDIS_HEAT_PATH_NONE; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives - || !sides || !nprimitives || fp_to_meter <= 0 || Tref < 0) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || !args->primitives || !args->sides || !args->nprimitives + || args->fp_to_meter <= 0) { res = RES_BAD_ARG; goto error; } @@ -117,8 +112,10 @@ XD(solve_boundary) goto error; } if(out_estimator) { - if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + if(args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } @@ -131,16 +128,23 @@ XD(solve_boundary) #endif SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims)); - FOR_EACH(i, 0, nprimitives) { - if(primitives[i] >= view_nprims) { + FOR_EACH(i, 0, args->nprimitives) { + if(args->primitives[i] >= view_nprims) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)primitives[i], + (unsigned long)args->primitives[i], (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; } + if((unsigned)args->sides[i] >= SDIS_SIDE_NULL__) { + log_err(scn->dev, + "%s: invalid side for the primitive `%lu'.\n", + FUNC_NAME, (unsigned long)args->primitives[i]); + res = RES_BAD_ARG; + goto error; + } } /* Create the Star-XD shape of the boundary */ @@ -153,18 +157,20 @@ XD(solve_boundary) /* Initialise the boundary shape with the triangles/segments of the * submitted primitives */ - ctx.primitives = primitives; + ctx.primitives = args->primitives; ctx.view = scn->sXd(view); vdata.usage = SXD_POSITION; vdata.get = XD(boundary_get_position); #if SDIS_XD_DIMENSION == 2 vdata.type = S2D_FLOAT2; - res = s2d_line_segments_setup_indexed_vertices(shape, (unsigned)nprimitives, - boundary_get_indices_2d, (unsigned)(nprimitives*2), &vdata, 1, &ctx); + res = s2d_line_segments_setup_indexed_vertices(shape, + (unsigned)args->nprimitives, boundary_get_indices_2d, + (unsigned)(args->nprimitives*2), &vdata, 1, &ctx); #else vdata.type = S3D_FLOAT3; - res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprimitives, - boundary_get_indices_3d, (unsigned)(nprimitives*3), &vdata, 1, &ctx); + res = s3d_mesh_setup_indexed_vertices(shape, + (unsigned)args->nprimitives, boundary_get_indices_3d, + (unsigned)(args->nprimitives*3), &vdata, 1, &ctx); #endif if(res != RES_OK) goto error; @@ -177,9 +183,15 @@ XD(solve_boundary) if(res != RES_OK) goto error; /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); @@ -213,6 +225,8 @@ XD(solve_boundary) if(res != RES_OK) goto error; } + nrealisations = args->nrealisations; + register_paths = out_estimator ? args->register_paths : SDIS_HEAT_PATH_NONE; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) { @@ -243,7 +257,7 @@ XD(solve_boundary) time_current(&t0); if(!out_green) { - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); if(register_paths) { heat_path_init(scn->dev->allocator, &heat_path); pheat_path = &heat_path; @@ -253,7 +267,10 @@ XD(solve_boundary) * function. Simply takes 0 as relative time */ time = 0; res_local = green_function_create_path(greens[ithread], &green_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } pgreen_path = &green_path; } @@ -274,21 +291,25 @@ XD(solve_boundary) &prim, st); d2_set_f2(uv, st); #endif - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } /* Map from boundary scene to sdis scene */ - ASSERT(prim.prim_id < nprimitives); - iprim = primitives[prim.prim_id]; - side = sides[prim.prim_id]; + ASSERT(prim.prim_id < args->nprimitives); + iprim = args->primitives[prim.prim_id]; + side = args->sides[prim.prim_id]; /* Invoke the boundary realisation */ res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side, - fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w); + args->fp_to_meter, args->ambient_radiative_temperature, + args->reference_temperature, pgreen_path, pheat_path, &w); /* Fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { ATOMIC_SET(&res, res_simul); - continue; + goto error_it; } /* Register heat path */ @@ -300,9 +321,14 @@ XD(solve_boundary) /* Check if the path must be saved regarding the register_paths mask */ if(!(register_paths & (int)pheat_path->status)) { heat_path_release(pheat_path); + pheat_path = NULL; } else { /* Register the sampled path */ res_local = estimator_add_and_release_heat_path(estimator, pheat_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } + pheat_path = NULL; } } @@ -324,6 +350,11 @@ XD(solve_boundary) progress = pcent; log_info(scn->dev, "Solving boundary temperature: %3d%%\r", progress); } + exit_it: + if(pheat_path) heat_path_release(pheat_path); + continue; + error_it: + goto exit_it; } if(res != RES_OK) goto error; @@ -342,6 +373,8 @@ XD(solve_boundary) estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; } /* Setup the green function */ @@ -397,13 +430,7 @@ error: static res_T XD(solve_boundary_flux) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t primitives[], /* List of boundary primitives to handle */ - const size_t nprimitives, /* #primitives */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_boundary_flux_args* args, struct sdis_estimator** out_estimator) { struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL); @@ -419,17 +446,24 @@ XD(solve_boundary_flux) struct accum* acc_fl = NULL; /* Per thread flux accumulator */ struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */ struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */ + size_t nrealisations = 0; + int64_t irealisation; size_t i; size_t view_nprims; - int64_t irealisation; int progress = 0; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives - || !time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1]) - || !nprimitives || fp_to_meter < 0 || Tref < 0 + if(!scn + || !args + || !args->nrealisations + || args->nrealisations > INT64_MAX + || !args->primitives + || args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || (args->time_range[1] > DBL_MAX && args->time_range[0] != args->time_range[1]) + || !args->nprimitives + || args->fp_to_meter < 0 || !out_estimator) { res = RES_BAD_ARG; goto error; @@ -442,12 +476,12 @@ XD(solve_boundary_flux) #endif SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims)); - FOR_EACH(i, 0, nprimitives) { - if(primitives[i] >= view_nprims) { + FOR_EACH(i, 0, args->nprimitives) { + if(args->primitives[i] >= view_nprims) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)primitives[i], + (unsigned long)args->primitives[i], (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; @@ -464,19 +498,20 @@ XD(solve_boundary_flux) /* Initialise the boundary shape with the triangles/segments of the * submitted primitives */ - ctx.primitives = primitives; + ctx.primitives = args->primitives; ctx.view = scn->sXd(view); vdata.get = XD(boundary_get_position); #if SDIS_XD_DIMENSION == 2 vdata.usage = S2D_POSITION; vdata.type = S2D_FLOAT2; - res = s2d_line_segments_setup_indexed_vertices(shape, (unsigned)nprimitives, - XD(boundary_get_indices), (unsigned)(nprimitives*2), &vdata, 1, &ctx); + res = s2d_line_segments_setup_indexed_vertices(shape, + (unsigned)args->nprimitives, XD(boundary_get_indices), + (unsigned)(args->nprimitives*2), &vdata, 1, &ctx); #else /* DIM == 3 */ vdata.usage = S3D_POSITION; vdata.type = S3D_FLOAT3; - res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprimitives, - XD(boundary_get_indices), (unsigned)(nprimitives*3), &vdata, 1, &ctx); + res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)args->nprimitives, + XD(boundary_get_indices), (unsigned)(args->nprimitives*3), &vdata, 1, &ctx); #endif if(res != RES_OK) goto error; @@ -489,9 +524,15 @@ XD(solve_boundary_flux) if(res != RES_OK) goto error; /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); @@ -517,6 +558,7 @@ XD(solve_boundary_flux) res = estimator_create(scn->dev, SDIS_ESTIMATOR_FLUX, &estimator); if(res != RES_OK) goto error; + nrealisations = args->nrealisations; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -534,6 +576,8 @@ XD(solve_boundary_flux) const struct sdis_medium *fmd, *bmd; enum sdis_side solid_side, fluid_side; double T_brf[3] = { 0, 0, 0 }; + const double Tref = args->reference_temperature; + const double Tarad = args->ambient_radiative_temperature; double epsilon, hc, hr; size_t iprim; double uv[DIM - 1]; @@ -550,7 +594,7 @@ XD(solve_boundary_flux) /* Begin time registration */ time_current(&t0); - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); /* Sample a position onto the boundary */ #if SDIS_XD_DIMENSION == 2 @@ -572,8 +616,8 @@ XD(solve_boundary_flux) if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } /* Map from boundary scene to sdis scene */ - ASSERT(prim.prim_id < nprimitives); - iprim = primitives[prim.prim_id]; + ASSERT(prim.prim_id < args->nprimitives); + iprim = args->primitives[prim.prim_id]; interf = scene_get_interface(scn, (unsigned)iprim); fmd = interface_get_medium(interf, SDIS_FRONT); @@ -605,7 +649,7 @@ XD(solve_boundary_flux) if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE; res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time, - solid_side, fp_to_meter, Tarad, Tref, flux_mask, T_brf); + solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); @@ -676,6 +720,9 @@ XD(solve_boundary_flux) estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2); estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; + exit: if(rngs) { FOR_EACH(i, 0, scn->dev->nthreads) { if(rngs[i]) SSP(rng_ref_put(rngs[i])); } diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h @@ -200,13 +200,7 @@ error: static res_T XD(solve_medium) (struct sdis_scene* scn, - const size_t nrealisations, - struct sdis_medium* mdm, - const double time_range[2], - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double Tarad, /* Ambient radiative temperature */ - const double Tref, /* Reference temperature */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_medium_args* args, struct sdis_green_function** out_green, /* May be NULL <=> No green func */ struct sdis_estimator** out_estimator) /* May be NULL <=> No estimator */ { @@ -218,15 +212,17 @@ XD(solve_medium) struct sdis_estimator* estimator = NULL; struct accum* acc_temps = NULL; struct accum* acc_times = NULL; + size_t nrealisations = 0; int64_t irealisation; int cumul_is_init = 0; size_t i; int progress = 0; + int register_paths = SDIS_HEAT_PATH_NONE; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !mdm || !nrealisations || nrealisations > INT64_MAX - || fp_to_meter <= 0 || Tref < 0) { + if(!scn || !args || !args->medium || !args->nrealisations + || args->nrealisations > INT64_MAX || args->fp_to_meter <= 0) { res = RES_BAD_ARG; goto error; } @@ -235,8 +231,10 @@ XD(solve_medium) goto error; } if(out_estimator) { - if(!time_range || time_range[0] < 0 || time_range[0] > time_range[1] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + if(args->time_range[0] < 0 + || args->time_range[0] > args->time_range[1] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } @@ -249,9 +247,15 @@ XD(solve_medium) #endif /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); @@ -272,7 +276,7 @@ XD(solve_medium) /* Compute the enclosure cumulative */ darray_enclosure_cumul_init(scn->dev->allocator, &cumul); cumul_is_init = 1; - res = compute_medium_enclosure_cumulative(scn, mdm, &cumul); + res = compute_medium_enclosure_cumulative(scn, args->medium, &cumul); if(res != RES_OK) goto error; if(out_green) { @@ -291,6 +295,8 @@ XD(solve_medium) if(res != RES_OK) goto error; } + nrealisations = args->nrealisations; + register_paths = out_estimator ? args->register_paths : SDIS_HEAT_PATH_NONE; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -318,7 +324,7 @@ XD(solve_medium) if(!out_green) { /* Sample the time */ - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); /* Prepare path registration if necessary */ if(register_paths) { @@ -330,7 +336,10 @@ XD(solve_medium) * function. Simply takes 0 as relative time */ time = 0; res_local = green_function_create_path(greens[ithread], &green_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } pgreen_path = &green_path; } @@ -342,16 +351,18 @@ XD(solve_medium) if(res_local != RES_OK) { log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME); ATOMIC_SET(&res, res_local); - continue; + goto error_it; } /* Run a probe realisation */ - res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, mdm, pos, - time, fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &weight); + res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, + args->medium, pos, time, args->fp_to_meter, + args->ambient_radiative_temperature, args->reference_temperature, + pgreen_path, pheat_path, &weight); if(res_simul != RES_OK && res_simul != RES_BAD_OP) { ATOMIC_SET(&res, res_simul); - continue; + goto error_it; } /* Finalize the registered path */ @@ -363,10 +374,15 @@ XD(solve_medium) /* Check if the path must be saved regarding the register_paths mask */ if(!(register_paths & (int)pheat_path->status)) { heat_path_release(pheat_path); + pheat_path = NULL; } else { /* Register the sampled path */ res_local = estimator_add_and_release_heat_path(estimator, pheat_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } } + pheat_path = NULL; } /* Stop time registration */ @@ -387,6 +403,11 @@ XD(solve_medium) progress = pcent; log_info(scn->dev, "Solving medium temperature: %3d%%\r", progress); } + exit_it: + if(pheat_path) heat_path_release(pheat_path); + continue; + error_it: + goto exit_it; } if(res != RES_OK) goto error; @@ -405,6 +426,8 @@ XD(solve_medium) estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; } if(out_green) { diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h @@ -33,13 +33,7 @@ static res_T XD(solve_probe) (struct sdis_scene* scn, - const size_t nrealisations, - const double position[3], - const double time_range[2], - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double Tarad, /* Ambient radiative temperature */ - const double Tref, /* Reference temperature */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_args* args, struct sdis_green_function** out_green, /* May be NULL <=> No green func */ struct sdis_estimator** out_estimator) /* May be NULL <=> No estimator */ { @@ -51,14 +45,16 @@ XD(solve_probe) struct ssp_rng** rngs = NULL; struct accum* acc_temps = NULL; struct accum* acc_times = NULL; + size_t nrealisations = 0; int64_t irealisation = 0; size_t i; int progress = 0; + int register_paths = SDIS_HEAT_PATH_NONE; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !position - || fp_to_meter <= 0 || Tref < 0) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || args->fp_to_meter <= 0) { res = RES_BAD_ARG; goto error; } @@ -67,13 +63,16 @@ XD(solve_probe) goto error; } if(out_estimator) { - if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + if(args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } } + #if SDIS_XD_DIMENSION == 2 if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; } #else @@ -81,9 +80,15 @@ XD(solve_probe) #endif /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs)); @@ -102,7 +107,7 @@ XD(solve_probe) if(!acc_times) { res = RES_MEM_ERR; goto error; } /* Retrieve the medium in which the submitted position lies */ - res = scene_get_medium(scn, position, NULL, &medium); + res = scene_get_medium(scn, args->position, NULL, &medium); if(res != RES_OK) goto error; /* Create the per thread green function */ @@ -122,6 +127,8 @@ XD(solve_probe) } /* Here we go! Launch the Monte Carlo estimation */ + nrealisations = args->nrealisations; + register_paths = out_estimator ? args->register_paths : SDIS_HEAT_PATH_NONE; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -147,7 +154,7 @@ XD(solve_probe) time_current(&t0); if(!out_green) { - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); if(register_paths) { heat_path_init(scn->dev->allocator, &heat_path); pheat_path = &heat_path; @@ -157,18 +164,22 @@ XD(solve_probe) * function. Simply takes 0 as relative time */ time = 0; res_local = green_function_create_path(greens[ithread], &green_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } - + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } pgreen_path = &green_path; } res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, medium, - position, time, fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w); + args->position, time, args->fp_to_meter, + args->ambient_radiative_temperature, args->reference_temperature, + pgreen_path, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { ATOMIC_SET(&res, res_simul); - continue; + goto error_it; } if(pheat_path) { @@ -179,9 +190,14 @@ XD(solve_probe) /* Check if the path must be saved regarding the register_paths mask */ if(!(register_paths & (int)pheat_path->status)) { heat_path_release(pheat_path); + pheat_path = NULL; } else { /* Register the sampled path */ res_local = estimator_add_and_release_heat_path(estimator, pheat_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } + pheat_path = NULL; } } @@ -203,6 +219,11 @@ XD(solve_probe) progress = pcent; log_info(scn->dev, "Solving probe temperature: %3d%%\r", progress); } + exit_it: + if(pheat_path) heat_path_release(pheat_path); + continue; + error_it: + goto exit_it; } if(res != RES_OK) goto error; @@ -221,6 +242,8 @@ XD(solve_probe) estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; } if(out_green) { diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -33,15 +33,7 @@ static res_T XD(solve_probe_boundary) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const enum sdis_side side, /* Side of iprim on which the probe lies */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ - const int register_paths, /* Combination of enum sdis_heat_path_flag */ + const struct sdis_solve_probe_boundary_args* args, struct sdis_green_function** out_green, struct sdis_estimator** out_estimator) { @@ -52,14 +44,16 @@ XD(solve_probe_boundary) struct ssp_rng** rngs = NULL; struct accum* acc_temps = NULL; struct accum* acc_times = NULL; + size_t nrealisations = 0; int64_t irealisation = 0; size_t i; int progress = 0; + int register_paths = SDIS_HEAT_PATH_NONE; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv - || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || args->fp_to_meter <= 0 || ((unsigned)args->side >= SDIS_SIDE_NULL__)) { res = RES_BAD_ARG; goto error; } @@ -68,8 +62,10 @@ XD(solve_probe_boundary) goto error; } if(out_estimator) { - if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) { + if(args->time_range[0] < 0 + || args->time_range[1] < args->time_range[0] + || ( args->time_range[1] > DBL_MAX + && args->time_range[0] != args->time_range[1])) { res = RES_BAD_ARG; goto error; } @@ -82,12 +78,12 @@ XD(solve_probe_boundary) #endif /* Check the primitive identifier */ - if(iprim >= scene_get_primitives_count(scn)) { + if(args->iprim >= scene_get_primitives_count(scn)) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. " "It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)iprim, + (unsigned long)args->iprim, (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; @@ -96,25 +92,25 @@ XD(solve_probe_boundary) /* Check parametric coordinates */ #if SDIS_XD_DIMENSION == 2 { - const double v = CLAMP(1.0 - uv[0], 0, 1); - if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) { + const double v = CLAMP(1.0 - args->uv[0], 0, 1); + if(args->uv[0] < 0 || args->uv[0] > 1 || !eq_eps(args->uv[0]+v, 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates %g." "u + (1-u) must be equal to 1 with u [0, 1].\n", - FUNC_NAME, uv[0]); + FUNC_NAME, args->uv[0]); res = RES_BAD_ARG; goto error; } } #else /* SDIS_XD_DIMENSION == 3 */ { - const double w = CLAMP(1 - uv[0] - uv[1], 0, 1); - if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1 - || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) { + const double w = CLAMP(1 - args->uv[0] - args->uv[1], 0, 1); + if(args->uv[0] < 0 || args->uv[1] < 0 || args->uv[0] > 1 || args->uv[1] > 1 + || !eq_eps(w + args->uv[0] + args->uv[1], 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates [%g, %g]. " "u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n", - FUNC_NAME, uv[0], uv[1]); + FUNC_NAME, args->uv[0], args->uv[1]); res = RES_BAD_ARG; goto error; } @@ -122,9 +118,15 @@ XD(solve_probe_boundary) #endif /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC @@ -160,6 +162,8 @@ XD(solve_probe_boundary) } /* Here we go! Launch the Monte Carlo estimation */ + nrealisations = args->nrealisations; + register_paths = out_estimator ? args->register_paths : SDIS_HEAT_PATH_NONE; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -185,7 +189,7 @@ XD(solve_probe_boundary) time_current(&t0); if(!out_green) { - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); if(register_paths) { heat_path_init(scn->dev->allocator, &heat_path); pheat_path = &heat_path; @@ -195,17 +199,21 @@ XD(solve_probe_boundary) * function. Simply takes 0 as relative time */ time = 0; res_local = green_function_create_path(greens[ithread], &green_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } pgreen_path = &green_path; } - res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side, - fp_to_meter, Tarad, Tref, pgreen_path, pheat_path, &w); + res_simul = XD(boundary_realisation)(scn, rng, args->iprim, args->uv, time, + args->side, args->fp_to_meter, args->ambient_radiative_temperature, + args->reference_temperature, pgreen_path, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { ATOMIC_SET(&res, res_simul); - continue; + goto error_it; } if(pheat_path) { @@ -216,9 +224,15 @@ XD(solve_probe_boundary) /* Check if the path must be saved regarding the register_paths mask */ if(!(register_paths & (int)pheat_path->status)) { heat_path_release(pheat_path); - } else { /* Register the sampled path */ + pheat_path = NULL; + } else { + /* Register the sampled path */ res_local = estimator_add_and_release_heat_path(estimator, pheat_path); - if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + if(res_local != RES_OK) { + ATOMIC_SET(&res, res_local); + goto error_it; + } + pheat_path = NULL; } } @@ -240,6 +254,12 @@ XD(solve_probe_boundary) progress = pcent; log_info(scn->dev, "Solving probe boundary temperature: %3d%%\r", progress); } + + exit_it: + if(pheat_path) heat_path_release(pheat_path); + continue; + error_it: + goto exit_it; } if(res != RES_OK) goto error; @@ -258,6 +278,8 @@ XD(solve_probe_boundary) estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count); estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; } if(out_green) { @@ -309,13 +331,7 @@ error: static res_T XD(solve_probe_boundary_flux) (struct sdis_scene* scn, - const size_t nrealisations, /* #realisations */ - const size_t iprim, /* Identifier of the primitive on which the probe lies */ - const double uv[2], /* Parametric coordinates of the probe onto the primitve */ - const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ + const struct sdis_solve_probe_boundary_flux_args* args, struct sdis_estimator** out_estimator) { struct sdis_estimator* estimator = NULL; @@ -330,17 +346,17 @@ XD(solve_probe_boundary_flux) struct accum* acc_fl = NULL; /* Per thread flux accumulator */ struct accum* acc_fc = NULL; /* Per thread convective flux accumulator */ struct accum* acc_fr = NULL; /* Per thread radiative flux accumulator */ + size_t nrealisations = 0; int64_t irealisation = 0; size_t i; int progress = 0; ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv - || !time_range || time_range[0] < 0 || time_range[1] < time_range[0] - || (time_range[1] > DBL_MAX && time_range[0] != time_range[1]) - || fp_to_meter <= 0 || Tref < 0 - || !out_estimator) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX + || args->time_range[0] < 0 || args->time_range[1] < args->time_range[0] + || (args->time_range[1]>DBL_MAX && args->time_range[0] != args->time_range[1]) + || args->fp_to_meter <= 0 || !out_estimator) { res = RES_BAD_ARG; goto error; } @@ -352,12 +368,12 @@ XD(solve_probe_boundary_flux) #endif /* Check the primitive identifier */ - if(iprim >= scene_get_primitives_count(scn)) { + if(args->iprim >= scene_get_primitives_count(scn)) { log_err(scn->dev, "%s: invalid primitive identifier `%lu'. " "It must be in the [0 %lu] range.\n", FUNC_NAME, - (unsigned long)iprim, + (unsigned long)args->iprim, (unsigned long)scene_get_primitives_count(scn)-1); res = RES_BAD_ARG; goto error; @@ -365,29 +381,33 @@ XD(solve_probe_boundary_flux) /* Check parametric coordinates */ if(scene_is_2d(scn)) { - const double v = CLAMP(1.0 - uv[0], 0, 1); - if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) { + const double v = CLAMP(1.0 - args->uv[0], 0, 1); + if(args->uv[0] < 0 || args->uv[0] > 1 + || !eq_eps(args->uv[0] + v, 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates %g. " "u + (1-u) must be equal to 1 with u [0, 1].\n", - FUNC_NAME, uv[0]); + FUNC_NAME, args->uv[0]); res = RES_BAD_ARG; goto error; } } else { - const double w = CLAMP(1 - uv[0] - uv[1], 0, 1); - if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1 - || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) { + const double w = CLAMP(1 - args->uv[0] - args->uv[1], 0, 1); + if(args->uv[0] < 0 + || args->uv[1] < 0 + || args->uv[0] > 1 + || args->uv[1] > 1 + || !eq_eps(w + args->uv[0] + args->uv[1], 1, 1.e-6)) { log_err(scn->dev, "%s: invalid parametric coordinates [%g, %g]. " "u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n", - FUNC_NAME, uv[0], uv[1]); + FUNC_NAME, args->uv[0], args->uv[1]); res = RES_BAD_ARG; goto error; } } /* Check medium is fluid on one side and solid on the other */ - interf = scene_get_interface(scn, (unsigned)iprim); + interf = scene_get_interface(scn, (unsigned)args->iprim); fmd = interface_get_medium(interf, SDIS_FRONT); bmd = interface_get_medium(interf, SDIS_BACK); if(!fmd || !bmd @@ -400,9 +420,15 @@ XD(solve_probe_boundary_flux) fluid_side = (fmd->type == SDIS_FLUID) ? SDIS_FRONT : SDIS_BACK; /* Create the proxy RNG */ - res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, - scn->dev->nthreads, &rng_proxy); - if(res != RES_OK) goto error; + if(args->rng_state) { + res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } else { + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + } /* Create the per thread RNG */ rngs = MEM_CALLOC @@ -427,7 +453,7 @@ XD(solve_probe_boundary_flux) /* Prebuild the interface fragment */ res = XD(build_interface_fragment) - (&frag, scn, (unsigned)iprim, uv, fluid_side); + (&frag, scn, (unsigned)args->iprim, args->uv, fluid_side); if(res != RES_OK) goto error; /* Create the estimator */ @@ -435,6 +461,7 @@ XD(solve_probe_boundary_flux) if(res != RES_OK) goto error; /* Here we go! Launch the Monte Carlo estimation */ + nrealisations = args->nrealisations; omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { @@ -449,6 +476,8 @@ XD(solve_probe_boundary_flux) double time, epsilon, hc, hr; int flux_mask = 0; double T_brf[3] = { 0, 0, 0 }; + const double Tref = args->reference_temperature; + const double Tarad = args->ambient_radiative_temperature; size_t n; int pcent; res_T res_simul = RES_OK; @@ -458,7 +487,7 @@ XD(solve_probe_boundary_flux) /* Begin time registration */ time_current(&t0); - time = sample_time(rng, time_range); + time = sample_time(rng, args->time_range); /* Compute hr and hc */ frag.time = time; @@ -470,8 +499,8 @@ XD(solve_probe_boundary_flux) flux_mask = 0; if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE; if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE; - res_simul = XD(boundary_flux_realisation)(scn, rng, iprim, uv, time, - solid_side, fp_to_meter, Tarad, Tref, flux_mask, T_brf); + res_simul = XD(boundary_flux_realisation)(scn, rng, args->iprim, args->uv, time, + solid_side, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); @@ -542,6 +571,9 @@ XD(solve_probe_boundary_flux) estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr[0].sum, acc_fr[0].sum2); estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl[0].sum, acc_fl[0].sum2); + res = estimator_save_rng_state(estimator, rng_proxy); + if(res != RES_OK) goto error; + exit: if(rngs) { FOR_EACH(i, 0, scn->dev->nthreads) {if(rngs[i]) SSP(rng_ref_put(rngs[i]));} diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -388,27 +388,30 @@ main(int argc, char** argv) struct sdis_estimator* estimator; struct sdis_estimator* estimator2; struct sdis_green_function* green; - double pos[3]; - double time_range[2] = { INF, INF }; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; double ref, u; size_t nreals = 0; size_t nfails = 0; const size_t N = 10000; - pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); - pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); - pos[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.nrealisations = N; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.position[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.reference_temperature = Tref; - OK(sdis_solve_probe(scn, N, pos, time_range, 1, -1, Tref, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); - u = (pos[0] + 1) / thickness; + u = (solve_args.position[0] + 1) / thickness; ref = u * Ts1 + (1-u) * Ts0; printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + SPLIT3(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -417,8 +420,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); /* Check green function */ - OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, Tref, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -426,8 +429,10 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); - OK(sdis_solve_probe(scn, 10, pos, time_range, 1, -1, Tref, - SDIS_HEAT_PATH_ALL, &estimator)); + solve_args.nrealisations = 10; + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); printf("\n"); diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -394,26 +394,29 @@ main(int argc, char** argv) struct sdis_estimator* estimator; struct sdis_estimator* estimator2; struct sdis_green_function* green; - double pos[2]; - double time_range[2] = { INF, INF }; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; double ref, u; size_t nreals = 0; size_t nfails = 0; const size_t N = 10000; - pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); - pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.nrealisations = N; + solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.reference_temperature = Tref; - OK(sdis_solve_probe(scn, 10000, pos, time_range, 1, -1, Tref, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); - u = (pos[0] + 1) / thickness; + u = (solve_args.position[0] + 1) / thickness; ref = u * Ts1 + (1-u) * Ts0; printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + SPLIT2(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -422,8 +425,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); /* Check green function */ - OK(sdis_solve_probe_green_function(scn, 10000, pos, 1, -1, Tref, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -431,8 +434,10 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); - OK(sdis_solve_probe - (scn, 10, pos, time_range, 1, -1, Tref, SDIS_HEAT_PATH_ALL, &estimator)); + solve_args.nrealisations = 10; + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); printf("\n"); diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -188,7 +188,7 @@ main(int argc, char** argv) struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; struct sdis_interface* box_interfaces[12/*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; - double pos[3]; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; double ref; double Tinf; double nu; @@ -262,23 +262,28 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf_T4)); OK(sdis_interface_ref_put(interf_T5)); - d3_splat(pos, 0.25); + solve_args.nrealisations = N; + solve_args.position[0] = 0.25; + solve_args.position[1] = 0.25; + solve_args.position[2] = 0.25; /* Test in 3D for various time values. */ nu = (6 * H) / (RHO*CP); Tinf = (H*(T0 + T1 + T2 + T3 + T4 + T5)) / (6 * H); - printf(">>> Temperature of the box at (%g %g %g)\n\n", SPLIT3(pos)); + printf(">>> Temperature of the box at (%g %g %g)\n\n", + SPLIT3(solve_args.position)); FOR_EACH(i, 0, 5) { double time = i ? (double)i / nu : INF; - double time_range[2]; - time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); + solve_args.time_range[0] = time; + solve_args.time_range[1] = time; + /* Setup stationary state */ *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(box_scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &mc_time)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -291,8 +296,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE * 3)); if(IS_INF(time)) { /* Check green function */ - OK(sdis_solve_probe_green_function(box_scn, N, pos, 1.0, 0, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); OK(sdis_estimator_ref_put(estimator2)); @@ -306,18 +311,20 @@ main(int argc, char** argv) /* Test in 2D for various time values. */ nu = (4 * H) / (RHO*CP); Tinf = (H * (T0 + T1 + T2 + T3)) / (4 * H); - printf(">>> Temperature of the square at (%g %g)\n\n", SPLIT2(pos)); + printf(">>> Temperature of the square at (%g %g)\n\n", + SPLIT2(solve_args.position)); FOR_EACH(i, 0, 5) { double time = i ? (double)i / nu : INF; - double time_range[2]; - time_range[0] = time_range[1] = time; ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); + solve_args.time_range[0] = time; + solve_args.time_range[1] = time; + /* Setup stationnary state */ *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 2D */ - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(square_scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -330,8 +337,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE * 3)); if(IS_INF(time)) { /* Check green function */ - OK(sdis_solve_probe_green_function(square_scn, N, pos, 1.0, 0, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); OK(sdis_estimator_ref_put(estimator2)); diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -196,9 +196,9 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct sdis_interface* box_interfaces[12/*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; - double pos[3]; double ref; double Tinf; double nu; @@ -278,23 +278,28 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf_T4)); OK(sdis_interface_ref_put(interf_T5)); - d3_splat(pos, 0.25); + solve_args.nrealisations = N; + solve_args.position[0] = 0.25; + solve_args.position[1] = 0.25; + solve_args.position[2] = 0.25; /* Test in 3D for various time values. */ nu = (HC0 + HC1 + HC2 + HC3 + HC4 + HC5) / (RHO * CP); Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3 + HC4 * T4 + HC5 * T5) / (HC0 + HC1 + HC2 + HC3 + HC4 + HC5); - printf(">>> Temperature of the box at (%g %g %g)\n\n", SPLIT3(pos)); + printf(">>> Temperature of the box at (%g %g %g)\n\n", + SPLIT3(solve_args.position)); FOR_EACH(i, 0, 5) { double time = i ? (double)i / nu : INF; - double time_range[2]; - time_range[0] = time_range[1] = time; + solve_args.time_range[0] = time; + solve_args.time_range[1] = time; + ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); *((int*)sdis_data_get(is_stationary)) = IS_INF(time); /* Solve in 3D */ - OK(sdis_solve_probe(box_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(box_scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -307,8 +312,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE * 3)); if(IS_INF(time)) { /* Check green function */ - OK(sdis_solve_probe_green_function(box_scn, N, pos, 1.0, 0, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); OK(sdis_estimator_ref_put(estimator2)); @@ -322,16 +327,19 @@ main(int argc, char** argv) /* Test in 2D for various time values. */ nu = (HC0 + HC1 + HC2 + HC3) / (RHO * CP); Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3) / (HC0 + HC1 + HC2 + HC3); - printf(">>> Temperature of the square at (%g %g)\n\n", SPLIT2(pos)); + printf(">>> Temperature of the square at (%g %g)\n\n", + SPLIT2(solve_args.position)); FOR_EACH(i, 0, 5) { double time = i ? (double)i / nu : INF; - double time_range[2]; - time_range[0] = time_range[1] = time; + + solve_args.time_range[0] = time; + solve_args.time_range[1] = time; + ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); *((int*)sdis_data_get(is_stationary)) = IS_INF(time); - OK(sdis_solve_probe(square_scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(square_scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); CHK(nfails + nreals == N); @@ -344,8 +352,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE * 3)); if(IS_INF(time)) { /* Check green function */ - OK(sdis_solve_probe_green_function(square_scn, N, pos, 1.0, 0, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); OK(sdis_estimator_ref_put(estimator2)); diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -136,6 +136,7 @@ solve(struct sdis_scene* scn, const double pos[]) struct sdis_green_function* green; struct sdis_mc T; struct sdis_mc time; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals; size_t nfails; double ref; @@ -145,8 +146,17 @@ solve(struct sdis_scene* scn, const double pos[]) ref = T0 + (1 - pos[0]) * PHI/LAMBDA; + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + solve_args.position[0] = pos[0]; + solve_args.position[1] = pos[1]; + solve_args.position[2] = dim == SDIS_SCENE_2D ? 0 : pos[2]; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + time_current(&t0); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); @@ -155,8 +165,6 @@ solve(struct sdis_scene* scn, const double pos[]) OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); - OK(sdis_scene_get_dimension(scn, &dim)); - switch(dim) { case SDIS_SCENE_2D: printf("Temperature at (%g %g) = %g ~ %g +/- %g\n", @@ -177,7 +185,7 @@ solve(struct sdis_scene* scn, const double pos[]) CHK(eq_eps(T.E, ref, T.SE*3)); time_current(&t0); - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); time_current(&t1); OK(sdis_green_function_solve(green, time_range, &estimator2)); time_current(&t2); diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c @@ -273,14 +273,14 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_solve_medium_args solve_mdm_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; struct interf* interf_param = NULL; struct solid* solid_param = NULL; struct context ctx; - double probe[3]; double lower[3]; double upper[3]; double size[3]; - const double time[2] = {DBL_MAX, DBL_MAX}; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -347,30 +347,38 @@ main(int argc, char** argv) interf_param->upper[2] = upper[2]; interf_param->h = 1; - probe[0] = 0; - probe[1] = 0; - probe[2] = 0; + solve_args.nrealisations = Nreals; + solve_args.position[0] = 0; + solve_args.position[1] = 0; + solve_args.position[2] = 0; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.register_paths = SDIS_HEAT_PATH_FAILURE; /* Launch probe estimation with trilinear profile set at interfaces */ interf_param->profile = PROFILE_TRILINEAR; - OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0, - SDIS_HEAT_PATH_FAILURE, &estimator)); - print_estimation_result(estimator, trilinear_temperature(probe)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + print_estimation_result + (estimator, trilinear_temperature(solve_args.position)); OK(sdis_estimator_ref_put(estimator)); /* Launch probe estimation with volumetric power profile set at interfaces */ interf_param->profile = PROFILE_VOLUMETRIC_POWER; solid_param->power = Pw; - OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0, - SDIS_HEAT_PATH_FAILURE, &estimator)); - print_estimation_result(estimator, volumetric_temperature(probe, upper)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + print_estimation_result + (estimator, volumetric_temperature(solve_args.position, upper)); solid_param->power = SDIS_VOLUMIC_POWER_NONE; OK(sdis_estimator_ref_put(estimator)); /* Launch medium integration */ interf_param->profile = PROFILE_UNKNOWN; - OK(sdis_solve_medium(scn, Nreals, solid, time, 1.0, -1, 0, - SDIS_HEAT_PATH_FAILURE, &estimator)); + solve_mdm_args.nrealisations = Nreals; + solve_mdm_args.medium = solid; + solve_mdm_args.time_range[0] = INF; + solve_mdm_args.time_range[1] = INF; + OK(sdis_solve_medium(scn, &solve_mdm_args, &estimator)); + print_estimation_result(estimator, Tfluid); /*dump_heat_paths(stdout, estimator);*/ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -191,16 +191,15 @@ main(int argc, char** argv) struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; + struct sdis_solve_probe_boundary_args probe_args = + SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; + struct sdis_solve_boundary_args bound_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; struct interf* interf_props = NULL; struct fluid* fluid_param; - double uv[2]; double pos[3]; - double time_range[2] = { INF, INF }; - double tr[2]; double ref; size_t prims[4]; enum sdis_side sides[4]; - size_t iprim; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -297,41 +296,57 @@ main(int argc, char** argv) #define SOLVE sdis_solve_probe_boundary #define GREEN sdis_solve_probe_boundary_green_function - #define F SDIS_FRONT - uv[0] = 0.3; - uv[1] = 0.3; - iprim = 6; - - BA(SOLVE(NULL, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, 0, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, 12, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, NULL, time_range, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, NULL, F, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, -1, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - tr[0] = 1; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - tr[1] = 0; - BA(SOLVE(box_scn, N, iprim, uv, tr, F, 1.0, 0, 0, 0, NULL)); - - OK(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos)); + + probe_args.nrealisations = N; + probe_args.uv[0] = 0.3; + probe_args.uv[1] = 0.3; + probe_args.iprim = 6; + probe_args.time_range[0] = INF; + probe_args.time_range[1] = INF; + probe_args.side = SDIS_FRONT; + + BA(SOLVE(NULL, &probe_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &probe_args, NULL)); + probe_args.nrealisations = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.nrealisations = N; + probe_args.iprim = 12; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.iprim = 6; + probe_args.side = SDIS_SIDE_NULL__; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.side = SDIS_FRONT; + probe_args.time_range[0] = probe_args.time_range[1] = -1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = 1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[1] = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = probe_args.time_range[1] = INF; + + OK(SOLVE(box_scn, &probe_args, &estimator)); + OK(sdis_scene_get_boundary_position + (box_scn, probe_args.iprim, probe_args.uv, pos)); printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos)); check_estimator(estimator, N, ref); - BA(GREEN(NULL, N, iprim, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, 0, iprim, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, 12, uv, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, NULL, F, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, -1, 1, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, F, 0, 0, 0, &green)); - BA(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, NULL)); + BA(GREEN(NULL, &probe_args, &green)); + BA(GREEN(box_scn, NULL, &green)); + BA(GREEN(box_scn, &probe_args, NULL)); + probe_args.nrealisations = 0; + BA(GREEN(box_scn, &probe_args, &green)); + probe_args.nrealisations = N; + probe_args.iprim = 12; + BA(GREEN(box_scn, &probe_args, &green)); + probe_args.iprim = 6; + probe_args.side = SDIS_SIDE_NULL__; + BA(GREEN(box_scn, &probe_args, &green)); + probe_args.side = SDIS_FRONT; + OK(GREEN(box_scn, &probe_args, &green)); - OK(GREEN(box_scn, N, iprim, uv, F, 1, 0, 0, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_green_function_ref_put(green)); @@ -339,27 +354,29 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); /* Dump paths */ - OK(SOLVE(box_scn, N_dump, iprim, uv, time_range, F, 1.0, 0, 0, - SDIS_HEAT_PATH_ALL, &estimator)); + probe_args.nrealisations = N_dump; + probe_args.register_paths = SDIS_HEAT_PATH_ALL; + OK(SOLVE(box_scn, &probe_args, &estimator)); dump_heat_paths(fp, estimator); OK(sdis_estimator_ref_put(estimator)); /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = UNKNOWN_TEMPERATURE; - BA(SOLVE(box_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); + BA(SOLVE(box_scn, &probe_args, &estimator)); fluid_param->temperature = Tf; - uv[0] = 0.5; - iprim = 3; - BA(SOLVE(square_scn, N, 4, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); - OK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos)); - printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos)); - check_estimator(estimator, N, ref); + probe_args.nrealisations = N; + probe_args.register_paths = SDIS_HEAT_PATH_NONE; + probe_args.uv[0] = 0.5; + probe_args.iprim = 4; + + BA(SOLVE(square_scn, &probe_args, &estimator)); + probe_args.iprim = 3; + OK(SOLVE(square_scn, &probe_args, &estimator)); - OK(GREEN(square_scn, N, iprim, uv, F, 1, 0, 0, &green)); + OK(GREEN(square_scn, &probe_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_estimator_ref_put(estimator)); @@ -368,7 +385,7 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = UNKNOWN_TEMPERATURE; - BA(SOLVE(square_scn, N, iprim, uv, time_range, F, 1.0, 0, 0, 0, &estimator)); + BA(SOLVE(square_scn, &probe_args, &estimator)); fluid_param->temperature = Tf; #undef F @@ -380,40 +397,78 @@ main(int argc, char** argv) sides[2] = SDIS_FRONT; sides[3] = SDIS_FRONT; + bound_args.nrealisations = N; + bound_args.sides = sides; + bound_args.primitives = prims; + bound_args.nprimitives = 2; + bound_args.time_range[0] = INF; + bound_args.time_range[1] = INF; + #define SOLVE sdis_solve_boundary #define GREEN sdis_solve_boundary_green_function prims[0] = 6; prims[1] = 7; - BA(SOLVE(NULL, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, 0, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, NULL, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, prims, NULL, 2, time_range, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, prims, sides, 0, time_range, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, prims, sides, 2, NULL, 1.0, 0, 0, 0, &estimator)); - BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL)); - tr[0] = 1; - BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL)); - tr[1] = 0; - BA(SOLVE(box_scn, N, prims, sides, 2, tr, 1.0, 0, 0, 0, NULL)); + + BA(SOLVE(NULL, &bound_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &bound_args, NULL)); + bound_args.nrealisations = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.nrealisations = N; + bound_args.primitives = NULL; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.primitives = prims; + bound_args.sides = NULL; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.sides = sides; + bound_args.nprimitives = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.nprimitives = 2; + prims[0] = 12; + BA(SOLVE(box_scn, &bound_args, &estimator)); + prims[0] = 6; + sides[0] = SDIS_SIDE_NULL__; + BA(SOLVE(box_scn, &bound_args, &estimator)); + sides[0] = SDIS_FRONT; + bound_args.time_range[0] = bound_args.time_range[1] = -1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = 1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[1] = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = bound_args.time_range[1] = INF; /* Average temperature on the right side of the box */ - OK(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); + OK(SOLVE(box_scn, &bound_args, &estimator)); printf("Average temperature of the right side of the box = "); check_estimator(estimator, N, ref); - BA(GREEN(NULL, N, prims, sides, 2, 1.0, 0, 0, &green)); - BA(GREEN(box_scn, 0, prims, sides, 2, 1.0, 0, 0, &green)); - BA(GREEN(box_scn, N, NULL, sides, 2, 1.0, 0, 0, &green)); - BA(GREEN(box_scn, N, prims, NULL, 2, 1.0, 0, 0, &green)); - BA(GREEN(box_scn, N, prims, sides, 0, 1.0, 0, 0, &green)); - BA(GREEN(box_scn, N, prims, sides, 2, 0.0, 0, 0, &green)); - BA(GREEN(box_scn, N, prims, sides, 2, 1.0, 0, 0, NULL)); + BA(GREEN(NULL, &bound_args, &green)); + BA(GREEN(box_scn, NULL, &green)); + BA(GREEN(box_scn, &bound_args, NULL)); + bound_args.nrealisations = 0; + BA(GREEN(box_scn, &bound_args, &green)); + bound_args.nrealisations = N; + bound_args.primitives = NULL; + BA(GREEN(box_scn, &bound_args, &green)); + bound_args.primitives = prims; + bound_args.sides = NULL; + BA(GREEN(box_scn, &bound_args, &green)); + bound_args.sides = sides; + bound_args.nprimitives = 0; + BA(GREEN(box_scn, &bound_args, &green)); + bound_args.nprimitives = 2; + prims[0] = 12; + BA(GREEN(box_scn, &bound_args, &green)); + prims[0] = 6; + sides[0] = SDIS_SIDE_NULL__; + BA(GREEN(box_scn, &bound_args, &green)); + sides[0] = SDIS_FRONT; + - OK(GREEN(box_scn, N, prims, sides, 2, 1.0, 0, 0, &green)); + OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_green_function_ref_put(green)); @@ -421,21 +476,35 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); /* Dump path */ - OK(SOLVE(box_scn, N_dump, prims, sides, 2, time_range, 1.0, 0, 0, - SDIS_HEAT_PATH_ALL, &estimator)); + bound_args.nrealisations = N_dump; + bound_args.register_paths = SDIS_HEAT_PATH_ALL; + + /* Check simulation error handling when paths are registered */ + fluid_param->temperature = UNKNOWN_TEMPERATURE; + BA(SOLVE(box_scn, &bound_args, &estimator)); + + /* Dump path */ + fluid_param->temperature = Tf; + OK(SOLVE(box_scn, &bound_args, &estimator)); dump_heat_paths(fp, estimator); OK(sdis_estimator_ref_put(estimator)); + /* Switch in 2D */ + bound_args.nrealisations = N; + bound_args.register_paths = SDIS_HEAT_PATH_NONE; + bound_args.nprimitives = 1; + prims[0] = 4; + BA(SOLVE(square_scn, &bound_args, &estimator)); + /* Average temperature on the right side of the square */ prims[0] = 3; - sides[0] = SDIS_FRONT; - OK(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, 0, &estimator)); + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average temperature of the right side of the square = "); check_estimator(estimator, N, ref); - OK(GREEN(square_scn, N, prims, sides, 1, 1.0, 0, 0, &green)); + OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_green_function_ref_put(green)); @@ -443,16 +512,14 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); /* Dump path */ - OK(SOLVE(square_scn, N_dump, prims, sides, 1, time_range, 1.0, 0, 0, - SDIS_HEAT_PATH_ALL, &estimator)); + bound_args.nrealisations = N_dump; + bound_args.register_paths = SDIS_HEAT_PATH_ALL; + OK(SOLVE(square_scn, &bound_args, &estimator)); dump_heat_paths(fp, estimator); OK(sdis_estimator_ref_put(estimator)); - /* Check out of bound prims */ - prims[0] = 12; - BA(SOLVE(box_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); - prims[0] = 4; - BA(SOLVE(square_scn, N, prims, sides, 1, time_range, 1.0, 0, 0, 0, &estimator)); + bound_args.register_paths = SDIS_HEAT_PATH_NONE; + bound_args.nrealisations = N; /* Average temperature on the left+right sides of the box */ prims[0] = 2; @@ -462,13 +529,14 @@ main(int argc, char** argv) ref = (ref + Tb) / 2; - OK(SOLVE(box_scn, N, prims, sides, 4, time_range, 1.0, 0, 0, 0, &estimator)); + bound_args.nprimitives = 4; + OK(SOLVE(box_scn, &bound_args, &estimator)); printf("Average temperature of the left+right sides of the box = "); check_estimator(estimator, N, ref); - OK(GREEN(box_scn, N, prims, sides, 4, 1.0, 0, 0, &green)); + OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_green_function_ref_put(green)); @@ -478,13 +546,14 @@ main(int argc, char** argv) /* Average temperature on the left+right sides of the square */ prims[0] = 1; prims[1] = 3; - OK(SOLVE(square_scn, N, prims, sides, 2, time_range, 1.0, 0, 0, 0, &estimator)); + bound_args.nprimitives = 2; + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average temperature of the left+right sides of the square = "); check_estimator(estimator, N, ref); - OK(GREEN(square_scn, N, prims, sides, 2, 1.0, 0, 0, &green)); + OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); check_estimator(estimator2, N, ref); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c @@ -234,16 +234,16 @@ main(int argc, char** argv) struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; + struct sdis_solve_probe_boundary_flux_args probe_args = + SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; + struct sdis_solve_boundary_flux_args bound_args = + SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT; struct interf* interf_props = NULL; struct fluid* fluid_param; enum sdis_estimator_type type; - double uv[2]; double pos[3]; - double time_range[2] = { INF, INF }; - double tr[2]; double analyticT, analyticCF, analyticRF, analyticTF; size_t prims[2]; - size_t iprim; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -345,81 +345,109 @@ main(int argc, char** argv) analyticTF = analyticCF + analyticRF; #define SOLVE sdis_solve_probe_boundary_flux - uv[0] = 0.3; - uv[1] = 0.3; - iprim = 6; - - BA(SOLVE(NULL, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, 0, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, 12, uv, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, NULL, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, NULL, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[0] = 1; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[1] = 0; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - tr[1] = INF; - BA(SOLVE(box_scn, N, iprim, uv, tr, 1.0, Trad, Tref, &estimator)); - - OK(SOLVE(box_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); + probe_args.nrealisations = N; + probe_args.iprim = 6; + probe_args.uv[0] = 0.3; + probe_args.uv[1] = 0.3; + probe_args.time_range[0] = INF; + probe_args.time_range[1] = INF; + probe_args.ambient_radiative_temperature = Trad; + probe_args.reference_temperature = Tref; + BA(SOLVE(NULL, &probe_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &probe_args, NULL)); + probe_args.nrealisations = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.nrealisations = N; + probe_args.iprim = 12; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.iprim = 6; + probe_args.uv[0] = probe_args.uv[1] = 1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.uv[0] = probe_args.uv[1] = 0.3; + probe_args.time_range[0] = probe_args.time_range[1] = -1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = 1; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[1] = 0; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[1] = INF; + BA(SOLVE(box_scn, &probe_args, &estimator)); + probe_args.time_range[0] = INF; + OK(SOLVE(box_scn, &probe_args, &estimator)); + OK(sdis_estimator_get_type(estimator, &type)); CHK(type == SDIS_ESTIMATOR_FLUX); - OK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos)); + OK(sdis_scene_get_boundary_position + (box_scn, probe_args.iprim, probe_args.uv, pos)); printf("Boundary values of the box at (%g %g %g) = ", SPLIT3(pos)); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); - uv[0] = 0.5; - iprim = 3; - BA(SOLVE(square_scn, N, 4, uv, time_range, 1.0, Trad, Tref, &estimator)); - OK(SOLVE(square_scn, N, iprim, uv, time_range, 1.0, Trad, Tref, &estimator)); - OK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos)); + probe_args.uv[0] = 0.5; + probe_args.iprim = 4; + BA(SOLVE(square_scn, &probe_args, &estimator)); + probe_args.iprim = 3; + OK(SOLVE(square_scn, &probe_args, &estimator)); + OK(sdis_scene_get_boundary_position + (square_scn, probe_args.iprim, probe_args.uv, pos)); printf("Boundary values of the square at (%g %g) = ", SPLIT2(pos)); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); #undef F #undef SOLVE - + #define SOLVE sdis_solve_boundary_flux prims[0] = 6; prims[1] = 7; - BA(SOLVE(NULL, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, 0, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, NULL, 2, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 0, time_range, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 2, NULL, 1.0, Trad, Tref, &estimator)); - BA(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, NULL)); - tr[0] = tr[1] = -1; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); - tr[0] = 1; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); - tr[1] = 0; - BA(SOLVE(box_scn, N, prims, 2, tr, 1.0, Trad, Tref, NULL)); + bound_args.nrealisations = N; + bound_args.primitives = prims; + bound_args.nprimitives = 2; + bound_args.time_range[0] = INF; + bound_args.time_range[1] = INF; + bound_args.ambient_radiative_temperature = Trad; + bound_args.reference_temperature = Tref; + + BA(SOLVE(NULL, &bound_args, &estimator)); + BA(SOLVE(box_scn, NULL, &estimator)); + BA(SOLVE(box_scn, &bound_args, NULL)); + bound_args.primitives = NULL; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.primitives = prims; + bound_args.nprimitives = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.nprimitives = 2; + bound_args.time_range[0] = bound_args.time_range[1] = -1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = 1; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[1] = 0; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[1] = INF; + BA(SOLVE(box_scn, &bound_args, &estimator)); + bound_args.time_range[0] = INF; + prims[0] = 12; + BA(SOLVE(box_scn, &bound_args, &estimator)); + prims[0] = 6; + OK(SOLVE(box_scn, &bound_args, &estimator)); /* Average temperature on the right side of the box */ - OK(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); printf("Average values of the right side of the box = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); /* Average temperature on the right side of the square */ + prims[0] = 4; + bound_args.nprimitives = 1; + BA(SOLVE(square_scn, &bound_args, &estimator)); prims[0] = 3; - OK(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average values of the right side of the square = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); - - /* Check out of bound prims */ - prims[0] = 12; - BA(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); - prims[0] = 4; - BA(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); - + /* Average temperature on the left side of the box */ prims[0] = 2; prims[1] = 3; @@ -429,14 +457,16 @@ main(int argc, char** argv) analyticRF = Hrad * (analyticT - Trad); analyticTF = analyticCF + analyticRF; - OK(SOLVE(box_scn, N, prims, 2, time_range, 1.0, Trad, Tref, &estimator)); + bound_args.nprimitives = 2; + OK(SOLVE(box_scn, &bound_args, &estimator)); printf("Average values of the left side of the box = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); /* Average temperature on the left/right side of the square */ prims[0] = 1; - OK(SOLVE(square_scn, N, prims, 1, time_range, 1.0, Trad, Tref, &estimator)); + bound_args.nprimitives = 1; + OK(SOLVE(square_scn, &bound_args, &estimator)); printf("Average values of the left side of the square = "); check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -30,12 +30,12 @@ #define UNKOWN_TEMPERATURE -1 #define IMG_WIDTH 640 #define IMG_HEIGHT 480 -#define SPP 1024 /* #Samples per pixel, i.e. #realisations per pixel */ +#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */ /* * The scene is composed of a solid cube whose temperature is unknown. The * emissivity of the cube is 1 and its convection coefficient with the - * surrounding fluid is null. At the center of the cube there is a spherical + * surrounding fluid at 300K is 0.1. At the center of the cube there is a spherical * fluid cavity whose temperature is 350K. The convection coefficient between * the solid and the cavity is 1 and the emissivity of this interface is null. * The ambient radiative temperature of the system is 300K. @@ -540,16 +540,18 @@ main(int argc, char** argv) struct sdis_interface* interf0 = NULL; struct sdis_interface* interf1 = NULL; struct sdis_scene* scn = NULL; + struct sdis_solve_camera_args solve_args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT; + struct ssp_rng* rng_state = NULL; struct fluid fluid_param = FLUID_NULL; struct solid solid_param = SOLID_NULL; struct interf interface_param = INTERF_NULL; + struct fluid* pfluid_param = NULL; size_t ntris, npos; size_t nreals, nfails; size_t definition[2]; double pos[3]; double tgt[3]; double up[3]; - double trange[2] = {INF, INF}; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -562,7 +564,7 @@ main(int argc, char** argv) create_fluid(dev, &fluid_param, &fluid0); /* Create the fluid1 */ - fluid_param.temperature = UNKOWN_TEMPERATURE; + fluid_param.temperature = 300; fluid_param.rho = 0; fluid_param.cp = 0; create_fluid(dev, &fluid_param, &fluid1); @@ -583,17 +585,12 @@ main(int argc, char** argv) create_interface(dev, solid, fluid0, &interface_param, &interf0); /* Create the fluid1/solid interface */ - interface_param.hc = 0; + interface_param.hc = 0.1; interface_param.epsilon = 1; interface_param.specular_fraction = 1; interface_param.temperature = UNKOWN_TEMPERATURE; create_interface(dev, fluid1, solid, &interface_param, &interf1); - /* Release the ownership onto the media */ - OK(sdis_medium_ref_put(solid)); - OK(sdis_medium_ref_put(fluid0)); - OK(sdis_medium_ref_put(fluid1)); - /* Setup the cube geometry */ OK(s3dut_create_cuboid(&allocator, 2, 2, 2, &msh)); OK(s3dut_mesh_get_data(msh, &msh_data)); @@ -633,36 +630,40 @@ main(int argc, char** argv) dump_mesh(stdout, geom.positions, npos, geom.indices, ntris); exit(0); #endif - BA(sdis_solve_camera(NULL, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, NULL, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, NULL, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 0, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, -1, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, 0, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, 0, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - 0, SDIS_HEAT_PATH_NONE, &buf)); - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, NULL)); - - trange[0] = -1; - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - trange[0] = 10; trange[1] = 1; - BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); - trange[0] = trange[1] = INF; + solve_args.cam = cam; + solve_args.time_range[0] = INF; + solve_args.time_range[0] = INF; + solve_args.image_resolution[0] = IMG_WIDTH; + solve_args.image_resolution[1] = IMG_HEIGHT; + solve_args.ambient_radiative_temperature = 300; + solve_args.reference_temperature = 300; + solve_args.spp = SPP; + + BA(sdis_solve_camera(NULL, &solve_args, &buf)); + BA(sdis_solve_camera(scn, NULL, &buf)); + BA(sdis_solve_camera(scn, &solve_args, NULL)); + solve_args.cam = NULL; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.cam = cam; + solve_args.fp_to_meter = 0; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.fp_to_meter = 1; + solve_args.ambient_radiative_temperature = -1; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.ambient_radiative_temperature = 300; + solve_args.reference_temperature = -1; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.reference_temperature = 300; + solve_args.time_range[0] = solve_args.time_range[1] = -1; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.time_range[0] = 1; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.time_range[1] = 0; + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.time_range[0] = solve_args.time_range[1] = INF; /* Launch the simulation */ - OK(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, - SPP, SDIS_HEAT_PATH_NONE, &buf)); + OK(sdis_solve_camera(scn, &solve_args, &buf)); BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals)); BA(sdis_estimator_buffer_get_realisation_count(buf, NULL)); @@ -680,6 +681,10 @@ main(int argc, char** argv) BA(sdis_estimator_buffer_get_realisation_time(buf, NULL)); OK(sdis_estimator_buffer_get_realisation_time(buf, &time)); + BA(sdis_estimator_buffer_get_rng_state(NULL, &rng_state)); + BA(sdis_estimator_buffer_get_rng_state(buf, NULL)); + OK(sdis_estimator_buffer_get_rng_state(buf, &rng_state)); + CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP); fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE); @@ -695,9 +700,20 @@ main(int argc, char** argv) /* Write the image */ dump_image(buf); + OK(sdis_estimator_buffer_ref_put(buf)); + + pfluid_param = sdis_data_get(sdis_medium_get_data(fluid1)); + pfluid_param->temperature = UNKOWN_TEMPERATURE; + + /* Check simulation error handling */ + BA(sdis_solve_camera(scn, &solve_args, &buf)); + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + BA(sdis_solve_camera(scn, &solve_args, &buf)); /* Release memory */ - OK(sdis_estimator_buffer_ref_put(buf)); + OK(sdis_medium_ref_put(solid)); + OK(sdis_medium_ref_put(fluid0)); + OK(sdis_medium_ref_put(fluid1)); OK(sdis_scene_ref_put(scn)); OK(sdis_camera_ref_put(cam)); OK(sdis_interface_ref_put(interf0)); diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -223,8 +223,8 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_solve_medium_args solve_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; struct context ctx; - const double trange[2] = {INF, INF}; double ref; double v, v0, v1; size_t nreals; @@ -352,13 +352,31 @@ main(int argc, char** argv) OK(sdis_scene_get_medium_spread(scn, fluid1, &v)); CHK(v == 0); - BA(sdis_solve_medium(NULL, N, solid0, trange, 1.f, -1, 0, 0, &estimator)); - BA(sdis_solve_medium(scn, 0, solid0, trange, 1.f, -1, 0, 0, &estimator)); - BA(sdis_solve_medium(scn, N, NULL, trange, 1.f, -1, 0, 0, &estimator)); - BA(sdis_solve_medium(scn, N, solid0, NULL, 1.f, -1, 0, 0, &estimator)); - BA(sdis_solve_medium(scn, N, solid0, trange, 0.f, -1, 0, 0, &estimator)); - BA(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, NULL)); - OK(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.medium = solid0; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + BA(sdis_solve_medium(NULL, &solve_args, &estimator)); + BA(sdis_solve_medium(scn, NULL, &estimator)); + BA(sdis_solve_medium(scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.nrealisations = N; + solve_args.medium = NULL; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.medium = solid0; + solve_args.fp_to_meter = 0; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.fp_to_meter = 1; + solve_args.time_range[0] = solve_args.time_range[1] = -1; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.time_range[0] = 1; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.time_range[1] = 0; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.time_range[0] = solve_args.time_range[1] = INF; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); @@ -371,7 +389,20 @@ main(int argc, char** argv) CHK(nreals + nfails == N); OK(sdis_estimator_ref_put(estimator)); - OK(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.medium = solid1; + + /* Check simulation error handling when paths are registered */ + solve_args.nrealisations = 10; + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + fluid_param->temperature = -1; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + fluid_param->temperature = Tf1; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); + OK(sdis_estimator_ref_put(estimator)); + solve_args.nrealisations = N; + solve_args.register_paths = SDIS_HEAT_PATH_NONE; + + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -383,13 +414,6 @@ main(int argc, char** argv) CHK(nreals + nfails == N); OK(sdis_estimator_ref_put(estimator)); -#if 0 - OK(sdis_solve_medium(scn, 1, solid1, trange, 1.f, -1, 0, - SDIS_HEAT_PATH_ALL, &estimator)); - dump_heat_paths(stderr, estimator); - OK(sdis_estimator_ref_put(estimator)); -#endif - /* Create a new scene with the same medium in the 2 super shapes */ OK(sdis_scene_ref_put(scn)); ctx.interf0 = solid0_fluid0; @@ -400,8 +424,11 @@ main(int argc, char** argv) OK(sdis_scene_get_medium_spread(scn, solid0, &v)); CHK(eq_eps(v, v0+v1, 1.e-6)); - BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator)); - OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.medium = solid1; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.medium = solid0; + solve_args.nrealisations = Np; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -413,16 +440,23 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE*3)); /* Solve green */ - BA(sdis_solve_medium_green_function(NULL, Np, solid0, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, 0, solid0, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, NULL, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 0.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, -1, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, NULL)); - BA(sdis_solve_medium_green_function(scn, Np, solid1, 1.0, 0, 0, &green)); - OK(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, &green)); - - OK(sdis_green_function_solve(green, trange, &estimator2)); + BA(sdis_solve_medium_green_function(NULL, &solve_args, &green)); + BA(sdis_solve_medium_green_function(scn, NULL, &green)); + BA(sdis_solve_medium_green_function(scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(sdis_solve_medium_green_function(scn, &solve_args, &green)); + solve_args.nrealisations = Np; + solve_args.medium = NULL; + BA(sdis_solve_medium_green_function(scn, &solve_args, &green)); + solve_args.medium = solid1; + BA(sdis_solve_medium_green_function(scn, &solve_args, &green)); + solve_args.medium = solid0; + solve_args.fp_to_meter = 0; + BA(sdis_solve_medium_green_function(scn, &solve_args, &green)); + solve_args.fp_to_meter = 1; + OK(sdis_solve_medium_green_function(scn, &solve_args, &green)); + + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -208,8 +208,8 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_solve_medium_args solve_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; struct context ctx; - const double trange[2] = {INF, INF}; double a, a0, a1; double ref; double* positions = NULL; @@ -334,8 +334,13 @@ main(int argc, char** argv) /* Rough estimation since the disk is coarsely discretized */ CHK(eq_eps(a1, PI, 1.e-1)); + solve_args.nrealisations = N; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + /* Estimate the temperature of the square */ - OK(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.medium = solid0; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -348,7 +353,8 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator)); /* Estimate the temperature of the disk */ - OK(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.medium = solid1; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -371,8 +377,11 @@ main(int argc, char** argv) CHK(eq_eps(a, a0+a1, 1.e-6)); /* Estimate the temperature of the square and disk shapes */ - BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator)); - OK(sdis_solve_medium(scn, Np, solid0, trange, 1.f, -1, 0, 0, &estimator)); + solve_args.medium = solid1; + solve_args.nrealisations = Np; + BA(sdis_solve_medium(scn, &solve_args, &estimator)); + solve_args.medium = solid0; + OK(sdis_solve_medium(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -385,16 +394,12 @@ main(int argc, char** argv) CHK(nreals + nfails == Np); /* Solve green */ - BA(sdis_solve_medium_green_function(NULL, Np, solid0, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, 0, solid0, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, NULL, 1.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 0.0, 0, 0, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, -1, &green)); - BA(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, NULL)); - BA(sdis_solve_medium_green_function(scn, Np, solid1, 1.0, 0, 0, &green)); - OK(sdis_solve_medium_green_function(scn, Np, solid0, 1.0, 0, 0, &green)); - - OK(sdis_green_function_solve(green, trange, &estimator2)); + BA(sdis_solve_medium_green_function(NULL, &solve_args, &green)); + BA(sdis_solve_medium_green_function(scn, NULL, &green)); + BA(sdis_solve_medium_green_function(scn, &solve_args, NULL)); + OK(sdis_solve_medium_green_function(scn, &solve_args, &green)); + + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -16,6 +16,7 @@ #include "sdis.h" #include "test_sdis_utils.h" +#include <star/ssp.h> #include <rsys/math.h> /* @@ -263,14 +264,14 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct dump_path_context dump_ctx = DUMP_PATH_CONTEXT_NULL; struct context ctx; struct fluid* fluid_param; struct solid* solid_param; struct interf* interface_param; + struct ssp_rng* rng_state = NULL; enum sdis_estimator_type type; - double pos[3]; - double time_range[2]; double ref; const size_t N = 1000; const size_t N_dump = 10; @@ -338,17 +339,30 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf)); /* Test the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - pos[2] = 0.5; - time_range[0] = time_range[1] = INF; - BA(sdis_solve_probe(NULL, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, 0, pos, time_range, 1.0, 0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, NULL, time_range, 1.0, 0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, 0, 0, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 0, 0, -1, 0, &estimator)); - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, NULL)); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.position[2] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + BA(sdis_solve_probe(NULL, &solve_args, &estimator)); + BA(sdis_solve_probe(scn, NULL, &estimator)); + BA(sdis_solve_probe(scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + solve_args.nrealisations = N; + solve_args.fp_to_meter = 0; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + solve_args.fp_to_meter = 1; + solve_args.time_range[0] = solve_args.time_range[1] = -1; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + solve_args.time_range[0] = 1; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + solve_args.time_range[1] = 0; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + solve_args.time_range[0] = solve_args.time_range[1] = INF; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); BA(sdis_estimator_get_type(estimator, NULL)); BA(sdis_estimator_get_type(NULL, &type)); @@ -384,9 +398,13 @@ main(int argc, char** argv) BA(sdis_estimator_get_realisation_time(NULL, &time)); OK(sdis_estimator_get_realisation_time(estimator, &time)); + BA(sdis_estimator_get_rng_state(NULL, &rng_state)); + BA(sdis_estimator_get_rng_state(estimator, NULL)); + OK(sdis_estimator_get_rng_state(estimator, &rng_state)); + ref = 300; printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + SPLIT3(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -402,23 +420,26 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, &solve_args, &estimator)); fluid_param->temperature = 300; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); - - BA(sdis_solve_probe_green_function(NULL, N, pos, 1.0, 0, 0, &green)); - BA(sdis_solve_probe_green_function(scn, 0, pos, 1.0, 0, 0, &green)); - BA(sdis_solve_probe_green_function(scn, N, NULL, 1.0, 0, 0, &green)); - BA(sdis_solve_probe_green_function(scn, N, pos, 0.0, 0, 0, &green)); - BA(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, -1, &green)); - BA(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, NULL)); - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); - - BA(sdis_green_function_solve(NULL, time_range, &estimator2)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + + BA(sdis_solve_probe_green_function(NULL, &solve_args, &green)); + BA(sdis_solve_probe_green_function(scn, NULL, &green)); + BA(sdis_solve_probe_green_function(scn, &solve_args, NULL)); + solve_args.nrealisations = 0; + BA(sdis_solve_probe_green_function(scn, &solve_args, &green)); + solve_args.nrealisations = N; + solve_args.fp_to_meter = 0; + BA(sdis_solve_probe_green_function(scn, &solve_args, &green)); + solve_args.fp_to_meter = 1; + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + + BA(sdis_green_function_solve(NULL, solve_args.time_range, &estimator2)); BA(sdis_green_function_solve(green, NULL, &estimator2)); - BA(sdis_green_function_solve(green, time_range, NULL)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + BA(sdis_green_function_solve(green, solve_args.time_range, NULL)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -432,15 +453,22 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); BA(sdis_estimator_get_paths_count(NULL, &n)); BA(sdis_estimator_get_paths_count(estimator, NULL)); OK(sdis_estimator_get_paths_count(estimator, &n)); CHK(n == 0); OK(sdis_estimator_ref_put(estimator)); - OK(sdis_solve_probe(scn, N_dump, pos, time_range, 1.0, 0, 0, - SDIS_HEAT_PATH_ALL, &estimator)); + solve_args.nrealisations = N_dump; + solve_args.register_paths = SDIS_HEAT_PATH_ALL; + + /* Check simulation error handling when paths are registered */ + fluid_param->temperature = -1; + BA(sdis_solve_probe(scn, &solve_args, &estimator)); + + fluid_param->temperature = 300; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_paths_count(estimator, &n)); CHK(n == N_dump); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -162,10 +162,9 @@ main(int argc, char** argv) struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; struct sdis_interface* interfaces[12]; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct context ctx; struct interf* interface_param = NULL; - double pos[3]; - double time_range[2] = { INF, INF }; double ref; const size_t N = 10000; size_t nreals; @@ -239,20 +238,23 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(T350)); /* Launch the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - pos[2] = 0.5; - time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.position[2] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); /* Print the estimation results */ - ref = 350 * pos[2] + (1-pos[2]) * 300; + ref = 350 * solve_args.position[2] + (1-solve_args.position[2]) * 300; printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + SPLIT3(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -262,8 +264,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE)); /* Check green */ - OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -159,10 +159,9 @@ main(int argc, char** argv) struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; struct sdis_interface* interfaces[4]; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct context ctx; struct interf* interface_param = NULL; - double pos[2]; - double time_range[2]; double ref; const size_t N = 10000; size_t nreals; @@ -186,14 +185,14 @@ main(int argc, char** argv) /* Create the fluid/solid interface with no limit conidition */ interface_shader.convection_coef = null_interface_value; - interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; - interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; OK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &Tnone)); interface_shader.convection_coef = null_interface_value; - interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; - interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; interface_shader.front.temperature = interface_get_temperature; /* Create the fluid/solid interface with a fixed temperature of 300K */ @@ -238,20 +237,23 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(T350)); /* Launch the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); - /* Print the estimation results */ - ref = 350 * pos[0] + (1-pos[0]) * 300; + ref = 350 * solve_args.position[0] + (1-solve_args.position[0]) * 300; printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + SPLIT2(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -261,8 +263,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE*2)); /* Check green function */ - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -184,12 +184,11 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct s3dut_mesh* msh = NULL; struct s3dut_mesh_data msh_data; struct context ctx = CONTEXT_NULL; struct interf* interface_param = NULL; - double pos[3]; - double time_range[2]; double ref; const size_t N = 10000; size_t ntris; @@ -295,20 +294,23 @@ main(int argc, char** argv) sa_release(ctx.indices); /* Launch the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - pos[2] = 0.5; - time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe( scn, N, pos, time_range, 1.0, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.position[2] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); /* Print the estimation results */ - ref = 350 * pos[2] + (1-pos[2]) * 300; + ref = 350 * solve_args.position[2] + (1-solve_args.position[2]) * 300; printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); + SPLIT3(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -318,8 +320,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE)); /* Check green function */ - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -181,10 +181,9 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct context ctx = CONTEXT_NULL; struct interf* interface_param = NULL; - double pos[2]; - double time_range[2]; double ref; const size_t N = 10000; size_t nsegs; @@ -289,19 +288,21 @@ main(int argc, char** argv) sa_release(ctx.indices); /* Launch the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_time(estimator, &time)); /* Print the estimation results */ - ref = 350 * pos[0] + (1-pos[0]) * 300; + ref = 350 * solve_args.position[0] + (1-solve_args.position[0]) * 300; printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + SPLIT2(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -311,8 +312,8 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, 3*T.SE)); /* Check green function */ - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -148,10 +148,9 @@ main(int argc, char** argv) struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct context ctx; struct fluid* fluid_param; - double pos[2]; - double time_range[2] = { INF, INF }; double ref; const size_t N = 1000; size_t nreals; @@ -199,10 +198,13 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf)); /* Test the solver */ - pos[0] = 0.5; - pos[1] = 0.5; - time_range[0] = time_range[1] = INF; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = 0.5; + solve_args.position[1] = 0.5; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); @@ -210,7 +212,7 @@ main(int argc, char** argv) ref = 300; printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); + SPLIT2(solve_args.position), ref, T.E, T.SE); printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); @@ -218,8 +220,8 @@ main(int argc, char** argv) CHK(nfails < N/1000); CHK(eq_eps(T.E, ref, T.SE)); - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); - OK(sdis_green_function_solve(green, time_range, &estimator2)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -230,7 +232,7 @@ main(int argc, char** argv) /* The external fluid cannot have an unknown temperature */ fluid_param->temperature = -1; - BA(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + BA(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -156,6 +156,7 @@ solve(struct sdis_scene* scn, const double pos[]) struct sdis_estimator* estimator; struct sdis_estimator* estimator2; struct sdis_green_function* green; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct sdis_mc T; size_t nreals; size_t nfails; @@ -168,8 +169,17 @@ solve(struct sdis_scene* scn, const double pos[]) x = pos[0] - 0.5; ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.position[0] = pos[0]; + solve_args.position[1] = pos[1]; + solve_args.position[2] = dim == SDIS_SCENE_2D ? 0 : pos[2]; + time_current(&t0); - OK(sdis_solve_probe(scn, N, pos, time_range, 1.0, 0, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); @@ -177,8 +187,6 @@ solve(struct sdis_scene* scn, const double pos[]) OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_scene_get_dimension(scn, &dim)); - switch(dim) { case SDIS_SCENE_2D: printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", @@ -198,7 +206,7 @@ solve(struct sdis_scene* scn, const double pos[]) CHK(eq_eps(T.E, ref, T.SE*3)); time_current(&t0); - OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, 0, 0, &green)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); time_current(&t1); OK(sdis_green_function_solve(green, time_range, &estimator2)); time_current(&t2); diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -229,20 +229,24 @@ static void check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) { struct sdis_estimator* estimator = NULL; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct sdis_mc T = SDIS_MC_NULL; size_t nreals; size_t nfails; double pos[3] = {0,0}; - double time_range[2] = { INF, INF }; size_t i; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + solve_args.nrealisations = N; + FOR_EACH(i, 0, nrefs) { double Tc; - pos[0] = refs[i].pos[0]; - pos[1] = refs[i].pos[1]; - pos[2] = refs[i].pos[2]; + solve_args.position[0] = refs[i].pos[0]; + solve_args.position[1] = refs[i].pos[1]; + solve_args.position[2] = refs[i].pos[2]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -250,18 +250,22 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) { struct sdis_estimator* estimator = NULL; struct sdis_mc T = SDIS_MC_NULL; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals; size_t nfails; double pos[2] = {0,0}; - double time_range[2] = { INF, INF }; size_t i; + solve_args.nrealisations = N; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + FOR_EACH(i, 0, nrefs) { double Tc; - pos[0] = refs[i].pos[0]; - pos[1] = refs[i].pos[1]; + solve_args.position[0] = refs[i].pos[0]; + solve_args.position[1] = refs[i].pos[1]; - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -261,6 +261,7 @@ main(int argc, char** argv) struct sdis_interface* interf_solid1_fluid = NULL; struct sdis_interface* interf_solid2_fluid = NULL; struct sdis_interface* interfaces[10/*#segment*/]; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct sdis_mc T = SDIS_MC_NULL; double Tref; double time_range[2] = { INF, INF }; @@ -441,7 +442,12 @@ main(int argc, char** argv) FATAL("Unreachable code.\n"); } - OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); + solve_args.nrealisations = N; + solve_args.position[0] = pos[0]; + solve_args.position[1] = pos[1]; + solve_args.time_range[0] = time_range[0]; + solve_args.time_range[1] = time_range[1]; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &T)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c @@ -225,9 +225,9 @@ main(int argc, char** argv) struct sdis_interface* interf_solid_fluid2 = NULL; struct sdis_interface* interfaces[12/*#max primitives*/]; struct sdis_mc T = SDIS_MC_NULL; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals, nfails; double pos[3]; - double time_range[2] = { INF, INF }; double Tref; double x; (void)argc, (void)argv; @@ -353,10 +353,17 @@ main(int argc, char** argv) x = pos[1]; Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA); + solve_args.nrealisations = N; + solve_args.position[0] = pos[0]; + solve_args.position[1] = pos[1]; + solve_args.position[2] = pos[2]; + solve_args.time_range[0] = INF; + solve_args.time_range[1] = INF; + printf(">>> 2D\n"); time_current(&t0); - OK(sdis_solve_probe(scn_2d, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); + OK(sdis_solve_probe(scn_2d, &solve_args, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); printf("Elapsed time = %s\n", dump); @@ -375,7 +382,7 @@ main(int argc, char** argv) printf("\n>>> 3D\n"); time_current(&t0); - OK(sdis_solve_probe(scn_3d, N, pos, time_range, 1.f, -1, 0, 0, &estimator)); + OK(sdis_solve_probe(scn_3d, &solve_args, &estimator)); time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); printf("Elapsed time = %s\n", dump);