stardis-solver

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

commit 814724c9eb2cc8c12148f321d9b5fa3f58f48e95
parent 97b6e62c27c7a9c7d7590612a5c4bb40b2be58a4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 20 Nov 2020 12:27:25 +0100

Merge branch 'release_0.11'

Diffstat:
MREADME.md | 137++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------
Mcmake/CMakeLists.txt | 8++++----
Msrc/sdis.h | 223++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------
Msrc/sdis_Xd_begin.h | 1-
Msrc/sdis_green.c | 239++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
Msrc/sdis_green.h | 7++++++-
Msrc/sdis_heat_path.h | 10----------
Msrc/sdis_heat_path_boundary_Xd.h | 34++++++++++++++--------------------
Msrc/sdis_heat_path_conductive_Xd.h | 11+++++------
Msrc/sdis_heat_path_convective_Xd.h | 5++---
Msrc/sdis_heat_path_radiative_Xd.h | 17+++++------------
Msrc/sdis_misc.h | 10++++------
Msrc/sdis_misc_Xd.h | 19++++++++++++-------
Msrc/sdis_realisation.c | 16+++++++---------
Msrc/sdis_realisation.h | 21---------------------
Msrc/sdis_realisation_Xd.h | 49+++++++++++++++++++++----------------------------
Msrc/sdis_scene.c | 82+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------
Msrc/sdis_scene_Xd.h | 104+++++++++++++++++++++++++++++++++----------------------------------------------
Msrc/sdis_scene_c.h | 2++
Msrc/sdis_solve.c | 21+++++----------------
Msrc/sdis_solve_boundary_Xd.h | 45+++++++++++++++++++--------------------------
Msrc/sdis_solve_medium_Xd.h | 30++++++++++--------------------
Msrc/sdis_solve_probe_Xd.h | 41+++++++++++++++++------------------------
Msrc/sdis_solve_probe_boundary_Xd.h | 45++++++++++++++++++++-------------------------
Dsrc/test_sdis_compute_mean_power.c | 343-------------------------------------------------------------------------------
Asrc/test_sdis_compute_power.c | 345+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_sdis_conducto_radiative.c | 96++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------
Msrc/test_sdis_conducto_radiative_2d.c | 18++++++++++++------
Msrc/test_sdis_convection.c | 29+++++++++++++++++++----------
Msrc/test_sdis_convection_non_uniform.c | 29+++++++++++++++++++----------
Msrc/test_sdis_flux.c | 326++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------
Msrc/test_sdis_scene.c | 155+++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------
Msrc/test_sdis_solid_random_walk_robustness.c | 18+++++++++++-------
Msrc/test_sdis_solve_boundary.c | 87++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------
Msrc/test_sdis_solve_boundary_flux.c | 31++++++++++++++++++++-----------
Msrc/test_sdis_solve_camera.c | 25++++++++++++-------------
Msrc/test_sdis_solve_medium.c | 23+++++++++++------------
Msrc/test_sdis_solve_medium_2d.c | 17+++++++++++------
Msrc/test_sdis_solve_probe.c | 66++++++++++++++++++++++++++++++++++++++++++++++++++----------------
Msrc/test_sdis_solve_probe2.c | 14++++++++++----
Msrc/test_sdis_solve_probe2_2d.c | 14++++++++++----
Msrc/test_sdis_solve_probe3.c | 14++++++++++----
Msrc/test_sdis_solve_probe3_2d.c | 14++++++++++----
Msrc/test_sdis_solve_probe_2d.c | 12+++++++++---
Msrc/test_sdis_transcient.c | 31+++++++++++++++++++++----------
Msrc/test_sdis_utils.c | 52++++++++++++++++++++++++++++++++++++++++------------
Msrc/test_sdis_utils.h | 4++--
Msrc/test_sdis_volumic_power.c | 346+++++++++++++++++++++++++++++++++++++++++++++++++------------------------------
Msrc/test_sdis_volumic_power2.c | 13+++++++++----
Msrc/test_sdis_volumic_power2_2d.c | 16++++++++++------
Msrc/test_sdis_volumic_power3_2d.c | 10++++++++--
Msrc/test_sdis_volumic_power4.c | 19+++++++++++++++----
52 files changed, 2005 insertions(+), 1339 deletions(-)

diff --git a/README.md b/README.md @@ -1,11 +1,97 @@ -# Stardis - -The purpose of this library is to solve coupled convecto - conducto - radiative -thermal problems in 2D and 3D environments. +# Stardis-Solver + +Stardis-Solver is *free software* that solves *coupled* convecto - conducto - +radiative *thermal problems* in *complex* 2D and 3D *environments*. This C89 +library internally relies on *Monte-Carlo* algorithms based on reformulations +of the main heat transfer phenomena as cross-recursive "thermal paths" that +explore space and time until a boundary condition or an initial condition is +found. The key concept here is that heat transfer phenomena are not considered +separately but naturally coupled via cross-recursive [Monte-Carlo +algorithms](https://hal.archives-ouvertes.fr/hal-02419604/). + +The hypothesis these algorithms are based upon are the following: + +- *conduction*: the discretization of thermal heat transfer in solids introduces + the notion of a conductive path length within the Monte-Carlo algorithm. + Solutions obtained using this algorithm are formally exact at the limit of a + null path length. In practice, this path length has to be adapted for a given + geometric configuration so that its value is small compared to the smallest + typical length of a solid. +- *convection*: fluid media are supposed to be isothermal, even if their + temperature may vary with time. This hypothesis relies on the assumption of + perfectly agitated fluids. +- *radiation*: local radiative transfer is linearised, i.e. instead of writing + the spectrally integrated net flux as a difference of temperatures to the + power 4, it is assumed of the same form as the convective flux (as a + difference of temperatures, multiplied by a radiative exchange coefficient). + In order to be valid, this representation of radiative transfer exchanges + requires that the temperature at any position and time is close to a known + reference temperature. + +In Stardis-Solver the system to simulate is represented by a *scene* whose +geometry defines the contour of the object only: in contrast to legacy thermal +solvers *no volumetric mesh* has to be provided. Each geometric primitive as an +associated *interface* that defines its physical properties (e.g. surface +emissivity) and reference the *media* defining the thermal properties on both +side of the primitive. The boundary and initial conditions are defined defined +on the interfaces (convection coefficient, fixed temperature/flux, etc.) and +the media (known temperature). Once fully and consistently described, +computations can be invoked on the resulting scene. + +The main features of the solver are currently: + +- *probe computation*: Stardis-Solver will compute the temperature at any given + position (spatial and temporal). The main idea is that thermal paths start + from this probe position, and scatter in space while going back in time, + until a (spatial) boundary condition or a (temporal) initial condition is + met. In addition to the value of temperature, using a Monte-Carlo method + makes possible to compute a numerical uncertainty (standard deviation of the + weight distribution) over each result. +- *flux computation*: Stardis-Solver can compute the flux over any surface (or + group of surfaces) at any time. Alternatively, it can also compute the total + energy output from a solid element where a internal source of power must be + taken into account. +- *green function*: the value of temperature computed at a probe position is + the average of the Monte-Carlo weight for every thermal path. In practice: + when no internal power sources have to be considered, the weight of any given + thermal path is the temperature of the boundary or initial condition that has + been reached; when internal power sources or imposed fluxes are taken into + account, additional contributions to the weight must be continuously + evaluated by the thermal conduction algorithm, but these contributions are + proportional to the local dissipated power/imposed flux. In any case, the + position and date at the end of each thermal path (and also accumulation + coefficients) can be stored during a first complete Monte-Carlo simulation. + This information, known as the Green function, can then be used in (very + fast) post-processing to compute all required results for different boundary + and initial conditions (and also different internal power sources/imposed + flux). Note that when using the Green function, only boundary and initial + conditions (as well as internal power sources) can be modified: in + particular, the geometry, thermal properties and exchange coefficients have + to remain identical. +- *path visualization*: Stardis-Solver can store the complete spatial and + temporal position along a set of thermal paths, for latter visualization. In + addition of their position and, each thermal path vertex register additional + data as the type of thermal phenomena it simulates, the accumulated + power/flux along the path, etc. + +Stardis-Solver is currently used in two frameworks. The +[Stardis](https://gitlab.com/meso-star/stardis.git) CLI and its associated +tools is the reference workflow of Stardis-Solver. It proposes a complete +toolchain from fileformats describing the scene (geometry, thermal properties, +limit and boundary conditions) to computations and post-treatments of the +results ([Stardis-Green](https://gitlab.com/meso-star/stardis-green.git). +Stardis-Solver is also integrated into +[SYRTHES](https://www.edf.fr/en/the-edf-group/world-s-largest-power-company/activities/research-and-development/scientific-communities/simulation-softwares?logiciel=10818), +the general thermal free software developed by Electricité De France (EDF). ## How to build -Stardis relies on the [CMake](http://www.cmake.org) and the +Stardis-Solver is compatible GNU/Linux as well as Microsoft Windows 7 and +later, both in 64-bits. It was successfully built with the [GNU Compiler +Collection](https://gcc.gnu.org) (versions 4.9.2 and later) as well as with +Microsoft Visual Studio 2015. + +It relies on the [CMake](http://www.cmake.org) and the [RCMake](https://gitlab.com/vaplv/rcmake/) package to build. It also depends on the [RSys](https://gitlab.com/vaplv/rsys/), @@ -25,6 +111,19 @@ variable the install directories of its dependencies. ## Release notes +### Version 0.11 + +- Add support of unsteady green evaluation. The resulting green function can + then be used to quickly evaluate the system at the same time bit with + different limit and initial conditions, volumetric powers and imposed fluxes. +- Add checks on green re-evaluation to ensure that the system remains unchanged + regarding its scale factor and its reference temperature. +- Remove the ambient radiative temperature, the reference temperature and the + geometry scale factor from the list of arguments submitted to the solve + functions. They become scene arguments defined on scene creation. +- Update the `sdis_scene_[2d_]create` function profile: its data are now + grouped into a variable of type `struct sdis_scene_create_args`. + ### Version 0.10.1 - In green function estimation, the time sent to the user callbacks is no more @@ -70,7 +169,7 @@ variable the install directories of its dependencies. data inconsistencies. - Fix a memory leak when the scene creation failed. - Enable parallelism on Star-Enclosure[2D] to improve the performances of the - enclosure extraction on the setup of the Stardis scene. + enclosure extraction on the setup of the Stardis-Solver scene. ### Version 0.8.1 @@ -107,12 +206,12 @@ evaluation. This means that one can estimate the Green function of a system only one time and then evaluate it with different limit conditions, boundary fluxes or power terms with negligible computation costs. -Currently, Stardis assumes that during the Green function estimation, the -properties of the system do not depend on time. In addition, it assumes that -the boundary fluxes and the volumic powers are constants in time and space. -Anyway, on Green function evaluation, the limit conditions of the system can -still vary in time and space; systems in steady state can be simulated with -Green functions. +Currently, Stardis-Solver assumes that during the Green function estimation, +the properties of the system do not depend on time. In addition, it assumes +that the boundary fluxes and the volumetric powers are constants in time and +space. Anyway, on Green function evaluation, the limit conditions of the +system can still vary in time and space; systems in steady state can be +simulated with Green functions. #### Add heat path registration @@ -158,7 +257,7 @@ the few paths that failed in order to diagnose what went wrong. not support time integration, yet. - Add support of an explicit initial time `t0` for the fluid. - Fix a bug in the estimation of unknown fluid temperatures: the associativity - between the internal Stardis data and the user defined data was wrong. + between the internal Stardis-Solver data and the user defined data was wrong. ### Version 0.5 @@ -242,17 +341,17 @@ Full rewrite of how the volumetric power is taken into account. ### Version 0.0 -First version and implementation of the Stardis solver API. +First version and implementation of the Stardis-Solver API. - Support fluid/solid and solid/solid interfaces. - Only conduction is currently fully supported: convection and radiative temperature are not computed yet. Fluid media can be added to the system but - currently, Stardis assumes that their temperature are known. + currently, Stardis-Solver assumes that their temperature are known. ## License -Copyright (C) 2016-2020 |Meso|Star> (<contact@meso-star.com>). Stardis is free -software released under the GPLv3+ license: GNU GPL version 3 or later. You are -welcome to redistribute it under certain conditions; refer to the COPYING files -for details. +Copyright (C) 2016-2020 |Meso|Star> (<contact@meso-star.com>). Stardis-Solver +is free software released under the GPLv3+ license: GNU GPL version 3 or later. +You are welcome to redistribute it under certain conditions; refer to the +COPYING files for details. diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -56,8 +56,8 @@ rcmake_append_runtime_dirs(_runtime_dirs # Configure and define targets ############################################################################### set(VERSION_MAJOR 0) -set(VERSION_MINOR 10) -set(VERSION_PATCH 1) +set(VERSION_MINOR 11) +set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SDIS_FILES_SRC @@ -175,7 +175,7 @@ if(NOT NO_TEST) endfunction() new_test(test_sdis_camera) - new_test(test_sdis_compute_mean_power) + new_test(test_sdis_compute_power) new_test(test_sdis_conducto_radiative) new_test(test_sdis_conducto_radiative_2d) new_test(test_sdis_convection) @@ -213,7 +213,7 @@ if(NOT NO_TEST) add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d) endif() - target_link_libraries(test_sdis_compute_mean_power Star3DUT) + target_link_libraries(test_sdis_compute_power Star3DUT) target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT) target_link_libraries(test_sdis_solve_medium Star3DUT) target_link_libraries(test_sdis_solve_probe3 Star3DUT) diff --git a/src/sdis.h b/src/sdis.h @@ -69,6 +69,7 @@ struct sdis_medium; struct sdis_scene; /* Forward declaration of non ref counted types */ +struct sdis_green_path; struct sdis_heat_path; /******************************************************************************* @@ -270,6 +271,14 @@ typedef res_T /******************************************************************************* * Green function data types ******************************************************************************/ +enum sdis_green_path_end_type { + SDIS_GREEN_PATH_END_AT_INTERFACE, + SDIS_GREEN_PATH_END_IN_VOLUME, + SDIS_GREEN_PATH_END_RADIATIVE, + SDIS_GREEN_PATH_END_TYPES_COUNT__, + SDIS_GREEN_PATH_END_ERROR = SDIS_GREEN_PATH_END_TYPES_COUNT__ +}; + enum sdis_point_type { SDIS_FRAGMENT, SDIS_VERTEX, @@ -277,16 +286,6 @@ enum sdis_point_type { SDIS_POINT_NONE = SDIS_POINT_TYPES_COUNT__ }; -/* Path used to estimate the green function */ -struct sdis_green_path { - /* Internal data. Should not be accessed */ - void* green__; - size_t id__; -}; -#define SDIS_GREEN_PATH_NULL__ {NULL, 0} -static const struct sdis_green_path SDIS_GREEN_PATH_NULL = - SDIS_GREEN_PATH_NULL__; - /* Spatio temporal point */ struct sdis_point { union { @@ -328,15 +327,68 @@ typedef res_T void* context); /******************************************************************************* + * Data types of scene creation + ******************************************************************************/ +/* Functor used to retrieve the indices toward the vertices of a geometric + * primitive. Geometric primitive means for segment in 2D and triangle in 3D */ +typedef void +(*sdis_get_primitive_indices_T) + (const size_t iprim, /* Index of the primitive */ + size_t ids[], /* Output list of primitive indices */ + void* ctx); /* User defined data */ + +/* Retrieve the interface, i.e. the physical properties, associated to a given + * geometric primitive */ +typedef void +(*sdis_get_primitive_interface_T) + (const size_t iprim, /* Index of the primitive */ + struct sdis_interface** interf, + void* ctx); + +/* Retrieve the coordinates of a vertex */ +typedef void +(*sdis_get_vertex_position_T) + (const size_t ivert, /* Index of the vertex */ + double pos[], /* Output list of vertex coordinates */ + void* ctx); + +struct sdis_scene_create_args { + /* Functors to retrieve the geometric description */ + sdis_get_primitive_indices_T get_indices; + sdis_get_primitive_interface_T get_interface; + sdis_get_vertex_position_T get_position; + + /* Pointer toward client side sent as the last argument of the callbacks */ + void* context; + + size_t nprimitives; /* #primitives, i.e. #segments or #triangles */ + size_t nvertices; /* #vertices */ + double fp_to_meter; /* Scale factor used to convert 1.0 in 1 meter */ + double trad; /* Ambiant radiative temperature */ + double tref; /* Temperature used to linearize the radiative temperature */ +}; + +#define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \ + NULL, /* Get indices */ \ + NULL, /* Get interfaces */ \ + NULL, /* Get position */ \ + NULL, /* Context */ \ + 0, /* #primitives */ \ + 0, /* #vertices */ \ + 1.0, /* #Floating point to meter scale factor */ \ + -1.0, /* Ambient radiative temperature */ \ + -1.0 /* Reference temperature */ \ +} +static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT = + SDIS_SCENE_CREATE_ARGS_DEFAULT__; + +/******************************************************************************* * 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 */ }; @@ -344,9 +396,6 @@ struct sdis_solve_probe_args { 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 */ \ } @@ -360,9 +409,6 @@ struct sdis_solve_probe_boundary_args { 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 */ }; @@ -372,9 +418,6 @@ struct sdis_solve_probe_boundary_args { {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 */ \ } @@ -388,9 +431,6 @@ struct sdis_solve_boundary_args { 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 */ }; @@ -400,9 +440,6 @@ struct sdis_solve_boundary_args { 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 */ \ } @@ -413,9 +450,6 @@ 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 */ }; @@ -423,9 +457,6 @@ struct sdis_solve_medium_args { 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 */ \ } @@ -437,9 +468,6 @@ struct sdis_solve_probe_boundary_flux_args { 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__ { \ @@ -447,9 +475,6 @@ struct sdis_solve_probe_boundary_flux_args { 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 @@ -461,9 +486,6 @@ struct sdis_solve_boundary_flux_args { 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__ { \ @@ -471,9 +493,6 @@ struct sdis_solve_boundary_flux_args { 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 @@ -483,9 +502,6 @@ 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 */ @@ -493,9 +509,6 @@ struct sdis_solve_camera_args { #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 \ @@ -507,14 +520,12 @@ struct sdis_compute_power_args { size_t nrealisations; struct sdis_medium* medium; /* Medium to solve */ double time_range[2]; /* Observation time */ - double fp_to_meter; /* Scale from floating point units to meters */ struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */ }; #define SDIS_COMPUTE_POWER_ARGS_DEFAULT__ { \ 10000, /* #realisations */ \ NULL, /* Medium */ \ {DBL_MAX,DBL_MAX}, /* Time range */ \ - 1, /* FP to meter */ \ NULL /* RNG state */ \ } static const struct sdis_compute_power_args @@ -747,7 +758,7 @@ sdis_interface_get_id * * Note that each triangle has 2 sides: a front and a back side. By convention, * the front side of a triangle is the side where its vertices are clock wise - * ordered. The back side of a triangle is the exact opposite: it is the side + * ordered. The back side of a triangle is the exact opposite: it is the side * where the triangle vertices are counter-clock wise ordered. The front and * back media of a triangle interface directly refer to this convention and * thus one has to take care of how the triangle vertices are defined to ensure @@ -755,15 +766,7 @@ sdis_interface_get_id SDIS_API res_T sdis_scene_create (struct sdis_device* dev, - const size_t ntris, /* #triangles */ - void (*indices) /* Retrieve the indices toward the vertices of `itri' */ - (const size_t itri, size_t ids[3], void* ctx), - void (*interf) /* Get the interface of the triangle `itri' */ - (const size_t itri, struct sdis_interface** bound, void* ctx), - const size_t nverts, /* #vertices */ - void (*position) /* Retrieve the position of the vertex `ivert' */ - (const size_t ivert, double pos[3], void* ctx), - void* ctx, /* Client side data sent as input of the previous callbacks */ + const struct sdis_scene_create_args* args, struct sdis_scene** scn); /* Create a 2D scene. The geometry of the 2D scene is defined by an indexed @@ -784,15 +787,7 @@ sdis_scene_create SDIS_API res_T sdis_scene_2d_create (struct sdis_device* dev, - const size_t nsegs, /* #segments */ - void (*indices) /* Retrieve the indices toward the vertices of `iseg' */ - (const size_t iseg, size_t ids[2], void* ctx), - void (*interf) /* Get the interface of the segment `iseg' */ - (const size_t iseg, struct sdis_interface** bound, void* ctx), - const size_t nverts, /* #vertices */ - void (*position) /* Retrieve the position of the vertex `ivert' */ - (const size_t ivert, double pos[2], void* ctx), - void* ctx, /* Client side data sent as input of the previous callbacks */ + const struct sdis_scene_create_args* args, struct sdis_scene** scn); SDIS_API res_T @@ -810,6 +805,44 @@ sdis_scene_get_aabb double lower[3], double upper[3]); +/* Get scene's fp_to_meter */ +SDIS_API res_T +sdis_scene_get_fp_to_meter + (const struct sdis_scene* scn, + double* fp_to_meter); + +/* Set scene's fp_to_meter */ +SDIS_API res_T +sdis_scene_set_fp_to_meter + (struct sdis_scene* scn, + const double fp_to_meter); + +/* Get scene's ambient radiative temperature */ +SDIS_API res_T +sdis_scene_get_ambient_radiative_temperature + (const struct sdis_scene* scn, + double* trad); + +/* Set scene's ambient radiative temperature. If set negative, any sample + * ending in ambient radiative temperature will fail */ +SDIS_API res_T +sdis_scene_set_ambient_radiative_temperature + (struct sdis_scene* scn, + const double trad); + +/* Get scene's reference temperature */ +SDIS_API res_T +sdis_scene_get_reference_temperature + (const struct sdis_scene* scn, + double* tref); + +/* Set scene's reference temperature. If set to 0, there is no radiative + * transfert in the whole system */ +SDIS_API res_T +sdis_scene_set_reference_temperature + (struct sdis_scene* scn, + const double tref); + /* Search the point onto the scene geometry that is the closest of `pos'. The * `radius' parameter controls the maximum search distance around `pos'. The * returned closest point is expressed locally to the geometric primitive onto @@ -993,7 +1026,6 @@ sdis_green_function_ref_put SDIS_API res_T sdis_green_function_solve (struct sdis_green_function* green, - const double time_range[2], /* Observation time */ struct sdis_estimator** estimator); SDIS_API res_T @@ -1007,6 +1039,12 @@ sdis_green_function_create_from_stream FILE* stream, /* Stream into which the green was serialized */ struct sdis_green_function** green); +/* Retrieve the scene used to compute the green function */ +SDIS_API res_T +sdis_green_function_get_scene + (const struct sdis_green_function* green, + struct sdis_scene** scn); + /* Retrieve the number of valid paths used to estimate the green function. It * is actually equal to the number of successful realisations. */ SDIS_API res_T @@ -1028,15 +1066,31 @@ sdis_green_function_for_each_path sdis_process_green_path_T func, void* context); +/* Retrieve the path's elapsed time */ +SDIS_API res_T +sdis_green_path_get_elapsed_time + (struct sdis_green_path* path_handle, + double* elapsed); + +/* Retrieve the path's end type. */ +SDIS_API res_T +sdis_green_path_get_end_type + (struct sdis_green_path* path, + enum sdis_green_path_end_type* type); + /* Retrieve the spatio-temporal end point of a path used to estimate the green - * function. Note that this point went back in time from the relative - * observation time 0. Its time is thus negative; its absolute value - * represents the time spent by the path into the system. */ + * function. Return RES_BAD_OP for paths ending radiative. */ SDIS_API res_T sdis_green_path_get_limit_point (struct sdis_green_path* path, struct sdis_point* pt); +/* Retrieve the green function the path belongs to */ +SDIS_API res_T +sdis_green_path_get_green_function + (struct sdis_green_path* path_handle, + struct sdis_green_function** green); + /* Retrieve the number of "power terms" associated to a path. */ SDIS_API res_T sdis_green_function_get_power_terms_count @@ -1149,18 +1203,16 @@ sdis_compute_power /******************************************************************************* * Green solvers. * - * Currently only steady computations are supported. As a consequence, the - * observation time is always fixed to infinity. - * - * In addition, the green solvers assumes that the interface fluxes are - * constants in time and space. In the same way the volumic power of the solid - * media must be constant in time and space too. Furthermore, note that only - * the interfaces/media that had a flux/volumic power during green estimation - * can update their flux/volumic power value for subsequent + * Note that only the interfaces/media with flux/volumic power defined during + * green estimation can update their flux/volumic power value for subsequent * sdis_green_function_solve invocations: others interfaces/media are * definitely registered against the green function as interfaces/media with no * flux/volumic power. * + * Also note that the green solvers assume that the interface fluxes are + * constant in time and space. The same applies to the volumic power of the + * solid media. + * * If these assumptions are not ensured by the caller, the behavior of the * estimated green function is undefined. ******************************************************************************/ @@ -1191,4 +1243,3 @@ sdis_solve_medium_green_function END_DECLS #endif /* SDIS_H */ - diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h @@ -105,7 +105,6 @@ static const struct XD(rwalk) XD(RWALK_NULL) = { struct XD(temperature) { res_T (*func)/* Next function to invoke in order to compute the temperature */ (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -28,12 +28,17 @@ #include <rsys/hash_table.h> #include <rsys/mem_allocator.h> #include <rsys/ref_count.h> - #include <rsys/dynamic_array.h> -#include <rsys/ref_count.h> #include <limits.h> +/* Path used to estimate the green function */ +struct sdis_green_path { + /* Internal data. Should not be accessed */ + void* green__; + size_t id__; +}; + struct power_term { double term; /* Power term computed during green estimation */ unsigned id; /* Identifier of the medium of the term */ @@ -77,6 +82,7 @@ flux_term_init(struct mem_allocator* allocator, struct flux_term* term) #include <rsys/dynamic_array.h> struct green_path { + double elapsed_time; struct darray_flux_term flux_terms; /* List of flux terms */ struct darray_power_term power_terms; /* List of volumic power terms */ union { @@ -84,7 +90,7 @@ struct green_path { struct sdis_interface_fragment fragment; } limit; unsigned limit_id; /* Identifier of the limit medium/interface */ - enum sdis_point_type limit_type; + enum sdis_green_path_end_type end_type; /* Indices of the last accessed medium/interface. Used to speed up the access * to the medium/interface. */ @@ -98,10 +104,11 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path) ASSERT(path); darray_flux_term_init(allocator, &path->flux_terms); darray_power_term_init(allocator, &path->power_terms); + path->elapsed_time = -INF; path->limit.vertex = SDIS_RWALK_VERTEX_NULL; path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL; path->limit_id = UINT_MAX; - path->limit_type = SDIS_POINT_NONE; + path->end_type = SDIS_GREEN_PATH_END_TYPES_COUNT__; path->ilast_medium = UINT16_MAX; path->ilast_interf = UINT16_MAX; } @@ -119,9 +126,10 @@ green_path_copy(struct green_path* dst, const struct green_path* src) { res_T res = RES_OK; ASSERT(dst && src); + dst->elapsed_time = src->elapsed_time; dst->limit = src->limit; dst->limit_id = src->limit_id; - dst->limit_type = src->limit_type; + dst->end_type = src->end_type; dst->ilast_medium = src->ilast_medium; dst->ilast_interf = src->ilast_interf; res = darray_flux_term_copy(&dst->flux_terms, &src->flux_terms); @@ -136,9 +144,10 @@ green_path_copy_and_clear(struct green_path* dst, struct green_path* src) { res_T res = RES_OK; ASSERT(dst && src); + dst->elapsed_time = src->elapsed_time; dst->limit = src->limit; dst->limit_id = src->limit_id; - dst->limit_type = src->limit_type; + dst->end_type = src->end_type; dst->ilast_medium = src->ilast_medium; dst->ilast_interf = src->ilast_interf; res = darray_flux_term_copy_and_clear(&dst->flux_terms, &src->flux_terms); @@ -154,9 +163,10 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src) { res_T res = RES_OK; ASSERT(dst && src); + dst->elapsed_time = src->elapsed_time; dst->limit = src->limit; dst->limit_id = src->limit_id; - dst->limit_type = src->limit_type; + dst->end_type = src->end_type; dst->ilast_medium = src->ilast_medium; dst->ilast_interf = src->ilast_interf; res = darray_flux_term_copy_and_release(&dst->flux_terms, &src->flux_terms); @@ -180,6 +190,9 @@ green_path_write(const struct green_path* path, FILE* stream) } \ } (void)0 + /* Write elapsed time */ + WRITE(&path->elapsed_time, 1); + /* Write the list of flux terms */ sz = darray_flux_term_size_get(&path->flux_terms); WRITE(&sz, 1); @@ -193,7 +206,7 @@ green_path_write(const struct green_path* path, FILE* stream) /* Write the limit point */ WRITE(&path->limit, 1); WRITE(&path->limit_id, 1); - WRITE(&path->limit_type, 1); + WRITE(&path->end_type, 1); /* Write miscellaneous data */ WRITE(&path->ilast_medium, 1); @@ -227,6 +240,9 @@ green_path_read(struct green_path* path, FILE* stream) } \ } (void)0 + /* Read elapsed time */ + READ(&path->elapsed_time, 1); + /* Read the list of flux terms */ READ(&sz, 1); res = darray_flux_term_resize(&path->flux_terms, sz); @@ -242,7 +258,7 @@ green_path_read(struct green_path* path, FILE* stream) /* Read the limit point */ READ(&path->limit, 1); READ(&path->limit_id, 1); - READ(&path->limit_type, 1); + READ(&path->end_type, 1); /* Read the miscellaneous data */ READ(&path->ilast_medium, 1); @@ -400,7 +416,6 @@ green_function_fetch_interf static res_T green_function_solve_path (struct sdis_green_function* green, - const double time, /* Sampled time */ const size_t ipath, double* weight) { @@ -409,19 +424,18 @@ green_function_solve_path const struct green_path* path = NULL; const struct sdis_medium* medium = NULL; const struct sdis_interface* interf = NULL; + struct sdis_scene* scn = NULL; struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; double power; double flux; - double temperature; - double time_curr; + double end_temperature; size_t i, n; res_T res = RES_OK; ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight); - ASSERT(time > 0); path = darray_green_path_cdata_get(&green->paths) + ipath; - if(path->limit_type == SDIS_POINT_NONE) { /* Rejected path */ + if(path->end_type == SDIS_GREEN_PATH_END_ERROR) { /* Rejected path */ res = RES_BAD_OP; goto error; } @@ -447,46 +461,31 @@ green_function_solve_path flux += flux_terms[i].term * interface_side_get_flux(interf, &frag); } - /* Setup time. */ - switch(path->limit_type) { - case SDIS_FRAGMENT: - time_curr = time + path->limit.fragment.time; + /* Compute path's end temperature */ + switch(path->end_type) { + case SDIS_GREEN_PATH_END_AT_INTERFACE: interf = green_function_fetch_interf(green, path->limit_id); - break; - case SDIS_VERTEX: - time_curr = time + path->limit.vertex.time; - medium = green_function_fetch_medium(green, path->limit_id); - break; - default: FATAL("Unreachable code.\n"); break; - } - - if(time_curr <= 0 - || (path->limit_type == SDIS_VERTEX && time_curr <= medium_get_t0(medium))) { - log_err(green->scn->dev, - "%s: invalid observation time \"%g\": the initial condition is reached " - "while instationary system are not supported by the green function.\n", - FUNC_NAME, time); - res = RES_BAD_ARG; - goto error; - } - - /* Compute limit condition */ - switch(path->limit_type) { - case SDIS_FRAGMENT: frag = path->limit.fragment; - frag.time = time_curr; - temperature = interface_side_get_temperature(interf, &frag); + end_temperature = interface_side_get_temperature(interf, &frag); break; - case SDIS_VERTEX: + case SDIS_GREEN_PATH_END_IN_VOLUME: + medium = green_function_fetch_medium(green, path->limit_id); vtx = path->limit.vertex; - vtx.time = time_curr; - temperature = medium_get_temperature(medium, &vtx); + end_temperature = medium_get_temperature(medium, &vtx); + break; + case SDIS_GREEN_PATH_END_RADIATIVE: + SDIS(green_function_get_scene(green, &scn)); + SDIS(scene_get_ambient_radiative_temperature(scn, &end_temperature)); + if(end_temperature < 0) { /* Cannot be negative if used */ + res = RES_BAD_ARG; + goto error; + } break; default: FATAL("Unreachable code.\n"); break; } /* Compute the path weight */ - *weight = power + flux + temperature; + *weight = power + flux + end_temperature; exit: return res; @@ -840,11 +839,9 @@ sdis_green_function_ref_put(struct sdis_green_function* green) res_T sdis_green_function_solve (struct sdis_green_function* green, - const double time_range[2], struct sdis_estimator** out_estimator) { struct sdis_estimator* estimator = NULL; - struct ssp_rng* rng = NULL; size_t npaths; size_t ipath; size_t N = 0; /* #realisations */ @@ -852,21 +849,11 @@ sdis_green_function_solve double accum2 = 0; res_T res = RES_OK; - if(!green || !time_range || time_range[0] < 0 - || time_range[1] < time_range[0] || !out_estimator) { + if(!green || !out_estimator) { res = RES_BAD_ARG; goto error; } - res = ssp_rng_create(green->scn->dev->allocator, &green->rng_type, &rng); - if(res != RES_OK) goto error; - - /* Avoid correlation by defining the RNG state from the final state of the - * RNG used to estimate the green function */ - rewind(green->rng_state); - res = ssp_rng_read(rng, green->rng_state); - if(res != RES_OK) goto error; - npaths = darray_green_path_size_get(&green->paths); /* Create the estimator */ @@ -875,10 +862,9 @@ sdis_green_function_solve /* Solve the green function */ FOR_EACH(ipath, 0, npaths) { - const double time = sample_time(rng, time_range); double w; - res = green_function_solve_path(green, time, ipath, &w); + res = green_function_solve_path(green, ipath, &w); if(res == RES_BAD_OP) continue; if(res != RES_OK) goto error; @@ -894,7 +880,6 @@ sdis_green_function_solve (estimator, green->realisation_time.sum, green->realisation_time.sum2); exit: - if(rng) SSP(rng_ref_put(rng)); if(out_estimator) *out_estimator = estimator; return res; error: @@ -1019,8 +1004,8 @@ sdis_green_function_create_from_stream READ(hash1, sizeof(hash256_T)); if(!hash256_eq(hash0, hash1)) { log_err(green->scn->dev, - "%s: the submitted scene does not match scene used to estimate the green " - "function.\n", FUNC_NAME); + "%s: the submitted scene does not match the scene used to estimate the " + "green function.\n", FUNC_NAME); res = RES_BAD_ARG; goto error; } @@ -1061,6 +1046,17 @@ error: } res_T +sdis_green_function_get_scene + (const struct sdis_green_function* green, + struct sdis_scene** scn) +{ + if(!green || !scn) return RES_BAD_ARG; + ASSERT(green->npaths_valid != SIZE_MAX); + *scn = green->scn; + return RES_OK; +} + +res_T sdis_green_function_get_paths_count (const struct sdis_green_function* green, size_t* npaths) { @@ -1100,7 +1096,7 @@ sdis_green_function_for_each_path struct sdis_green_path path_handle; const struct green_path* path = darray_green_path_cdata_get(&green->paths)+ipath; - if(path->limit_type == SDIS_POINT_NONE) continue; + if(path->end_type == SDIS_GREEN_PATH_END_ERROR) continue; path_handle.green__ = green; path_handle.id__ = ipath; @@ -1116,6 +1112,56 @@ error: } res_T +sdis_green_path_get_elapsed_time + (struct sdis_green_path* path_handle, double* elapsed) +{ + const struct green_path* path = NULL; + struct sdis_green_function* green = NULL; + res_T res = RES_OK; + + if(!path_handle || !elapsed) { + res = RES_BAD_ARG; + goto error; + } + + green = path_handle->green__; + ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); + + path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; + *elapsed = path->elapsed_time; + +exit: + return res; +error: + goto exit; +} + +res_T +sdis_green_path_get_end_type + (struct sdis_green_path* path_handle, enum sdis_green_path_end_type* type) +{ + const struct green_path* path = NULL; + struct sdis_green_function* green = NULL; + res_T res = RES_OK; + + if(!path_handle || !type) { + res = RES_BAD_ARG; + goto error; + } + + green = path_handle->green__; + ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); + + path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; + *type = path->end_type; + +exit: + return res; +error: + goto exit; +} + +res_T sdis_green_path_get_limit_point (struct sdis_green_path* path_handle, struct sdis_point* pt) { @@ -1132,16 +1178,21 @@ sdis_green_path_get_limit_point ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; - pt->type = path->limit_type; - switch(path->limit_type) { - case SDIS_FRAGMENT: + switch(path->end_type) { + case SDIS_GREEN_PATH_END_AT_INTERFACE: pt->data.itfrag.intface = green_function_fetch_interf(green, path->limit_id); pt->data.itfrag.fragment = path->limit.fragment; + pt->type = SDIS_FRAGMENT; break; - case SDIS_VERTEX: + case SDIS_GREEN_PATH_END_IN_VOLUME: pt->data.mdmvert.medium = green_function_fetch_medium(green, path->limit_id); pt->data.mdmvert.vertex = path->limit.vertex; + pt->type = SDIS_VERTEX; + break; + case SDIS_GREEN_PATH_END_RADIATIVE: + res = RES_BAD_OP; + goto error; break; default: FATAL("Unreachable code.\n"); break; } @@ -1153,6 +1204,31 @@ error: } res_T +sdis_green_path_get_green_function + (struct sdis_green_path* path_handle, + struct sdis_green_function** out_green) + +{ + struct sdis_green_function* green = NULL; + res_T res = RES_OK; + + if(!path_handle || !out_green) { + res = RES_BAD_ARG; + goto error; + } + + green = path_handle->green__; + ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); + + *out_green = green; + +exit: + return res; +error: + goto exit; +} + +res_T sdis_green_function_get_power_terms_count (const struct sdis_green_path* path_handle, size_t* nterms) @@ -1427,7 +1503,7 @@ green_function_finalize n = darray_green_path_size_get(&green->paths); FOR_EACH(i, 0, n) { const struct green_path* path = darray_green_path_cdata_get(&green->paths)+i; - green->npaths_valid += path->limit_type != SDIS_POINT_NONE; + green->npaths_valid += (path->end_type != SDIS_GREEN_PATH_END_ERROR); } green->npaths_invalid = n - green->npaths_valid; @@ -1467,13 +1543,13 @@ green_path_set_limit_interface_fragment { res_T res = RES_OK; ASSERT(handle && interf && frag); - ASSERT(handle->path->limit_type == SDIS_POINT_NONE); + ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); res = ensure_interface_registration(handle->green, interf); if(res != RES_OK) return res; + handle->path->elapsed_time = elapsed_time; handle->path->limit.fragment = *frag; - handle->path->limit.fragment.time = -elapsed_time; handle->path->limit_id = interface_get_id(interf); - handle->path->limit_type = SDIS_FRAGMENT; + handle->path->end_type = SDIS_GREEN_PATH_END_AT_INTERFACE; return RES_OK; } @@ -1486,13 +1562,25 @@ green_path_set_limit_vertex { res_T res = RES_OK; ASSERT(handle && mdm && vert); - ASSERT(handle->path->limit_type == SDIS_POINT_NONE); + ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); res = ensure_medium_registration(handle->green, mdm); if(res != RES_OK) return res; + handle->path->elapsed_time = elapsed_time; handle->path->limit.vertex = *vert; - handle->path->limit.vertex.time = -elapsed_time; handle->path->limit_id = medium_get_id(mdm); - handle->path->limit_type = SDIS_VERTEX; + handle->path->end_type = SDIS_GREEN_PATH_END_IN_VOLUME; + return RES_OK; +} + +res_T +green_path_set_limit_radiative + (struct green_path_handle* handle, + const double elapsed_time) +{ + ASSERT(handle); + ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); + handle->path->elapsed_time = elapsed_time; + handle->path->end_type = SDIS_GREEN_PATH_END_RADIATIVE; return RES_OK; } @@ -1612,4 +1700,3 @@ exit: error: goto exit; } - diff --git a/src/sdis_green.h b/src/sdis_green.h @@ -21,7 +21,7 @@ /* Current version the green function data structure. One should increment it * and perform a version management onto serialized data when the gren function * data structure is updated. */ -static const int SDIS_GREEN_FUNCTION_VERSION = 0; +static const int SDIS_GREEN_FUNCTION_VERSION = 1; /* Forward declaration */ struct accum; @@ -82,6 +82,11 @@ green_path_set_limit_vertex const double elapsed_time); extern LOCAL_SYM res_T +green_path_set_limit_radiative + (struct green_path_handle* handle, + const double elapsed_time); + +extern LOCAL_SYM res_T green_path_add_power_term (struct green_path_handle* path, struct sdis_medium* mdm, diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -115,7 +115,6 @@ extern LOCAL_SYM res_T trace_radiative_path_2d (struct sdis_scene* scn, const float ray_dir[3], - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct ssp_rng* rng, @@ -125,7 +124,6 @@ extern LOCAL_SYM res_T trace_radiative_path_3d (struct sdis_scene* scn, const float ray_dir[3], - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct ssp_rng* rng, @@ -134,7 +132,6 @@ trace_radiative_path_3d extern LOCAL_SYM res_T radiative_path_2d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct ssp_rng* rng, @@ -143,7 +140,6 @@ radiative_path_2d extern LOCAL_SYM res_T radiative_path_3d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct ssp_rng* rng, @@ -155,7 +151,6 @@ radiative_path_3d extern LOCAL_SYM res_T convective_path_2d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct ssp_rng* rng, @@ -164,7 +159,6 @@ convective_path_2d extern LOCAL_SYM res_T convective_path_3d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct ssp_rng* rng, @@ -176,7 +170,6 @@ convective_path_3d extern LOCAL_SYM res_T conductive_path_2d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct ssp_rng* rng, @@ -185,7 +178,6 @@ conductive_path_2d extern LOCAL_SYM res_T conductive_path_3d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct ssp_rng* rng, @@ -197,7 +189,6 @@ conductive_path_3d extern LOCAL_SYM res_T boundary_path_2d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct ssp_rng* rng, @@ -206,7 +197,6 @@ boundary_path_2d extern LOCAL_SYM res_T boundary_path_3d (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct ssp_rng* rng, diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -458,7 +458,6 @@ XD(check_rwalk_fragment_consistency) static res_T XD(solid_solid_boundary_path) (const struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, const struct sdis_interface_fragment* frag, struct XD(rwalk)* rwalk, @@ -491,7 +490,7 @@ XD(solid_solid_boundary_path) int move; int reinjection_is_valid; res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T); + ASSERT(scn && ctx && frag && rwalk && rng && T); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); (void)frag, (void)ctx; @@ -589,7 +588,7 @@ XD(solid_solid_boundary_path) /* Handle the volumic power */ power = solid_get_volumic_power(mdm, &rwalk->vtx); if(power != SDIS_VOLUMIC_POWER_NONE) { - const double delta_in_meter = reinject_dst * fp_to_meter; + const double delta_in_meter = reinject_dst * scn->fp_to_meter; const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx); tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); T->value += power * tmp; @@ -601,7 +600,7 @@ XD(solid_solid_boundary_path) } /* Time rewind */ - res = XD(time_rewind)(mdm, rng, reinject_dst, fp_to_meter, ctx, rwalk, T); + res = XD(time_rewind)(mdm, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T); if(res != RES_OK) goto error; if(T->done) goto exit; /* Limit condition was reached */ @@ -633,7 +632,6 @@ error: static res_T XD(solid_fluid_boundary_path) (const struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, const struct sdis_interface_fragment* frag, struct XD(rwalk)* rwalk, @@ -666,7 +664,7 @@ XD(solid_fluid_boundary_path) int iattempt; int reinjection_is_valid = 0; res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx); + ASSERT(scn && rwalk && rng && T && ctx); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); /* Retrieve the solid and the fluid split by the boundary */ @@ -741,7 +739,7 @@ XD(solid_fluid_boundary_path) hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon; /* Compute the probas to switch in solid, fluid or radiative random walk */ - tmp = lambda / (delta*fp_to_meter); + tmp = lambda / (delta * scn->fp_to_meter); fluid_proba = hc / (tmp + hr + hc); radia_proba = hr / (tmp + hr + hc); /*solid_proba = tmp / (tmp + hr + hc);*/ @@ -759,7 +757,7 @@ XD(solid_fluid_boundary_path) /* Handle the volumic power */ const double power = solid_get_volumic_power(solid, &rwalk->vtx); if(power != SDIS_VOLUMIC_POWER_NONE) { - const double delta_in_meter = reinject_dst * fp_to_meter; + const double delta_in_meter = reinject_dst * scn->fp_to_meter; tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); T->value += power * tmp; @@ -770,7 +768,7 @@ XD(solid_fluid_boundary_path) } /* Time rewind */ - res = XD(time_rewind)(solid, rng, reinject_dst, fp_to_meter, ctx, rwalk, T); + res = XD(time_rewind)(solid, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T); if(res != RES_OK) goto error; if(T->done) goto exit; /* Limit condition was reached */ @@ -803,7 +801,6 @@ error: static res_T XD(solid_boundary_with_flux_path) (const struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, const struct sdis_interface_fragment* frag, const double phi, @@ -888,7 +885,7 @@ XD(solid_boundary_with_flux_path) delta = reinject_dst / sqrt(DIM); /* Handle the flux */ - delta_in_meter = delta*fp_to_meter; + delta_in_meter = delta * scn->fp_to_meter; tmp = delta_in_meter / lambda; T->value += phi * tmp; if(ctx->green_path) { @@ -899,7 +896,7 @@ XD(solid_boundary_with_flux_path) /* Handle the volumic power */ power = solid_get_volumic_power(mdm, &rwalk->vtx); if(power != SDIS_VOLUMIC_POWER_NONE) { - delta_in_meter = reinject_dst * fp_to_meter; + delta_in_meter = reinject_dst * scn->fp_to_meter; tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); T->value += power * tmp; if(ctx->green_path) { @@ -909,7 +906,7 @@ XD(solid_boundary_with_flux_path) } /* Time rewind */ - res = XD(time_rewind)(mdm, rng, reinject_dst, fp_to_meter, ctx, rwalk, T); + res = XD(time_rewind)(mdm, rng, reinject_dst * scn->fp_to_meter, ctx, rwalk, T); if(res != RES_OK) goto error; if(T->done) goto exit; /* Limit condition was reached */ @@ -945,7 +942,6 @@ error: res_T XD(boundary_path) (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -958,7 +954,7 @@ XD(boundary_path) struct sdis_medium* mdm = NULL; double tmp; res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); + ASSERT(scn && ctx && rwalk && rng && T); ASSERT(rwalk->mdm == NULL); ASSERT(!SXD_HIT_NONE(&rwalk->hit)); @@ -993,7 +989,7 @@ XD(boundary_path) const double phi = interface_side_get_flux(interf, &frag); if(phi != SDIS_FLUX_NONE) { res = XD(solid_boundary_with_flux_path) - (scn, fp_to_meter, ctx, &frag, phi, rwalk, rng, T); + (scn, ctx, &frag, phi, rwalk, rng, T); if(res != RES_OK) goto error; goto exit; @@ -1004,11 +1000,9 @@ XD(boundary_path) mdm_back = interface_get_medium(interf, SDIS_BACK); if(mdm_front->type == mdm_back->type) { - res = XD(solid_solid_boundary_path) - (scn, fp_to_meter, ctx, &frag, rwalk, rng, T); + res = XD(solid_solid_boundary_path)(scn, ctx, &frag, rwalk, rng, T); } else { - res = XD(solid_fluid_boundary_path) - (scn, fp_to_meter, ctx, &frag, rwalk, rng, T); + res = XD(solid_fluid_boundary_path)(scn, ctx, &frag, rwalk, rng, T); } if(res != RES_OK) goto error; diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -199,7 +199,6 @@ error: res_T XD(conductive_path) (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -211,7 +210,7 @@ XD(conductive_path) struct sdis_medium* mdm; size_t istep = 0; /* Help for debug */ res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T); + ASSERT(scn && rwalk && rng && T); ASSERT(rwalk->mdm->type == SDIS_SOLID); (void)ctx, (void)istep; @@ -286,12 +285,12 @@ XD(conductive_path) /* Add the volumic power density to the measured temperature */ if(power != SDIS_VOLUMIC_POWER_NONE) { if((S3D_HIT_NONE(&hit0) && S3D_HIT_NONE(&hit1))) { /* Hit nothing */ - const double delta_in_meter = delta * fp_to_meter; + const double delta_in_meter = delta * scn->fp_to_meter; power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); T->value += power * power_factor; } else { const double delta_s_adjusted = delta_solid * RAY_RANGE_MAX_SCALE; - const double delta_s_in_meter = delta_solid * fp_to_meter; + const double delta_s_in_meter = delta_solid * scn->fp_to_meter; double h; double h_in_meter; double cos_U_N; @@ -307,7 +306,7 @@ XD(conductive_path) } h = delta * fabs(cos_U_N); - h_in_meter = h * fp_to_meter; + h_in_meter = h * scn->fp_to_meter; /* The regular power term at wall */ tmp = h_in_meter * h_in_meter / (2.0 * lambda); @@ -343,7 +342,7 @@ XD(conductive_path) } /* Rewind the time */ - res = XD(time_rewind)(rwalk->mdm, rng, delta, fp_to_meter, ctx, rwalk, T); + res = XD(time_rewind)(rwalk->mdm, rng, delta * scn->fp_to_meter, ctx, rwalk, T); if(res != RES_OK) goto error; if(T->done) break; /* Limit condition was reached */ diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -70,7 +70,6 @@ XD(register_heat_vertex_in_fluid) res_T XD(convective_path) (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -93,8 +92,8 @@ XD(convective_path) float st[2]; #endif res_T res = RES_OK; - (void)rng, (void)fp_to_meter, (void)ctx; - ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); + (void)rng, (void)ctx; + ASSERT(scn && ctx && rwalk && rng && T); ASSERT(rwalk->mdm->type == SDIS_FLUID); tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx); diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -33,7 +33,6 @@ res_T XD(trace_radiative_path) (struct sdis_scene* scn, const float ray_dir[3], - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -45,8 +44,7 @@ XD(trace_radiative_path) float dir[3] = {0, 0, 0}; res_T res = RES_OK; - ASSERT(scn && ray_dir && fp_to_meter > 0 && ctx && rwalk && rng && T); - (void)fp_to_meter; + ASSERT(scn && ray_dir && ctx && rwalk && rng && T); f3_set(dir, ray_dir); @@ -81,11 +79,8 @@ XD(trace_radiative_path) T->done = 1; if(ctx->green_path) { - struct sdis_rwalk_vertex vtx; - d3_splat(vtx.P, INF); - vtx.time = rwalk->vtx.time; - res = green_path_set_limit_vertex - (ctx->green_path, rwalk->mdm, &vtx, rwalk->elapsed_time); + res = green_path_set_limit_radiative + (ctx->green_path, rwalk->elapsed_time); if(res != RES_OK) goto error; } if(ctx->heat_path) { @@ -193,7 +188,6 @@ error: res_T XD(radiative_path) (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -205,9 +199,8 @@ XD(radiative_path) float dir[3] = {0, 0, 0}; res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); + ASSERT(scn && ctx && rwalk && rng && T); ASSERT(!SXD_HIT_NONE(&rwalk->hit)); - (void)fp_to_meter; /* Normalize the normal of the interface and ensure that it points toward the * current medium */ @@ -220,7 +213,7 @@ XD(radiative_path) ssp_ran_hemisphere_cos_float(rng, N, dir, NULL); /* Launch the radiative random walk */ - res = XD(trace_radiative_path)(scn, dir, fp_to_meter, ctx, rwalk, rng, T); + res = XD(trace_radiative_path)(scn, dir, ctx, rwalk, rng, T); if(res != RES_OK) goto error; exit: diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -136,20 +136,18 @@ register_heat_vertex extern LOCAL_SYM res_T time_rewind_2d - (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */ + (struct sdis_medium* mdm, /* Medium into which the time is rewinded */ struct ssp_rng* rng, - const double delta, - const double fp_to_meter, + const double dist_in_meter, const struct rwalk_context* ctx, struct rwalk_2d* rwalk, struct temperature_2d* T); extern LOCAL_SYM res_T time_rewind_3d - (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */ + (struct sdis_medium* mdm, /* Medium into which the time is rewinded */ struct ssp_rng* rng, - const double delta, - const double fp_to_meter, + const double dist_in_meter, const struct rwalk_context* ctx, struct rwalk_3d* rwalk, struct temperature_3d* T); diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h @@ -17,6 +17,7 @@ #include "sdis_log.h" #include "sdis_medium_c.h" #include "sdis_misc.h" +#include "sdis_green.h" #include <star/ssp.h> @@ -24,20 +25,18 @@ res_T XD(time_rewind) - (const struct sdis_medium* mdm, + (struct sdis_medium* mdm, struct ssp_rng* rng, - const double delta, - const double fp_to_meter, + const double dist_in_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct XD(temperature)* T) { - const double delta_in_meter = delta * fp_to_meter; double temperature; double lambda, rho, cp; double tau, mu, t0; res_T res = RES_OK; - ASSERT(mdm && rng && delta && fp_to_meter && ctx && rwalk); + ASSERT(mdm && rng && ctx && rwalk && dist_in_meter > 0); ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID); ASSERT(T->done == 0); @@ -48,11 +47,11 @@ XD(time_rewind) t0 = solid_get_t0(mdm); /* Limit time */ /* Sample the time to reroll */ - mu = (2*DIM*lambda)/(rho*cp*delta_in_meter*delta_in_meter); + mu = (2*DIM*lambda)/(rho*cp*dist_in_meter*dist_in_meter); tau = ssp_ran_exp(rng, mu); /* Increment the elapsed time */ - ASSERT(rwalk->vtx.time > t0); + ASSERT(rwalk->vtx.time >= t0); rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0); if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */ @@ -84,6 +83,12 @@ XD(time_rewind) vtx->weight = T->value; } + if(ctx->green_path) { + res = green_path_set_limit_vertex(ctx->green_path, mdm, &rwalk->vtx, + rwalk->elapsed_time); + if(res != RES_OK) goto error; + } + exit: return res; error: diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c @@ -29,9 +29,6 @@ ray_realisation_3d const double position[], const double direction[], const double time, - const double fp_to_meter, - const double Tarad, - const double Tref, struct sdis_heat_path* heat_path, /* May be NULL */ double* weight) { @@ -40,8 +37,8 @@ ray_realisation_3d struct temperature_3d T = TEMPERATURE_NULL_3d; float dir[3]; res_T res = RES_OK; - ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight); - ASSERT(Tref >= 0 && medium && medium->type == SDIS_FLUID); + ASSERT(scn && position && direction && time>=0 && weight); + ASSERT(medium && medium->type == SDIS_FLUID); d3_set(rwalk.vtx.P, position); rwalk.vtx.time = time; @@ -49,8 +46,9 @@ ray_realisation_3d rwalk.hit_side = SDIS_SIDE_NULL__; rwalk.mdm = medium; - ctx.Tarad = Tarad; - ctx.Tref3 = Tref*Tref*Tref; + ctx.Tarad = scn->ambient_radiative_temperature; + ctx.Tref3 = scn->reference_temperature * scn->reference_temperature + * scn->reference_temperature; ctx.heat_path = heat_path; f3_set_d3(dir, direction); @@ -59,11 +57,11 @@ ray_realisation_3d res = register_heat_vertex(heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE); if(res != RES_OK) goto error; - res = trace_radiative_path_3d(scn, dir, fp_to_meter, &ctx, &rwalk, rng, &T); + res = trace_radiative_path_3d(scn, dir, &ctx, &rwalk, rng, &T); if(res != RES_OK) goto error; if(!T.done) { - res = compute_temperature_3d(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = compute_temperature_3d(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) goto error; } diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h @@ -43,9 +43,6 @@ probe_realisation_2d struct sdis_medium* medium, const double position[2], const double time, - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double ambient_radiative_temperature, - const double reference_temperature, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight); @@ -58,9 +55,6 @@ probe_realisation_3d struct sdis_medium* medium, const double position[3], const double time, - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double ambient_radiative_temperature, - const double reference_temperature, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight); @@ -76,9 +70,6 @@ boundary_realisation_2d const double u[1], const double time, const enum sdis_side side, - const double fp_to_meter, - const double ambient_radiative_temperature, - const double reference_temperature, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight); @@ -91,9 +82,6 @@ boundary_realisation_3d const double uv[2], const double time, const enum sdis_side side, - const double fp_to_meter, - const double ambient_radiative_temperature, - const double reference_temperature, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight); @@ -106,9 +94,6 @@ boundary_flux_realisation_2d const double uv[1], const double time, const enum sdis_side solid_side, - const double fp_to_meter, - const double ambient_radiative_temperature, - const double reference_temperature, const int flux_mask, /* Combination of enum flux_flag */ double weight[FLUX_NAMES_COUNT__]); @@ -120,9 +105,6 @@ boundary_flux_realisation_3d const double uv[2], const double time, const enum sdis_side solid_side, - const double fp_to_meter, - const double ambient_radiative_temperature, - const double reference_temperature, const int flux_mask, /* Combination of enum flux_flag */ double weight[FLUX_NAMES_COUNT__]); @@ -137,9 +119,6 @@ ray_realisation_3d const double position[3], const double direction[3], const double time, - const double fp_to_meter, - const double ambient_radiative_temperature, - const double reference_temperature, struct sdis_heat_path* heat_path, /* May be NULL */ double* weight); diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -32,7 +32,6 @@ static res_T XD(compute_temperature) (struct sdis_scene* scn, - const double fp_to_meter, const struct rwalk_context* ctx, struct XD(rwalk)* rwalk, struct ssp_rng* rng, @@ -50,7 +49,7 @@ XD(compute_temperature) /* Maximum accepted #failures before stopping the realisation */ const size_t MAX_FAILS = 1; res_T res = RES_OK; - ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); + ASSERT(scn && ctx && rwalk && rng && T); if(ctx->heat_path && T->func == XD(boundary_path)) { heat_vtx = heat_path_get_last_vertex(ctx->heat_path); @@ -72,7 +71,7 @@ XD(compute_temperature) /* Reject the step if a BAD_OP occurs and retry up to MAX_FAILS times */ do { - res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T); + res = T->func(scn, ctx, rwalk, rng, T); if(res == RES_BAD_OP) { *rwalk = rwalk_bkp; *T = T_bkp; } } while(res == RES_BAD_OP && ++nfails < MAX_FAILS); if(res != RES_OK) goto error; @@ -119,9 +118,6 @@ XD(probe_realisation) struct sdis_medium* medium, const double position[], const double time, - const double fp_to_meter,/* Scale factor from floating point unit to meter */ - const double ambient_radiative_temperature, - const double reference_temperature, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight) @@ -135,7 +131,7 @@ XD(probe_realisation) (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx); res_T res = RES_OK; - ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0); + ASSERT(medium && position && weight && time >= 0); (void)irealisation; switch(medium->type) { @@ -186,15 +182,16 @@ XD(probe_realisation) ctx.green_path = green_path; ctx.heat_path = heat_path; - ctx.Tarad = ambient_radiative_temperature; + ctx.Tarad = scn->ambient_radiative_temperature; ctx.Tref3 = - reference_temperature - * reference_temperature - * reference_temperature; + scn->reference_temperature + * scn->reference_temperature + * scn->reference_temperature; - res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) goto error; + ASSERT(T.value >= 0); *weight = T.value; exit: @@ -211,9 +208,6 @@ XD(boundary_realisation) const double uv[2], const double time, const enum sdis_side side, - const double fp_to_meter, - const double Tarad, - const double Tref, struct green_path_handle* green_path, /* May be NULL */ struct sdis_heat_path* heat_path, /* May be NULL */ double* weight) @@ -228,7 +222,7 @@ XD(boundary_realisation) float st[2]; #endif res_T res = RES_OK; - ASSERT(uv && fp_to_meter > 0 && weight && time >= 0); + ASSERT(uv && weight && time >= 0); T.func = XD(boundary_path); rwalk.hit_side = side; @@ -266,10 +260,11 @@ XD(boundary_realisation) ctx.green_path = green_path; ctx.heat_path = heat_path; - ctx.Tarad = Tarad; - ctx.Tref3 = Tref*Tref*Tref; + ctx.Tarad = scn->ambient_radiative_temperature; + ctx.Tref3 = scn->reference_temperature * scn->reference_temperature + * scn->reference_temperature; - res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) goto error; *weight = T.value; @@ -288,9 +283,6 @@ XD(boundary_flux_realisation) const double uv[DIM], const double time, const enum sdis_side solid_side, - const double fp_to_meter, - const double Tarad, - const double Tref, const int flux_mask, double weight[3]) { @@ -309,13 +301,14 @@ XD(boundary_flux_realisation) #endif double P[SDIS_XD_DIMENSION]; float N[SDIS_XD_DIMENSION]; - const double Tr3 = Tref * Tref * Tref; + const double Tr3 = scn->reference_temperature * scn->reference_temperature + * scn->reference_temperature; const enum sdis_side fluid_side = (solid_side == SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT; res_T res = RES_OK; const char compute_radiative = (flux_mask & FLUX_FLAG_RADIATIVE) != 0; const char compute_convective = (flux_mask & FLUX_FLAG_CONVECTIVE) != 0; - ASSERT(uv && fp_to_meter > 0 && weight && time >= 0 && Tref >= 0); + ASSERT(uv && weight && time >= 0 ); #if SDIS_XD_DIMENSION == 2 #define SET_PARAM(Dest, Src) (Dest).u = (Src); @@ -344,7 +337,7 @@ XD(boundary_flux_realisation) rwalk.mdm = (Mdm); \ rwalk.hit.prim = prim; \ SET_PARAM(rwalk.hit, st); \ - ctx.Tarad = Tarad; \ + ctx.Tarad = scn->ambient_radiative_temperature; \ ctx.Tref3 = Tr3; \ dX(set)(rwalk.vtx.P, P); \ fX(set)(rwalk.hit.normal, N); \ @@ -354,7 +347,7 @@ XD(boundary_flux_realisation) /* Compute boundary temperature */ RESET_WALK(solid_side, NULL); T.func = XD(boundary_path); - res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; weight[0] = T.value; @@ -366,7 +359,7 @@ XD(boundary_flux_realisation) if(compute_radiative) { RESET_WALK(fluid_side, fluid_mdm); T.func = XD(radiative_path); - res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; weight[1] = T.value; } @@ -375,7 +368,7 @@ XD(boundary_flux_realisation) if(compute_convective) { RESET_WALK(fluid_side, fluid_mdm); T.func = XD(convective_path); - res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; weight[2] = T.value; } diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -125,31 +125,19 @@ scene_release(ref_T * ref) res_T sdis_scene_create (struct sdis_device* dev, - const size_t ntris, /* #triangles */ - void (*indices)(const size_t itri, size_t ids[3], void*), - void (*interf)(const size_t itri, struct sdis_interface** bound, void*), - const size_t nverts, /* #vertices */ - void (*position)(const size_t ivert, double pos[3], void* ctx), - void* ctx, + const struct sdis_scene_create_args* args, struct sdis_scene** out_scn) { - return scene_create_3d - (dev, ntris, indices, interf, nverts, position, ctx, out_scn); + return scene_create_3d(dev, args, out_scn); } res_T sdis_scene_2d_create (struct sdis_device* dev, - const size_t nsegs, /* #segments */ - void (*indices)(const size_t iseg, size_t ids[2], void*), - void (*interf)(const size_t iseg, struct sdis_interface** bound, void*), - const size_t nverts, /* #vertices */ - void (*position)(const size_t ivert, double pos[2], void* ctx), - void* ctx, + const struct sdis_scene_create_args* args, struct sdis_scene** out_scn) { - return scene_create_2d - (dev, nsegs, indices, interf, nverts, position, ctx, out_scn); + return scene_create_2d(dev, args, out_scn); } res_T @@ -191,6 +179,66 @@ sdis_scene_get_aabb } res_T +sdis_scene_get_fp_to_meter + (const struct sdis_scene* scn, + double* fp_to_meter) +{ + if(!scn || !fp_to_meter) return RES_BAD_ARG; + *fp_to_meter = scn->fp_to_meter; + return RES_OK; +} + +res_T +sdis_scene_set_fp_to_meter + (struct sdis_scene* scn, + const double fp_to_meter) +{ + if(!scn || fp_to_meter <= 0) return RES_BAD_ARG; + scn->fp_to_meter = fp_to_meter; + return RES_OK; +} + +res_T +sdis_scene_get_ambient_radiative_temperature + (const struct sdis_scene* scn, + double* trad) +{ + if(!scn || !trad) return RES_BAD_ARG; + *trad = scn->ambient_radiative_temperature; + return RES_OK; +} + +res_T +sdis_scene_set_reference_temperature + (struct sdis_scene* scn, + const double tref) +{ + if(!scn || tref < 0) return RES_BAD_ARG; + scn->reference_temperature = tref; + return RES_OK; +} + +res_T +sdis_scene_get_reference_temperature + (const struct sdis_scene* scn, + double* tref) +{ + if(!scn || !tref) return RES_BAD_ARG; + *tref = scn->reference_temperature; + return RES_OK; +} + +res_T +sdis_scene_set_ambient_radiative_temperature + (struct sdis_scene* scn, + const double trad) +{ + if(!scn) return RES_BAD_ARG; + scn->ambient_radiative_temperature = trad; + return RES_OK; +} + +res_T sdis_scene_find_closest_point (const struct sdis_scene* scn, const double pos[3], @@ -424,6 +472,8 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash) } else { S3D(scene_view_primitives_count(scn->s3d_view, &nprims)); } + WRITE(&scn->reference_temperature, 1); + WRITE(&scn->fp_to_meter, 1); FOR_EACH(iprim, 0, nprims) { struct sdis_interface* interf = NULL; size_t ivert; diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -106,6 +106,20 @@ clear_properties(struct sdis_scene* scn) darray_prim_prop_clear(&scn->prim_props); } +static INLINE int +check_sdis_scene_create_args(const struct sdis_scene_create_args* args) +{ + return args + && args->get_indices + && args->get_interface + && args->get_position + && args->nprimitives + && args->nprimitives < UINT_MAX + && args->nvertices + && args->nvertices < UINT_MAX + && args->fp_to_meter > 0; +} + #endif /* SDIS_SCENE_XD_H */ #else /* !SDIS_SCENE_DIMENSION */ @@ -524,10 +538,10 @@ static res_T XD(run_analyze) (struct sdis_scene* scn, const size_t nprims, /* #primitives */ - void (*indices)(const size_t iprim, size_t ids[], void*), - void (interf)(const size_t iprim, struct sdis_interface**, void*), + sdis_get_primitive_indices_T indices, + sdis_get_primitive_interface_T interf, const size_t nverts, /* #vertices */ - void (*position)(const size_t ivert, double pos[], void*), + sdis_get_vertex_position_T position, void* ctx, struct sencXd(scene)** out_scn) { @@ -804,11 +818,9 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_scn) { struct sencXd(enclosure)* enc = NULL; unsigned ienc, nencs; - unsigned enclosed_medium; - int outer_found = 0; + int inner_multi = 0; res_T res = RES_OK; ASSERT(scn && senc3d_scn); - (void)outer_found; SENCXD(scene_get_enclosure_count(senc3d_scn, &nencs)); FOR_EACH(ienc, 0, nencs) { @@ -818,53 +830,12 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_scn) SENCXD(enclosure_get_header(enc, &header)); if(header.is_infinite) { - ASSERT(!outer_found); - outer_found = 1; + ASSERT(scn->outer_enclosure_id == UINT_MAX); /* Not set yet */ scn->outer_enclosure_id = ienc; } - /* As paths don't go in infinite enclosures we can accept models are broken - * there. But nowhere else. */ - if(header.enclosed_media_count != 1 && !header.is_infinite) { - /* Dump the problematic enclosure. */ - double tmp[DIM]; - unsigned indices[DIM]; - unsigned i; - log_warn(scn->dev, "# Found internal enclosure with %u materials:\n", - header.enclosed_media_count); - FOR_EACH(i, 0, header.enclosed_media_count) { - unsigned imed; - const struct sdis_medium* med; - SENCXD(enclosure_get_medium(enc, i, &imed)); - med = darray_medium_cdata_get(&scn->media)[imed]; - log_warn(scn->dev, "# %u (%s)\n", - imed, (med->type == SDIS_SOLID ? "solid" : "fluid")); - } - FOR_EACH(i, 0, header.vertices_count) { - SENCXD(enclosure_get_vertex(enc, i, tmp)); - #if DIM == 2 - log_warn(scn->dev, "v %g %g\n", SPLIT2(tmp)); - #else - log_warn(scn->dev, "v %g %g %g\n", SPLIT3(tmp)); - #endif - } - FOR_EACH(i, 0, header.primitives_count) { - SENCXD(enclosure_get_primitive(enc, i, indices)); - #if DIM == 2 - log_warn(scn->dev, "f %u %u\n", indices[0]+1, indices[1]+1); - #else - log_warn(scn->dev, "f %u %u %u\n", - indices[0]+1, indices[1]+1, indices[2]+1); - #endif - } - SENCXD(enclosure_ref_put(enc)); - enc = NULL; - res = RES_BAD_ARG; - goto error; - } - - SENCXD(enclosure_get_medium(enc, 0, &enclosed_medium)); - ASSERT(enclosed_medium < darray_medium_size_get(&scn->media)); + if(header.enclosed_media_count != 1 && !header.is_infinite) + inner_multi++; /* Silently discard infinite enclosures */ if(!header.is_infinite) { @@ -875,6 +846,12 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_scn) enc = NULL; } + if(inner_multi) { + log_info(scn->dev, + "# Found %d internal enclosure(s) with more than 1 medium.\n", + inner_multi); + } + /* tmp table no more useful */ htable_d_purge(&scn->tmp_hc_ub); exit: @@ -888,20 +865,14 @@ error: static res_T XD(scene_create) (struct sdis_device* dev, - const size_t nprims, /* #primitives */ - void (*indices)(const size_t iprim, size_t ids[], void*), - void (*interf)(const size_t iprim, struct sdis_interface** bound, void*), - const size_t nverts, /* #vertices */ - void (*position)(const size_t ivert, double pos[], void* ctx), - void* ctx, + const struct sdis_scene_create_args* args, struct sdis_scene** out_scn) { struct sencXd(scene)* senc3d_scn = NULL; struct sdis_scene* scn = NULL; res_T res = RES_OK; - if(!dev || !out_scn || !nprims || !indices || !interf || !nverts - || !position || nprims > UINT_MAX || nverts > UINT_MAX) { + if(!dev || !check_sdis_scene_create_args(args) || !out_scn) { res = RES_BAD_ARG; goto error; } @@ -912,10 +883,13 @@ XD(scene_create) res = RES_MEM_ERR; goto error; } + ref_init(&scn->ref); SDIS(device_ref_get(dev)); scn->dev = dev; - scn->ambient_radiative_temperature = -1; + scn->fp_to_meter = args->fp_to_meter; + scn->ambient_radiative_temperature = args->trad; + scn->reference_temperature = args->tref; scn->outer_enclosure_id = UINT_MAX; darray_interf_init(dev->allocator, &scn->interfaces); darray_medium_init(dev->allocator, &scn->media); @@ -923,12 +897,20 @@ XD(scene_create) htable_enclosure_init(dev->allocator, &scn->enclosures); htable_d_init(dev->allocator, &scn->tmp_hc_ub); - res = XD(run_analyze)(scn, nprims, indices, interf, nverts, position, ctx, &senc3d_scn); + res = XD(run_analyze) + (scn, + args->nprimitives, + args->get_indices, + args->get_interface, + args->nvertices, + args->get_position, + args->context, + &senc3d_scn); if(res != RES_OK) { log_err(dev, "%s: error during the scene analysis.\n", FUNC_NAME); goto error; } - res = XD(setup_properties)(scn, senc3d_scn, interf, ctx); + res = XD(setup_properties)(scn, senc3d_scn, args->get_interface, args->context); if(res != RES_OK) { log_err(dev, "%s: could not setup the scene interfaces and their media.\n", FUNC_NAME); diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -208,7 +208,9 @@ struct sdis_scene { struct htable_enclosure enclosures; /* Map an enclosure id to its data */ unsigned outer_enclosure_id; + double fp_to_meter; double ambient_radiative_temperature; /* In Kelvin */ + double reference_temperature; ref_T ref; struct sdis_device* dev; diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -68,9 +68,6 @@ solve_pixel struct sdis_medium* mdm, const struct sdis_camera* cam, const double time_range[2], /* Observation time */ - const double fp_to_meter, /* Scale from floating point units to meters */ - const double Tarad, /* In Kelvin */ - const double Tref, /* In Kelvin */ const size_t ipix[2], /* Pixel coordinate in the image plane */ const size_t nrealisations, const int register_paths, /* Combination of enum sdis_heat_path_flag */ @@ -82,7 +79,7 @@ solve_pixel struct sdis_heat_path* pheat_path = NULL; size_t irealisation; res_T res = RES_OK; - ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0); + ASSERT(scn && mdm && rng && cam && ipix && nrealisations); ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0); ASSERT(estimator && time_range); @@ -115,7 +112,7 @@ solve_pixel /* Launch the realisation */ res_simul = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir, - time, fp_to_meter, Tarad, Tref, pheat_path, &w); + time, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { @@ -170,9 +167,6 @@ solve_tile struct sdis_medium* mdm, const struct sdis_camera* cam, const double time_range[2], - const double fp_to_meter, - const double Tarad, - const double Tref, const size_t origin[2], /* Tile origin in image plane */ const size_t size[2], /* #pixels in X and Y */ const size_t spp, /* #samples per pixel */ @@ -183,7 +177,7 @@ solve_tile size_t mcode; /* Morton code of the tile pixel */ size_t npixels; res_T res = RES_OK; - ASSERT(scn && rng && mdm && cam && spp && origin && Tref >= 0); + ASSERT(scn && rng && mdm && cam && spp && origin); ASSERT(size &&size[0] && size[1] && buf); ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range); @@ -206,7 +200,7 @@ solve_tile /* Fetch the pixel estimator */ estimator = estimator_buffer_grab(buf, ipix[0], ipix[1]); - res = solve_pixel(scn, rng, mdm, cam, time_range, fp_to_meter, Tarad, Tref, + res = solve_pixel(scn, rng, mdm, cam, time_range, ipix, spp, register_paths, pix_sz, estimator); if(res != RES_OK) goto error; } @@ -362,12 +356,9 @@ sdis_solve_camera || !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 @@ -448,9 +439,7 @@ sdis_solve_camera /* Draw the tile */ 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); + tile_org, tile_sz, args->spp, args->register_paths, pix_sz, buf); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -102,8 +102,7 @@ XD(solve_boundary) ATOMIC res = RES_OK; if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX - || !args->primitives || !args->sides || !args->nprimitives - || args->fp_to_meter <= 0) { + || !args->primitives || !args->sides || !args->nprimitives) { res = RES_BAD_ARG; goto error; } @@ -111,14 +110,14 @@ XD(solve_boundary) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - 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(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { + res = RES_BAD_ARG; + goto error; + } + if(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 @@ -256,16 +255,8 @@ XD(solve_boundary) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported yet */ - time = INF; + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -274,6 +265,11 @@ XD(solve_boundary) pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + /* Sample a position onto the boundary */ #if SDIS_XD_DIMENSION == 2 res_local = s2d_scene_view_sample @@ -303,8 +299,7 @@ XD(solve_boundary) /* Invoke the boundary realisation */ res_simul = XD(boundary_realisation)(scn, rng, iprim, uv, time, side, - args->fp_to_meter, args->ambient_radiative_temperature, - args->reference_temperature, pgreen_path, pheat_path, &w); + pgreen_path, pheat_path, &w); /* Fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { @@ -464,7 +459,6 @@ XD(solve_boundary_flux) || 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; @@ -579,8 +573,7 @@ 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; + const double Tref = scn->reference_temperature; double epsilon, hc, hr, imposed_flux, imposed_temp; size_t iprim; double uv[DIM - 1]; @@ -663,7 +656,7 @@ XD(solve_boundary_flux) 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, args->fp_to_meter, Tarad, Tref, flux_mask, T_brf); + solid_side, flux_mask, T_brf); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h @@ -222,7 +222,7 @@ XD(solve_medium) ATOMIC res = RES_OK; if(!scn || !args || !args->medium || !args->nrealisations - || args->nrealisations > INT64_MAX || args->fp_to_meter <= 0) { + || args->nrealisations > INT64_MAX) { res = RES_BAD_ARG; goto error; } @@ -321,29 +321,22 @@ XD(solve_medium) if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ time_current(&t0); - - if(!out_green) { - /* Sample the time */ - time = sample_time(rng, args->time_range); - - /* Prepare path registration if necessary */ - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported yet */ - time = INF; + + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); goto error_it; } - pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + /* Uniformly Sample an enclosure that surround the submitted medium and * uniformly sample a position into it */ enc = sample_medium_enclosure(&cumul, rng); @@ -356,9 +349,7 @@ XD(solve_medium) /* Run a probe realisation */ 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); + args->medium, pos, time, pgreen_path, pheat_path, &weight); if(res_simul != RES_OK && res_simul != RES_BAD_OP) { ATOMIC_SET(&res, res_simul); @@ -504,7 +495,6 @@ XD(compute_power) || !args->medium || !args->nrealisations || args->nrealisations > INT64_MAX - || args->fp_to_meter <= 0 || args->time_range[0] < 0 || args->time_range[0] > args->time_range[1] || ( args->time_range[1] > DBL_MAX diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h @@ -53,8 +53,7 @@ XD(solve_probe) ATOMIC nsolved_realisations = 0; ATOMIC res = RES_OK; - if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX - || args->fp_to_meter <= 0) { + if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX) { res = RES_BAD_ARG; goto error; } @@ -62,16 +61,15 @@ XD(solve_probe) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - 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(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { + res = RES_BAD_ARG; + goto error; + } + if(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; } @@ -153,16 +151,9 @@ XD(solve_probe) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported */ - time = INF; + time = sample_time(rng, args->time_range); + + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -170,11 +161,13 @@ XD(solve_probe) } pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } res_simul = XD(probe_realisation)((size_t)irealisation, scn, rng, medium, - args->position, time, args->fp_to_meter, - args->ambient_radiative_temperature, args->reference_temperature, - pgreen_path, pheat_path, &w); + args->position, time, pgreen_path, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -53,7 +53,7 @@ XD(solve_probe_boundary) ATOMIC res = RES_OK; if(!scn || !args || !args->nrealisations || args->nrealisations > INT64_MAX - || args->fp_to_meter <= 0 || ((unsigned)args->side >= SDIS_SIDE_NULL__)) { + || ((unsigned)args->side >= SDIS_SIDE_NULL__)) { res = RES_BAD_ARG; goto error; } @@ -61,14 +61,14 @@ XD(solve_probe_boundary) res = RES_BAD_ARG; goto error; } - if(out_estimator) { - 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(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { + res = RES_BAD_ARG; + goto error; + } + if(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 @@ -188,16 +188,8 @@ XD(solve_probe_boundary) /* Begin time registration */ time_current(&t0); - if(!out_green) { - time = sample_time(rng, args->time_range); - if(register_paths) { - heat_path_init(scn->dev->allocator, &heat_path); - pheat_path = &heat_path; - } - } else { - /* Do not take care of the submitted time when registering the green - * function. Only steady systems are supported */ - time = INF; + time = sample_time(rng, args->time_range); + if(out_green) { res_local = green_function_create_path(greens[ithread], &green_path); if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); @@ -206,9 +198,13 @@ XD(solve_probe_boundary) pgreen_path = &green_path; } + if(register_paths) { + heat_path_init(scn->dev->allocator, &heat_path); + pheat_path = &heat_path; + } + 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); + args->side, pgreen_path, pheat_path, &w); /* Handle fatal error */ if(res_simul != RES_OK && res_simul != RES_BAD_OP) { @@ -357,7 +353,7 @@ XD(solve_probe_boundary_flux) 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) { + || !out_estimator) { res = RES_BAD_ARG; goto error; } @@ -484,8 +480,7 @@ XD(solve_probe_boundary_flux) double time, epsilon, hc, hr, imposed_flux, imposed_temp; int flux_mask = 0; double T_brf[3] = { 0, 0, 0 }; - const double Tref = args->reference_temperature; - const double Tarad = args->ambient_radiative_temperature; + const double Tref = scn->reference_temperature; size_t n; int pcent; res_T res_simul = RES_OK; @@ -519,7 +514,7 @@ XD(solve_probe_boundary_flux) if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE; if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE; 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); + time, solid_side, flux_mask, T_brf); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); diff --git a/src/test_sdis_compute_mean_power.c b/src/test_sdis_compute_mean_power.c @@ -1,343 +0,0 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. */ - -#include "sdis.h" -#include "test_sdis_utils.h" - -#include <rsys/stretchy_array.h> -#include <star/s3dut.h> - -#define UNKOWN_TEMPERATURE -1 -#define N 100000ul /* #realisations */ -#define POWER0 10 -#define POWER1 5 - -static INLINE void -check_intersection - (const double val0, - const double eps0, - const double val1, - const double eps1) -{ - double interval0[2], interval1[2]; - double intersection[2]; - interval0[0] = val0 - eps0; - interval0[1] = val0 + eps0; - interval1[0] = val1 - eps1; - interval1[1] = val1 + eps1; - intersection[0] = MMAX(interval0[0], interval1[0]); - intersection[1] = MMIN(interval0[1], interval1[1]); - CHK(intersection[0] <= intersection[1]); -} - -/******************************************************************************* - * Geometry - ******************************************************************************/ -struct context { - struct s3dut_mesh_data msh0; - struct s3dut_mesh_data msh1; - struct sdis_interface* interf0; - struct sdis_interface* interf1; -}; - -static void -get_indices(const size_t itri, size_t ids[3], void* context) -{ - const struct context* ctx = context; - /* Note that we swap the indices to ensure that the triangle normals point - * inward the super shape */ - if(itri < ctx->msh0.nprimitives) { - ids[0] = ctx->msh0.indices[itri*3+0]; - ids[2] = ctx->msh0.indices[itri*3+1]; - ids[1] = ctx->msh0.indices[itri*3+2]; - } else { - const size_t itri2 = itri - ctx->msh0.nprimitives; - ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices; - ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices; - ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices; - } -} - -static void -get_position(const size_t ivert, double pos[3], void* context) -{ - const struct context* ctx = context; - if(ivert < ctx->msh0.nvertices) { - pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0; - pos[1] = ctx->msh0.positions[ivert*3+1]; - pos[2] = ctx->msh0.positions[ivert*3+2]; - } else { - const size_t ivert2 = ivert - ctx->msh0.nvertices; - pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0; - pos[1] = ctx->msh1.positions[ivert2*3+1]; - pos[2] = ctx->msh1.positions[ivert2*3+2]; - } -} - -static void -get_interface(const size_t itri, struct sdis_interface** bound, void* context) -{ - const struct context* ctx = context; - *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1; -} - -/******************************************************************************* - * Interface - ******************************************************************************/ -static double -interface_get_convection_coef - (const struct sdis_interface_fragment* frag, struct sdis_data* data) -{ - CHK(frag != NULL); (void)data; - return 1.0; -} - -/******************************************************************************* - * Fluid medium - ******************************************************************************/ -static double -fluid_get_temperature - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx == NULL); (void)data; - return 300; -} - -/******************************************************************************* - * Solid medium - ******************************************************************************/ -static double -solid_get_calorific_capacity - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx != NULL); (void)data; - return 1.0; -} - -static double -solid_get_thermal_conductivity - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx != NULL); (void)data; - return 1.0; -} - -static double -solid_get_volumic_mass - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx != NULL); (void)data; - return 1.0; -} - -static double -solid_get_delta - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx != NULL); (void)data; - return 1.0 / 20.0; -} - -static double -solid_get_volumic_power - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - CHK(vtx != NULL); (void)data; - return vtx->P[0] < 0 ? POWER0 : POWER1; -} - -/******************************************************************************* - * Test - ******************************************************************************/ -int -main(int argc, char** argv) -{ - struct mem_allocator allocator; - struct context ctx; - struct s3dut_mesh* sphere = NULL; - struct s3dut_mesh* cylinder = NULL; - struct sdis_data* data = NULL; - struct sdis_device* dev = NULL; - struct sdis_estimator* estimator = NULL; - struct sdis_medium* fluid = NULL; - struct sdis_medium* solid0 = NULL; - struct sdis_medium* solid1 = NULL; - struct sdis_interface* interf0 = NULL; - struct sdis_interface* interf1 = NULL; - struct sdis_scene* scn = NULL; - struct sdis_mc mpow = SDIS_MC_NULL; - struct sdis_mc time = SDIS_MC_NULL; - struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; - struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; - struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; - struct sdis_compute_power_args args = SDIS_COMPUTE_POWER_ARGS_DEFAULT; - size_t nverts = 0; - size_t ntris = 0; - double ref = 0; - (void)argc, (void) argv; - - OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); - OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev)); - - /* Setup the interface shader */ - interf_shader.convection_coef = interface_get_convection_coef; - - /* Setup the fluid shader */ - fluid_shader.temperature = fluid_get_temperature; - - /* Setup the solid shader */ - solid_shader.calorific_capacity = solid_get_calorific_capacity; - solid_shader.thermal_conductivity = solid_get_thermal_conductivity; - solid_shader.volumic_mass = solid_get_volumic_mass; - solid_shader.delta_solid = solid_get_delta; - solid_shader.volumic_power = solid_get_volumic_power; - - /* Create the fluid */ - OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); - - /* Create the solids */ - OK(sdis_solid_create(dev, &solid_shader, data, &solid0)); - OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); - - /* Create the interface0 */ - OK(sdis_interface_create(dev, solid0, fluid, &interf_shader, NULL, &interf0)); - OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, NULL, &interf1)); - ctx.interf0 = interf0; - ctx.interf1 = interf1; - - /* Create the geometry */ - OK(s3dut_create_sphere(&allocator, 1, 512, 256, &sphere)); - OK(s3dut_create_cylinder(&allocator, 1, 10, 512, 8, &cylinder)); - OK(s3dut_mesh_get_data(sphere, &ctx.msh0)); - OK(s3dut_mesh_get_data(cylinder, &ctx.msh1)); - - /* Create the scene */ - ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives; - nverts = ctx.msh0.nvertices + ctx.msh1.nvertices; - OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); - - /* Test sdis_compute_power function */ - args.nrealisations = N; - args.medium = solid0; - args.time_range[0] = INF; - args.time_range[1] = INF; - BA(sdis_compute_power(NULL, &args, &estimator)); - BA(sdis_compute_power(scn, NULL, &estimator)); - BA(sdis_compute_power(scn, &args, NULL)); - args.nrealisations = 0; - BA(sdis_compute_power(scn, &args, &estimator)); - args.nrealisations = N; - args.medium = NULL; - BA(sdis_compute_power(scn, &args, &estimator)); - args.medium = solid0; - args.fp_to_meter = 0; - BA(sdis_compute_power(scn, &args, &estimator)); - args.fp_to_meter = 1; - args.time_range[0] = args.time_range[1] = -1; - BA(sdis_compute_power(scn, &args, &estimator)); - args.time_range[0] = 1; - BA(sdis_compute_power(scn, &args, &estimator)); - args.time_range[1] = 0; - BA(sdis_compute_power(scn, &args, &estimator)); - args.time_range[0] = args.time_range[1] = INF; - OK(sdis_compute_power(scn, &args, &estimator)); - - BA(sdis_estimator_get_power(NULL, &mpow)); - BA(sdis_estimator_get_power(estimator, NULL)); - OK(sdis_estimator_get_power(estimator, &mpow)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - /* Check results for solid 0 */ - ref = 4.0/3.0 * PI * POWER0; - printf("Mean power of the solid0 = %g ~ %g +/- %g\n", - ref, mpow.E, mpow.SE); - check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); - OK(sdis_estimator_ref_put(estimator)); - - /* Check results for solid 1 */ - args.medium = solid1; - OK(sdis_compute_power(scn, &args, &estimator)); - OK(sdis_estimator_get_power(estimator, &mpow)); - ref = PI * 10 * POWER1; - printf("Mean power of the solid1 = %g ~ %g +/- %g\n", - ref, mpow.E, mpow.SE); - check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); - OK(sdis_estimator_ref_put(estimator)); - - /* Check for a not null time range */ - args.time_range[0] = 0; - args.time_range[1] = 10; - OK(sdis_compute_power(scn, &args, &estimator)); - OK(sdis_estimator_get_power(estimator, &mpow)); - ref = PI * 10 * POWER1 / 10; - printf("Mean power of the solid1 in [0, 10] s = %g ~ %g +/- %g\n", - ref, mpow.E, mpow.SE); - check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); - OK(sdis_estimator_ref_put(estimator)); - - /* Reset the scene with only one solid medium */ - OK(sdis_scene_ref_put(scn)); - ctx.interf0 = interf0; - ctx.interf1 = interf0; - OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); - - /* Check invalid medium */ - args.time_range[0] = args.time_range[1] = 1; - args.medium = solid1; - BA(sdis_compute_power(scn, &args, &estimator)); - - /* Check non constant volumic power */ - args.medium = solid0; - OK(sdis_compute_power(scn, &args, &estimator)); - OK(sdis_estimator_get_power(estimator, &mpow)); - ref = 4.0/3.0*PI*POWER0 + PI*10*POWER1; - printf("Mean power of the sphere+cylinder = %g ~ %g +/- %g\n", - ref, mpow.E, mpow.SE); - check_intersection(ref, 1.5e-1, mpow.E, 3*mpow.SE); - OK(sdis_estimator_ref_put(estimator)); - -#if 0 - { - double* vertices = NULL; - size_t* indices = NULL; - size_t i; - CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices))); - CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices))); - FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx); - FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx); - dump_mesh(stdout, vertices, nverts, indices, ntris); - MEM_RM(&allocator, vertices); - MEM_RM(&allocator, indices); - } -#endif - - /* Clean up memory */ - OK(sdis_device_ref_put(dev)); - OK(sdis_medium_ref_put(fluid)); - OK(sdis_medium_ref_put(solid0)); - OK(sdis_medium_ref_put(solid1)); - OK(sdis_interface_ref_put(interf0)); - OK(sdis_interface_ref_put(interf1)); - OK(sdis_scene_ref_put(scn)); - OK(s3dut_mesh_ref_put(sphere)); - OK(s3dut_mesh_ref_put(cylinder)); - - check_memory_allocator(&allocator); - mem_shutdown_proxy_allocator(&allocator); - CHK(mem_allocated_size() == 0); - return 0; -} diff --git a/src/test_sdis_compute_power.c b/src/test_sdis_compute_power.c @@ -0,0 +1,345 @@ +/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sdis.h" +#include "test_sdis_utils.h" + +#include <rsys/stretchy_array.h> +#include <star/s3dut.h> + +#define UNKOWN_TEMPERATURE -1 +#define N 100000ul /* #realisations */ +#define POWER0 10 +#define POWER1 5 + +static INLINE void +check_intersection + (const double val0, + const double eps0, + const double val1, + const double eps1) +{ + double interval0[2], interval1[2]; + double intersection[2]; + interval0[0] = val0 - eps0; + interval0[1] = val0 + eps0; + interval1[0] = val1 - eps1; + interval1[1] = val1 + eps1; + intersection[0] = MMAX(interval0[0], interval1[0]); + intersection[1] = MMIN(interval0[1], interval1[1]); + CHK(intersection[0] <= intersection[1]); +} + +/******************************************************************************* + * Geometry + ******************************************************************************/ +struct context { + struct s3dut_mesh_data msh0; + struct s3dut_mesh_data msh1; + struct sdis_interface* interf0; + struct sdis_interface* interf1; +}; + +static void +get_indices(const size_t itri, size_t ids[3], void* context) +{ + const struct context* ctx = context; + /* Note that we swap the indices to ensure that the triangle normals point + * inward the super shape */ + if(itri < ctx->msh0.nprimitives) { + ids[0] = ctx->msh0.indices[itri*3+0]; + ids[2] = ctx->msh0.indices[itri*3+1]; + ids[1] = ctx->msh0.indices[itri*3+2]; + } else { + const size_t itri2 = itri - ctx->msh0.nprimitives; + ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices; + ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices; + ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices; + } +} + +static void +get_position(const size_t ivert, double pos[3], void* context) +{ + const struct context* ctx = context; + if(ivert < ctx->msh0.nvertices) { + pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0; + pos[1] = ctx->msh0.positions[ivert*3+1]; + pos[2] = ctx->msh0.positions[ivert*3+2]; + } else { + const size_t ivert2 = ivert - ctx->msh0.nvertices; + pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0; + pos[1] = ctx->msh1.positions[ivert2*3+1]; + pos[2] = ctx->msh1.positions[ivert2*3+2]; + } +} + +static void +get_interface(const size_t itri, struct sdis_interface** bound, void* context) +{ + const struct context* ctx = context; + *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1; +} + +/******************************************************************************* + * Interface + ******************************************************************************/ +static double +interface_get_convection_coef + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + CHK(frag != NULL); (void)data; + return 1.0; +} + +/******************************************************************************* + * Fluid medium + ******************************************************************************/ +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx == NULL); (void)data; + return 300; +} + +/******************************************************************************* + * Solid medium + ******************************************************************************/ +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); (void)data; + return 1.0; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); (void)data; + return 1.0; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); (void)data; + return 1.0; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); (void)data; + return 1.0 / 20.0; +} + +static double +solid_get_volumic_power + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx != NULL); (void)data; + return vtx->P[0] < 0 ? POWER0 : POWER1; +} + +/******************************************************************************* + * Test + ******************************************************************************/ +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct context ctx; + struct s3dut_mesh* sphere = NULL; + struct s3dut_mesh* cylinder = NULL; + struct sdis_data* data = NULL; + struct sdis_device* dev = NULL; + struct sdis_estimator* estimator = NULL; + struct sdis_medium* fluid = NULL; + struct sdis_medium* solid0 = NULL; + struct sdis_medium* solid1 = NULL; + struct sdis_interface* interf0 = NULL; + struct sdis_interface* interf1 = NULL; + struct sdis_scene* scn = NULL; + struct sdis_mc mpow = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; + struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; + struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; + struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_compute_power_args args = SDIS_COMPUTE_POWER_ARGS_DEFAULT; + size_t nverts = 0; + size_t ntris = 0; + double ref = 0; + (void)argc, (void) argv; + + OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); + OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev)); + + /* Setup the interface shader */ + interf_shader.convection_coef = interface_get_convection_coef; + + /* Setup the fluid shader */ + fluid_shader.temperature = fluid_get_temperature; + + /* Setup the solid shader */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.volumic_power = solid_get_volumic_power; + + /* Create the fluid */ + OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + + /* Create the solids */ + OK(sdis_solid_create(dev, &solid_shader, data, &solid0)); + OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); + + /* Create the interface0 */ + OK(sdis_interface_create(dev, solid0, fluid, &interf_shader, NULL, &interf0)); + OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, NULL, &interf1)); + ctx.interf0 = interf0; + ctx.interf1 = interf1; + + /* Create the geometry */ + OK(s3dut_create_sphere(&allocator, 1, 512, 256, &sphere)); + OK(s3dut_create_cylinder(&allocator, 1, 10, 512, 8, &cylinder)); + OK(s3dut_mesh_get_data(sphere, &ctx.msh0)); + OK(s3dut_mesh_get_data(cylinder, &ctx.msh1)); + + /* Create the scene */ + ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives; + nverts = ctx.msh0.nvertices + ctx.msh1.nvertices; + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ntris; + scn_args.nvertices = nverts; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); + + /* Test sdis_compute_power function */ + args.nrealisations = N; + args.medium = solid0; + args.time_range[0] = INF; + args.time_range[1] = INF; + BA(sdis_compute_power(NULL, &args, &estimator)); + BA(sdis_compute_power(scn, NULL, &estimator)); + BA(sdis_compute_power(scn, &args, NULL)); + args.nrealisations = 0; + BA(sdis_compute_power(scn, &args, &estimator)); + args.nrealisations = N; + args.medium = NULL; + BA(sdis_compute_power(scn, &args, &estimator)); + args.medium = solid0; + args.time_range[0] = args.time_range[1] = -1; + BA(sdis_compute_power(scn, &args, &estimator)); + args.time_range[0] = 1; + BA(sdis_compute_power(scn, &args, &estimator)); + args.time_range[1] = 0; + BA(sdis_compute_power(scn, &args, &estimator)); + args.time_range[0] = args.time_range[1] = INF; + OK(sdis_compute_power(scn, &args, &estimator)); + + BA(sdis_estimator_get_power(NULL, &mpow)); + BA(sdis_estimator_get_power(estimator, NULL)); + OK(sdis_estimator_get_power(estimator, &mpow)); + OK(sdis_estimator_get_realisation_time(estimator, &time)); + + /* Check results for solid 0 */ + ref = 4.0/3.0 * PI * POWER0; + printf("Mean power of the solid0 = %g ~ %g +/- %g\n", + ref, mpow.E, mpow.SE); + check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); + OK(sdis_estimator_ref_put(estimator)); + + /* Check results for solid 1 */ + args.medium = solid1; + OK(sdis_compute_power(scn, &args, &estimator)); + OK(sdis_estimator_get_power(estimator, &mpow)); + ref = PI * 10 * POWER1; + printf("Mean power of the solid1 = %g ~ %g +/- %g\n", + ref, mpow.E, mpow.SE); + check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); + OK(sdis_estimator_ref_put(estimator)); + + /* Check for a not null time range */ + args.time_range[0] = 0; + args.time_range[1] = 10; + OK(sdis_compute_power(scn, &args, &estimator)); + OK(sdis_estimator_get_power(estimator, &mpow)); + ref = PI * 10 * POWER1 / 10; + printf("Mean power of the solid1 in [0, 10] s = %g ~ %g +/- %g\n", + ref, mpow.E, mpow.SE); + check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE); + OK(sdis_estimator_ref_put(estimator)); + + /* Reset the scene with only one solid medium */ + OK(sdis_scene_ref_put(scn)); + ctx.interf0 = interf0; + ctx.interf1 = interf0; + OK(sdis_scene_create(dev, &scn_args, &scn)); + + /* Check invalid medium */ + args.time_range[0] = args.time_range[1] = 1; + args.medium = solid1; + BA(sdis_compute_power(scn, &args, &estimator)); + + /* Check non constant volumic power */ + args.medium = solid0; + OK(sdis_compute_power(scn, &args, &estimator)); + OK(sdis_estimator_get_power(estimator, &mpow)); + ref = 4.0/3.0*PI*POWER0 + PI*10*POWER1; + printf("Mean power of the sphere+cylinder = %g ~ %g +/- %g\n", + ref, mpow.E, mpow.SE); + check_intersection(ref, 1.5e-1, mpow.E, 3*mpow.SE); + OK(sdis_estimator_ref_put(estimator)); + +#if 0 + { + double* vertices = NULL; + size_t* indices = NULL; + size_t i; + CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices))); + CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices))); + FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx); + FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx); + dump_mesh(stdout, vertices, nverts, indices, ntris); + MEM_RM(&allocator, vertices); + MEM_RM(&allocator, indices); + } +#endif + + /* Clean up memory */ + OK(sdis_device_ref_put(dev)); + OK(sdis_medium_ref_put(fluid)); + OK(sdis_medium_ref_put(solid0)); + OK(sdis_medium_ref_put(solid1)); + OK(sdis_interface_ref_put(interf0)); + OK(sdis_interface_ref_put(interf1)); + OK(sdis_scene_ref_put(scn)); + OK(s3dut_mesh_ref_put(sphere)); + OK(s3dut_mesh_ref_put(cylinder)); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHK(mem_allocated_size() == 0); + return 0; +} diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -126,6 +126,8 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* ctx) ******************************************************************************/ struct solid { double lambda; + double initial_temperature; + double t0; }; static double @@ -168,6 +170,21 @@ solid_get_delta return 1.0/10.0; } +static double +solid_get_temperature +(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + double t0; + CHK(vtx != NULL); + CHK(data != NULL); + t0 = ((const struct solid*)sdis_data_cget(data))->t0; + if(vtx->time > t0) { + return UNKNOWN_TEMPERATURE; + } else { + return ((const struct solid*)sdis_data_cget(data))->initial_temperature; + } +} + /******************************************************************************* * Interface ******************************************************************************/ @@ -265,12 +282,13 @@ main(int argc, char** argv) struct sdis_medium* solid2 = NULL; struct sdis_interface* interfaces[5] = {NULL}; struct sdis_interface* prim_interfaces[32/*#triangles*/]; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_scene* scn = NULL; struct ssp_rng* rng = NULL; - const size_t nsimuls = 4; - size_t isimul; + const int nsimuls = 4; + int isimul; const double emissivity = 1;/* Emissivity of the side +/-X of the solid */ const double lambda = 0.1; /* Conductivity of the solid */ const double Tref = 300; /* Reference temperature */ @@ -292,11 +310,13 @@ main(int argc, char** argv) OK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); ((struct solid*)sdis_data_get(data))->lambda = lambda; + ((struct solid*)sdis_data_get(data))->t0 = 0; + ((struct solid*)sdis_data_get(data))->initial_temperature = (T0 + T1) / 2; solid_shader.calorific_capacity = solid_get_calorific_capacity; solid_shader.thermal_conductivity = solid_get_thermal_conductivity; solid_shader.volumic_mass = solid_get_volumic_mass; solid_shader.delta_solid = solid_get_delta; - solid_shader.temperature = temperature_unknown; + solid_shader.temperature = solid_get_temperature; OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); @@ -373,8 +393,14 @@ main(int argc, char** argv) geom.positions = vertices; geom.indices = indices; geom.interfaces = prim_interfaces; - OK(sdis_scene_create(dev, ntriangles, get_indices, get_interface, nvertices, - get_position, &geom, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ntriangles; + scn_args.nvertices = nvertices; + scn_args.tref = Tref; + scn_args.context = &geom; + OK(sdis_scene_create(dev, &scn_args, &scn)); hr = 4.0 * BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity; @@ -389,23 +415,26 @@ main(int argc, char** argv) struct sdis_estimator* estimator2; struct sdis_green_function* green; struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; - double ref, u; + double ref = -1; size_t nreals = 0; size_t nfails = 0; const size_t N = 10000; double T1b; + int steady = (isimul % 2) == 0; /* Reset temperature */ p_intface->temperature = T1; solve_args.nrealisations = N; - solve_args.time_range[0] = INF; - solve_args.time_range[1] = INF; + if(steady) { + solve_args.time_range[0] = solve_args.time_range[1] = INF; + } else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } 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); - u = (solve_args.position[0] + 1) / thickness; - solve_args.reference_temperature = Tref; OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -416,28 +445,36 @@ main(int argc, char** argv) tmp = lambda / (2 * lambda + thickness * hr) * (T1 - T0); Ts0 = T0 + tmp; Ts1 = T1 - tmp; - ref = u * Ts1 + (1-u) * Ts0; - printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + if(steady) { + double u = (solve_args.position[0] + 1) / thickness; + ref = u * Ts1 + (1 - u) * Ts0; + printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" + " and T1=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + p_intface->temperature, T.E, T.SE); + } printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); printf("\n"); - /* Check green used at a different temperature */ + /* Check same green used at a different temperature */ p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10; OK(sdis_solve_probe(scn, &solve_args, &estimator)); @@ -449,17 +486,26 @@ main(int argc, char** argv) tmp = lambda / (2 * lambda + thickness * hr) * (T1b - T0); Ts0 = T0 + tmp; Ts1 = T1b - tmp; - ref = u * Ts1 + (1 - u) * Ts0; - printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", - SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + + if(steady) { + double u = (solve_args.position[0] + 1) / thickness; + ref = u * Ts1 + (1 - u) * Ts0; + printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" + " and T1=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + p_intface->temperature, T.E, T.SE); + } printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, 3*T.SE) == 1); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -473,7 +519,7 @@ main(int argc, char** argv) OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); - printf("\n"); + printf("\n\n"); } /* Release memory */ diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -283,6 +283,7 @@ main(int argc, char** argv) struct sdis_medium* solid2 = NULL; struct sdis_interface* interfaces[5] = {NULL}; struct sdis_interface* prim_interfaces[10/*#segment*/]; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; const size_t nsimuls = 4; @@ -378,8 +379,14 @@ main(int argc, char** argv) geom.positions = vertices; geom.indices = indices; geom.interfaces = prim_interfaces; - OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, nvertices, - get_position, &geom, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = nsegments; + scn_args.nvertices = nvertices; + scn_args.tref = Tref; + scn_args.context = &geom; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity; tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0); @@ -405,7 +412,6 @@ main(int argc, char** argv) 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, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); @@ -426,10 +432,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); @@ -441,7 +447,7 @@ main(int argc, char** argv) OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_ref_put(estimator)); - printf("\n"); + printf("\n\n"); } /* Release memory */ diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -183,6 +183,7 @@ main(int argc, char** argv) struct sdis_estimator* estimator = NULL; struct sdis_estimator* estimator2 = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; @@ -245,14 +246,22 @@ main(int argc, char** argv) square_interfaces[3] = interf_T1; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_T0)); @@ -297,10 +306,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, box_scn, solve_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } @@ -339,10 +348,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, square_scn, solve_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -193,6 +193,7 @@ main(int argc, char** argv) struct sdis_estimator* estimator = NULL; struct sdis_estimator* estimator2 = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; @@ -261,14 +262,22 @@ main(int argc, char** argv) square_interfaces[3] = interf_T1; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_T0)); @@ -313,10 +322,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, box_scn, solve_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } @@ -354,10 +363,10 @@ main(int argc, char** argv) if(IS_INF(time)) { /* Check green function */ OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, square_scn, solve_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator2)); OK(sdis_green_function_ref_put(green)); } diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -18,13 +18,14 @@ #include <rsys/clock_time.h> #include <rsys/double3.h> +#include <star/ssp.h> /* * The scene is composed of a solid cube/square whose temperature is unknown. * The temperature is fixed at T0 on the +X face. The Flux of the -X face is * fixed to PHI. The flux on the other faces is null (i.e. adiabatic). This * test computes the temperature of a probe position pos into the solid and - * check that it is equal to: + * check that at t=inf it is equal to: * * T(pos) = T0 + (A-pos) * PHI/LAMBDA * @@ -94,7 +95,10 @@ solid_get_temperature { (void)data; CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; + if(vtx->time > 0) + return UNKNOWN_TEMPERATURE; + else + return T0; } /******************************************************************************* @@ -129,7 +133,7 @@ interface_get_flux static void solve (struct sdis_scene* scn, - const double pos[], + struct ssp_rng* rng, struct interf* interf) { char dump[128]; @@ -142,139 +146,195 @@ solve struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals; size_t nfails; - double ref; - const double time_range[2] = {INF, INF}; + double ref = -1; + const int nsimuls = 4; + int isimul; enum sdis_scene_dimension dim; - ASSERT(scn && pos && interf); - - /* Restore phi value */ - interf->phi = PHI; - - ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; + ASSERT(scn && rng && interf); 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, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - 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)); - - switch(dim) { + FOR_EACH(isimul, 0, nsimuls) { + int steady = (isimul % 2) == 0; + + /* Restore phi value */ + interf->phi = PHI; + + solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[2] = + dim == SDIS_SCENE_2D ? 0 : ssp_rng_uniform_double(rng, 0.1, 0.9); + + solve_args.nrealisations = N; + if(steady) + solve_args.time_range[0] = solve_args.time_range[1] = INF; + else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + 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)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, T.SE*3)); - - time_current(&t0); - 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); - - OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); - OK(sdis_estimator_get_failure_count(estimator2, &nfails)); - OK(sdis_estimator_get_temperature(estimator2, &T)); - - switch(dim) { + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); + + time_current(&t0); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + time_current(&t1); + OK(sdis_green_function_solve(green, &estimator2)); + time_current(&t2); + + OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); + OK(sdis_estimator_get_failure_count(estimator2, &nfails)); + OK(sdis_estimator_get_temperature(estimator2, &T)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady Green temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean Green temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady Green temperature at (%g, %g, %g) with Phi=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean Green temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - time_sub(&t0, &t1, &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green estimation time = %s\n", dump); - time_sub(&t1, &t2, &t1); - time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); - - check_green_function(green); - check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, time_range); - - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - printf("\n"); - - /* Check green used at a different phi */ - interf->phi = 3 * PHI; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA; - - switch (dim) { - case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), interf->phi, ref, T.E, T.SE); - break; - case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), interf->phi, ref, T.E, T.SE); - break; - default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + time_sub(&t0, &t1, &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green estimation time = %s\n", dump); + time_sub(&t1, &t2, &t1); + time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + check_green_serialization(green, scn); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check same green used at a different flux value */ + interf->phi = 3 * PHI; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + 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)); + + switch(dim) { + case SDIS_SCENE_2D: + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g) with Phi=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } + break; + case SDIS_SCENE_3D: + if(steady) { + ref = T0 + (1 - solve_args.position[0]) * interf->phi / LAMBDA; + printf("Steady temperature at (%g, %g, %g) with Phi=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), interf->phi, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Phi=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + interf->phi, T.E, T.SE); + } + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, T.SE * 3)); + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); - time_current(&t0); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); + time_current(&t0); + OK(sdis_green_function_solve(green, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); - check_green_function(green); - check_estimator_eq(estimator, estimator2); + check_green_function(green); + check_estimator_eq(estimator, estimator2); - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - OK(sdis_green_function_ref_put(green)); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); - printf("\n"); + printf("\n\n"); + } } /******************************************************************************* @@ -293,13 +353,14 @@ main(int argc, char** argv) struct sdis_interface* interf_phi = NULL; struct sdis_scene* box_scn = NULL; struct sdis_scene* square_scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; - double pos[3]; + struct ssp_rng* rng = NULL; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -365,14 +426,22 @@ main(int argc, char** argv) square_interfaces[3] = interf_T0; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_adiabatic)); @@ -380,15 +449,16 @@ main(int argc, char** argv) OK(sdis_interface_ref_put(interf_phi)); /* Solve */ - d3_splat(pos, 0.25); + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); printf(">> Box scene\n"); - solve(box_scn, pos, interf_props); + solve(box_scn, rng, interf_props); printf(">> Square Scene\n"); - solve(square_scn, pos, interf_props); + solve(square_scn, rng, interf_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(ssp_rng_ref_put(rng)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c @@ -97,6 +97,7 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf) struct context ctx; struct senc2d_scene* scn2d; struct senc3d_scene* scn3d; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_scene* scn = NULL; size_t ntris, npos; size_t iprim; @@ -111,37 +112,51 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf) ntris = box_ntriangles; npos = box_nvertices; - #define CREATE sdis_scene_create - #define IDS get_indices_3d - #define POS get_position_3d - #define IFA get_interface - - BA(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL)); - BA(CREATE(dev, 0, IDS, IFA, npos, POS, &ctx, &scn)); - BA(CREATE(dev, ntris, NULL, IFA, npos, POS, &ctx, &scn)); - BA(CREATE(dev, ntris, IDS, NULL, npos, POS, &ctx, &scn)); - BA(CREATE(dev, ntris, IDS, IFA, 0, POS, &ctx, &scn)); - BA(CREATE(dev, ntris, IDS, IFA, npos, NULL, &ctx, &scn)); + scn_args.get_indices = get_indices_3d; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position_3d; + scn_args.nprimitives = ntris; + scn_args.nvertices = npos; + scn_args.context = &ctx; + BA(sdis_scene_create(NULL, &SDIS_SCENE_CREATE_ARGS_DEFAULT, NULL)); + BA(sdis_scene_create(NULL, &SDIS_SCENE_CREATE_ARGS_DEFAULT, NULL)); + + scn_args.nprimitives = 0; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.nprimitives = ntris; + scn_args.get_indices = NULL; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.get_indices = get_indices_3d; + scn_args.get_interface = NULL; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.get_interface = get_interface; + scn_args.get_position = NULL; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.get_position = get_position_3d; + scn_args.nvertices = 0; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.nvertices = npos; + scn_args.fp_to_meter = 0; + BA(sdis_scene_create(dev, &scn_args, &scn)); + scn_args.fp_to_meter = 1; /* Duplicated vertex */ ctx.positions = duplicated_vertices; ctx.indices = dup_vrtx_indices; - BA(CREATE(dev, 1, IDS, IFA, npos, POS, &ctx, &scn)); + BA(sdis_scene_create(dev, &scn_args, &scn)); /* Duplicated triangle */ ctx.positions = box_vertices; ctx.indices = duplicated_indices; - BA(CREATE(dev, 2, IDS, IFA, npos, POS, &ctx, &scn)); + BA(sdis_scene_create(dev, &scn_args, &scn)); /* Degenerated triangle */ ctx.indices = degenerated_indices; - BA(CREATE(dev, 1, IDS, IFA, npos, POS, &ctx, &scn)); - + BA(sdis_scene_create(dev, &scn_args, &scn)); ctx.positions = box_vertices; ctx.indices = box_indices; - OK(CREATE(dev, ntris, IDS, IFA, npos, POS, &ctx, &scn)); - - #undef CREATE - #undef IDS - #undef POS - #undef IFA + BA(sdis_scene_create(dev, &SDIS_SCENE_CREATE_ARGS_DEFAULT, &scn)); + BA(sdis_scene_create(NULL, &scn_args, &scn)); + BA(sdis_scene_create(dev, NULL, &scn)); + BA(sdis_scene_create(dev, &scn_args, NULL)); + OK(sdis_scene_create(dev, &scn_args, &scn)); BA(sdis_scene_get_dimension(NULL, &dim)); BA(sdis_scene_get_dimension(scn, NULL)); @@ -242,9 +257,10 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf) size_t degenerated_indices[] = { 0, 0 }; double duplicated_vertices[] = { 0, 0, 0, 0 }; struct sdis_scene* scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; double lower[2], upper[2]; double u0, u1, u2, pos[2], pos1[2]; - double dst; + double dst, fp, t; struct context ctx; struct senc2d_scene* scn2d; struct senc3d_scene* scn3d; @@ -260,37 +276,50 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf) nsegs = square_nsegments; npos = square_nvertices; - #define CREATE sdis_scene_2d_create - #define IDS get_indices_2d - #define POS get_position_2d - #define IFA get_interface - - BA(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL)); - BA(CREATE(dev, 0, IDS, IFA, npos, POS, &ctx, &scn)); - BA(CREATE(dev, nsegs, NULL, IFA, npos, POS, &ctx, &scn)); - BA(CREATE(dev, nsegs, IDS, NULL, npos, POS, &ctx, &scn)); - BA(CREATE(dev, nsegs, IDS, IFA, 0, POS, &ctx, &scn)); - BA(CREATE(dev, nsegs, IDS, IFA, npos, NULL, &ctx, &scn)); + scn_args.get_indices = get_indices_2d; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position_2d; + scn_args.nprimitives = nsegs; + scn_args.nvertices = npos; + scn_args.context = &ctx; + + BA(sdis_scene_2d_create(NULL, &SDIS_SCENE_CREATE_ARGS_DEFAULT, NULL)); + scn_args.nprimitives = 0; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.nprimitives = nsegs; + scn_args.get_indices = NULL; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.get_indices = get_indices_2d; + scn_args.get_interface = NULL; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.get_interface = get_interface; + scn_args.get_position = NULL; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.get_position = get_position_2d; + scn_args.nvertices = 0; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.nvertices = npos; + scn_args.fp_to_meter = 0; + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); + scn_args.fp_to_meter = 1; /* Duplicated vertex */ ctx.positions = duplicated_vertices; ctx.indices = dup_vrtx_indices; - BA(CREATE(dev, 1, IDS, IFA, npos, POS, &ctx, &scn)); + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); /* Duplicated segment */ ctx.positions = square_vertices; ctx.indices = duplicated_indices; - BA(CREATE(dev, 2, IDS, IFA, npos, POS, &ctx, &scn)); + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); /* Degenerated segment */ ctx.indices = degenerated_indices; - BA(CREATE(dev, 1, IDS, IFA, npos, POS, &ctx, &scn)); - + BA(sdis_scene_2d_create(dev, &scn_args, &scn)); ctx.positions = square_vertices; ctx.indices = square_indices; - OK(CREATE(dev, nsegs, IDS, IFA, npos, POS, &ctx, &scn)); - - #undef CREATE - #undef IDS - #undef POS - #undef IFA + BA(sdis_scene_2d_create(dev, &SDIS_SCENE_CREATE_ARGS_DEFAULT, &scn)); + BA(sdis_scene_2d_create(NULL, &scn_args, &scn)); + BA(sdis_scene_2d_create(dev, NULL, &scn)); + BA(sdis_scene_2d_create(dev, &scn_args, NULL)); + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); BA(sdis_scene_get_dimension(NULL, &dim)); BA(sdis_scene_get_dimension(scn, NULL)); @@ -308,6 +337,46 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf) u0 = 0.5; + BA(sdis_scene_get_fp_to_meter(NULL, NULL)); + BA(sdis_scene_get_fp_to_meter(scn, NULL)); + BA(sdis_scene_get_fp_to_meter(NULL, &fp)); + OK(sdis_scene_get_fp_to_meter(scn, &fp)); + CHK(fp == 1); + + fp = 0; + BA(sdis_scene_set_fp_to_meter(NULL, fp)); + BA(sdis_scene_set_fp_to_meter(scn, fp)); + fp = 2; + OK(sdis_scene_set_fp_to_meter(scn, fp)); + OK(sdis_scene_get_fp_to_meter(scn, &fp)); + CHK(fp == 2); + + BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL)); + BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL)); + BA(sdis_scene_get_ambient_radiative_temperature(NULL, &t)); + OK(sdis_scene_get_ambient_radiative_temperature(scn, &t)); + CHK(t == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad); + + t = 100; + BA(sdis_scene_set_ambient_radiative_temperature(NULL, t)); + OK(sdis_scene_set_ambient_radiative_temperature(scn, t)); + OK(sdis_scene_get_ambient_radiative_temperature(scn, &t)); + CHK(t == 100); + + BA(sdis_scene_get_reference_temperature(NULL, NULL)); + BA(sdis_scene_get_reference_temperature(scn, NULL)); + BA(sdis_scene_get_reference_temperature(NULL, &t)); + OK(sdis_scene_get_reference_temperature(scn, &t)); + CHK(t == SDIS_SCENE_CREATE_ARGS_DEFAULT.tref); + + t = -1; + BA(sdis_scene_set_reference_temperature(NULL, t)); + BA(sdis_scene_set_reference_temperature(scn, t)); + t = 100; + OK(sdis_scene_set_reference_temperature(scn, t)); + OK(sdis_scene_get_reference_temperature(scn, &t)); + CHK(t == 100); + BA(sdis_scene_get_boundary_position(NULL, 1, &u0, pos)); BA(sdis_scene_get_boundary_position(scn, 4, &u0, pos)); BA(sdis_scene_get_boundary_position(scn, 1, NULL, pos)); diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c @@ -270,6 +270,7 @@ main(int argc, char** argv) struct sdis_medium* solid = NULL; struct sdis_interface* interf = NULL; struct sdis_scene* scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -280,7 +281,7 @@ main(int argc, char** argv) struct context ctx; double lower[3]; double upper[3]; - double size[3]; + double spread; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -331,16 +332,19 @@ main(int argc, char** argv) /* Create the scene */ ctx.interf = interf; - OK(sdis_scene_create(dev, ctx.msh.nprimitives, get_indices, get_interface, - ctx.msh.nvertices, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ctx.msh.nprimitives; + scn_args.nvertices = ctx.msh.nvertices; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); /*dump_mesh(stdout, ctx.msh.positions, ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/ /* Compute the delta of the solid random walk */ - size[0] = upper[0] - lower[0]; - size[1] = upper[1] - lower[1]; - size[2] = upper[2] - lower[2]; - solid_param->delta = MMIN(MMIN(size[0], size[1]), size[2]) / 20.0; + OK(sdis_scene_get_medium_spread(scn, solid, &spread)); + solid_param->delta = 0.4 / spread; /* (4V/S) / 10 */ interf_param->upper[0] = upper[0]; interf_param->upper[1] = upper[1]; diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -66,7 +66,6 @@ fluid_get_temperature return ((const struct fluid*)sdis_data_cget(data))->temperature; } - static double solid_get_calorific_capacity (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) @@ -109,7 +108,8 @@ solid_get_temperature { (void)data; CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; + if(vtx->time > 0) return UNKNOWN_TEMPERATURE; + return Tf; } /******************************************************************************* @@ -186,6 +186,7 @@ main(int argc, char** argv) struct sdis_estimator* estimator = NULL; struct sdis_estimator* estimator2 = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -217,7 +218,7 @@ main(int argc, char** argv) OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); OK(sdis_data_ref_put(data)); - /* Create the solid_medium */ + /* Create the solid medium */ solid_shader.calorific_capacity = solid_get_calorific_capacity; solid_shader.thermal_conductivity = solid_get_thermal_conductivity; solid_shader.volumic_mass = solid_get_volumic_mass; @@ -278,14 +279,22 @@ main(int argc, char** argv) square_interfaces[3] = interf_H; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_adiabatic)); @@ -346,9 +355,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &probe_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, probe_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -377,9 +386,9 @@ main(int argc, char** argv) OK(GREEN(square_scn, &probe_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, probe_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); @@ -389,6 +398,22 @@ main(int argc, char** argv) fluid_param->temperature = UNKNOWN_TEMPERATURE; BA(SOLVE(square_scn, &probe_args, &estimator)); fluid_param->temperature = Tf; + + /* Right-side temperature at initial time */ + probe_args.time_range[0] = 0; + probe_args.time_range[1] = 0; + + probe_args.iprim = 6; + OK(SOLVE(box_scn, &probe_args, &estimator)); + check_estimator(estimator, N, Tf); + + OK(sdis_estimator_ref_put(estimator)); + + probe_args.iprim = 3; + OK(SOLVE(square_scn, &probe_args, &estimator)); + check_estimator(estimator, N, Tf); + + OK(sdis_estimator_ref_put(estimator)); #undef F #undef SOLVE @@ -469,9 +494,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, bound_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -506,9 +531,9 @@ main(int argc, char** argv) OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, bound_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -539,9 +564,9 @@ main(int argc, char** argv) OK(GREEN(box_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, box_scn, bound_args.time_range); + check_green_serialization(green, box_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); @@ -557,13 +582,33 @@ main(int argc, char** argv) OK(GREEN(square_scn, &bound_args, &green)); check_green_function(green); - OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_estimator(estimator2, N, ref); - check_green_serialization(green, square_scn, bound_args.time_range); + check_green_serialization(green, square_scn); OK(sdis_green_function_ref_put(green)); OK(sdis_estimator_ref_put(estimator)); OK(sdis_estimator_ref_put(estimator2)); + + /* Right-side temperature at initial time */ + bound_args.time_range[0] = 0; + bound_args.time_range[1] = 0; + + prims[0] = 6; + prims[1] = 7; + bound_args.nprimitives = 2; + OK(SOLVE(box_scn, &bound_args, &estimator)); + check_estimator(estimator, N, Tf); + + OK(sdis_estimator_ref_put(estimator)); + + prims[0] = 3; + bound_args.nprimitives = 1; + OK(SOLVE(square_scn, &bound_args, &estimator)); + check_estimator(estimator, N, Tf); + + OK(sdis_estimator_ref_put(estimator)); + #undef SOLVE #undef GREEN diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c @@ -22,7 +22,7 @@ * The scene is composed of a solid cube/square whose temperature is unknown. * The convection coefficient with the surrounding fluid is null excepted for * the X faces whose value is 'H'. The Temperature T of the -X face is fixed - * to Tb. The ambiant radiative temperature is 0 excepted for the X faces + * to Tb. The ambient radiative temperature is 0 excepted for the X faces * whose value is 'Trad'. * This test computes temperature and fluxes on the X faces and check that * they are equal to: @@ -229,6 +229,7 @@ main(int argc, char** argv) struct sdis_scene* box_scn = NULL; struct sdis_scene* square_scn = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -325,14 +326,26 @@ main(int argc, char** argv) square_interfaces[3] = interf_H; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.trad = Trad; + scn_args.tref = Tref; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.trad = Trad; + scn_args.tref = Tref; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_adiabatic)); @@ -351,8 +364,6 @@ main(int argc, char** argv) 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)); @@ -407,8 +418,6 @@ main(int argc, char** argv) 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)); diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -540,6 +540,7 @@ main(int argc, char** argv) struct sdis_interface* interf0 = NULL; struct sdis_interface* interf1 = NULL; struct sdis_scene* scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_solve_camera_args solve_args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT; struct ssp_rng* rng_state = NULL; struct fluid fluid_param = FLUID_NULL; @@ -608,9 +609,15 @@ main(int argc, char** argv) /* Setup the scene */ ntris = sa_size(geom.indices) / 3; /* #primitives */ npos = sa_size(geom.positions) / 3; /* #positions */ - OK(sdis_scene_create(dev, ntris, geometry_get_indices, - geometry_get_interface, npos, geometry_get_position, - &geom, &scn)); + scn_args.get_indices = geometry_get_indices; + scn_args.get_interface = geometry_get_interface; + scn_args.get_position = geometry_get_position; + scn_args.nprimitives = ntris; + scn_args.nvertices = npos; + scn_args.trad = 300; + scn_args.tref = 300; + scn_args.context = &geom; + OK(sdis_scene_create(dev, &scn_args, &scn)); #if 0 dump_mesh(stdout, geom.positions, sa_size(geom.positions)/3, geom.indices, @@ -635,8 +642,6 @@ main(int argc, char** argv) 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)); @@ -645,15 +650,9 @@ main(int argc, char** argv) solve_args.cam = NULL; BA(sdis_solve_camera(scn, &solve_args, &buf)); solve_args.cam = cam; - solve_args.fp_to_meter = 0; + OK(sdis_scene_set_ambient_radiative_temperature(scn, -1)); 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; + OK(sdis_scene_set_ambient_radiative_temperature(scn, 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; diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -220,6 +220,7 @@ main(int argc, char** argv) struct fluid* fluid_param = NULL; struct solid* solid_param = NULL; struct interf* interface_param = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; @@ -337,8 +338,13 @@ main(int argc, char** argv) } #endif - OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ntris; + scn_args.nvertices = nverts; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); BA(sdis_scene_get_medium_spread(NULL, solid0, &v0)); BA(sdis_scene_get_medium_spread(scn, NULL, &v0)); @@ -366,9 +372,6 @@ main(int argc, char** argv) 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; @@ -418,8 +421,7 @@ main(int argc, char** argv) OK(sdis_scene_ref_put(scn)); ctx.interf0 = solid0_fluid0; ctx.interf1 = solid0_fluid1; - OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); + OK(sdis_scene_create(dev, &scn_args, &scn)); OK(sdis_scene_get_medium_spread(scn, solid0, &v)); CHK(eq_eps(v, v0+v1, 1.e-6)); @@ -451,15 +453,12 @@ main(int argc, char** argv) 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)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -205,6 +205,7 @@ main(int argc, char** argv) struct fluid* fluid_param = NULL; struct solid* solid_param = NULL; struct interf* interface_param = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; @@ -325,8 +326,13 @@ main(int argc, char** argv) ctx.nsegments_interf0 = square_nsegments; ctx.interf0 = solid0_fluid0; ctx.interf1 = solid1_fluid1; - OK(sdis_scene_2d_create(dev, sa_size(indices)/2, get_indices, get_interface, - sa_size(positions)/2, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = sa_size(indices)/2; + scn_args.nvertices = sa_size(positions)/2; + scn_args.context = &ctx; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); OK(sdis_scene_get_medium_spread(scn, solid0, &a0)); CHK(eq_eps(a0, 1.0, 1.e-6)); @@ -370,8 +376,7 @@ main(int argc, char** argv) OK(sdis_scene_ref_put(scn)); ctx.interf0 = solid0_fluid0; ctx.interf1 = solid0_fluid1; - OK(sdis_scene_2d_create(dev, sa_size(indices)/2, get_indices, get_interface, - sa_size(positions)/2, get_position, &ctx, &scn)); + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); OK(sdis_scene_get_medium_spread(scn, solid0, &a)); CHK(eq_eps(a, a0+a1, 1.e-6)); @@ -399,10 +404,10 @@ main(int argc, char** argv) 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)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); OK(sdis_green_function_ref_put(green)); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -262,6 +262,7 @@ main(int argc, char** argv) struct sdis_estimator* estimator3 = NULL; struct sdis_green_function* green = NULL; const struct sdis_heat_path* path = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; @@ -335,8 +336,13 @@ main(int argc, char** argv) ctx.positions = box_vertices; ctx.indices = box_indices; ctx.interf = interf; - OK(sdis_scene_create(dev, box_ntriangles, get_indices, get_interface, - box_nvertices, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); OK(sdis_interface_ref_put(interf)); @@ -354,9 +360,6 @@ main(int argc, char** argv) 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; @@ -433,15 +436,12 @@ main(int argc, char** argv) 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, solve_args.time_range, NULL)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + BA(sdis_green_function_solve(NULL, &estimator2)); + BA(sdis_green_function_solve(green, NULL)); + BA(sdis_green_function_solve(NULL, NULL)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); @@ -450,8 +450,9 @@ main(int argc, char** argv) OK(sdis_estimator_ref_put(estimator2)); printf("\n"); - /* Check green used at a different temperature */ + /* Check same green used at a different temperature */ fluid_param->temperature = 500; + OK(sdis_solve_probe(scn, &solve_args, &estimator)); OK(sdis_estimator_get_realisation_count(estimator, &nreals)); OK(sdis_estimator_get_failure_count(estimator, &nfails)); @@ -468,10 +469,10 @@ main(int argc, char** argv) CHK(nfails < N / 1000); CHK(eq_eps(T.E, ref, T.SE)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - + stream = tmpfile(); CHK(stream); BA(sdis_green_function_write(NULL, stream)); @@ -491,7 +492,7 @@ main(int argc, char** argv) OK(sdis_green_function_create_from_stream(scn, stream, &green)); CHK(!fclose(stream)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator3)); + OK(sdis_green_function_solve(green, &estimator3)); check_green_function(green); check_estimator_eq_strict(estimator2, estimator3); @@ -532,6 +533,39 @@ main(int argc, char** argv) OK(sdis_estimator_for_each_path(estimator, process_heat_path, &dump_ctx)); OK(sdis_estimator_ref_put(estimator)); + + printf("\n"); + + /* Green and ambient radiative temperature */ + solve_args.nrealisations = N; + OK(sdis_scene_set_ambient_radiative_temperature(scn, 300)); + OK(sdis_scene_set_reference_temperature(scn, 300)); + + interface_param->epsilon = 1; + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + OK(sdis_green_function_solve(green, &estimator2)); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + + /* Check same green used at different ambient radiative temperature */ + OK(sdis_scene_set_ambient_radiative_temperature(scn, 600)); + + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + OK(sdis_green_function_solve(green, &estimator2)); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); + OK(sdis_scene_ref_put(scn)); OK(sdis_device_ref_put(dev)); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -159,6 +159,7 @@ main(int argc, char** argv) struct sdis_interface* T350 = NULL; struct sdis_scene* scn = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -230,8 +231,13 @@ main(int argc, char** argv) ctx.positions = box_vertices; ctx.indices = box_indices; ctx.interfaces = interfaces; - OK(sdis_scene_create(dev, box_ntriangles, get_indices, get_interface, - box_nvertices, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(Tnone)); @@ -266,10 +272,10 @@ main(int argc, char** argv) /* Check green */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -155,6 +155,7 @@ main(int argc, char** argv) struct sdis_interface* T350 = NULL; struct sdis_scene* scn = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -228,8 +229,13 @@ main(int argc, char** argv) ctx.positions = square_vertices; ctx.indices = square_indices; ctx.interfaces = interfaces; - OK(sdis_scene_2d_create(dev, square_nsegments, get_indices, get_interface, - square_nvertices, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = &ctx; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(Tnone)); @@ -264,10 +270,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -181,6 +181,7 @@ main(int argc, char** argv) struct sdis_interface* solid_solid = NULL; struct sdis_scene* scn = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -282,8 +283,13 @@ main(int argc, char** argv) ctx.solid_solid = solid_solid; nverts = sa_size(ctx.positions) / 3; ntris = sa_size(ctx.indices) / 3; - OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ntris; + scn_args.nvertices = nverts; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); /* Release the scene data */ OK(sdis_interface_ref_put(Tnone)); @@ -321,10 +327,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -178,6 +178,7 @@ main(int argc, char** argv) struct sdis_interface* solid_solid = NULL; struct sdis_scene* scn = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -276,8 +277,13 @@ main(int argc, char** argv) ctx.solid_solid = solid_solid; nverts = sa_size(ctx.positions) / 2; nsegs = sa_size(ctx.indices) / 2; - OK(sdis_scene_2d_create(dev, nsegs, get_indices, get_interface, nverts, - get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = nsegs; + scn_args.nvertices = nverts; + scn_args.context = &ctx; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); /* Release the scene data */ OK(sdis_interface_ref_put(Tnone)); @@ -313,10 +319,10 @@ main(int argc, char** argv) /* Check green function */ OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, solve_args.time_range); + check_green_serialization(green, scn); /* Release data */ OK(sdis_estimator_ref_put(estimator)); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -145,6 +145,7 @@ main(int argc, char** argv) struct sdis_estimator* estimator = NULL; struct sdis_estimator* estimator2 = NULL; struct sdis_green_function* green = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; @@ -192,8 +193,13 @@ main(int argc, char** argv) ctx.positions = square_vertices; ctx.indices = square_indices; ctx.interf = interf; - OK(sdis_scene_2d_create(dev, square_nsegments, get_indices, get_interface, - square_nvertices, get_position, &ctx, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = &ctx; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); OK(sdis_interface_ref_put(interf)); @@ -221,7 +227,7 @@ main(int argc, char** argv) CHK(eq_eps(T.E, ref, T.SE)); OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); - OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2)); + OK(sdis_green_function_solve(green, &estimator2)); check_green_function(green); check_estimator_eq(estimator, estimator2); diff --git a/src/test_sdis_transcient.c b/src/test_sdis_transcient.c @@ -144,7 +144,7 @@ matriochka_indices(const size_t itri, size_t ids[3], void* context) } static void -matriocka_interface +matriochka_interface (const size_t itri, struct sdis_interface** bound, void* context) { struct matriochka_context* ctx = context; @@ -471,6 +471,7 @@ main(int argc, char** argv) struct sdis_medium* fluid = NULL; struct sdis_medium* solid = NULL; struct sdis_data* data = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -557,8 +558,13 @@ main(int argc, char** argv) ctx.scale = boxsz; /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, get_indices, get_interface, - box_nvertices, get_position, &ctx, &box_scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Setup the box2 scene context */ ctx.indices = indices; @@ -576,9 +582,10 @@ main(int argc, char** argv) ctx.interfs[20] = ctx.interfs[21] = interfs[5]; /* Zmax */ ctx.scale = boxsz; - /* Create the box scene */ - OK(sdis_scene_create(dev, ntriangles, get_indices, get_interface, - nvertices, get_position, &ctx, &box2_scn)); + /* Create the box2 scene */ + scn_args.nprimitives = ntriangles; + scn_args.nvertices = nvertices; + OK(sdis_scene_create(dev, &scn_args, &box2_scn)); /* Setup the matriochka context */ matriochka_ctx.interfs[0] = matriochka_ctx.interfs[1] = interfs[4]; /* Zmin */ @@ -589,12 +596,16 @@ main(int argc, char** argv) matriochka_ctx.interfs[10] = matriochka_ctx.interfs[11] = interfs[2]; /* Ymin */ matriochka_ctx.interfs[12] = interfs[6]; /* The remaining internal triangles */ matriochka_ctx.scale = boxsz; - matriochka_ctx.nboxes = nmatriochkas; + matriochka_ctx.nboxes = nmatriochkas; /* Create the matriochka scene */ - OK(sdis_scene_create(dev, box_ntriangles*nmatriochkas, matriochka_indices, - matriocka_interface, box_nvertices*nmatriochkas, matriochka_position, - &matriochka_ctx, &box_matriochka_scn)); + scn_args.get_indices = matriochka_indices; + scn_args.get_interface = matriochka_interface; + scn_args.get_position = matriochka_position; + scn_args.nprimitives = box_ntriangles*nmatriochkas; + scn_args.nvertices = box_nvertices*nmatriochkas; + scn_args.context = &matriochka_ctx; + OK(sdis_scene_create(dev, &scn_args, &box_matriochka_scn)); /* Setup and run the simulation */ probe[0] = 0.1; diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c @@ -94,9 +94,10 @@ solve_green_path(struct sdis_green_path* path, void* ctx) struct green_accum* acc = NULL; struct sdis_data* data = NULL; enum sdis_medium_type type; + enum sdis_green_path_end_type end_type; double power = 0; double flux = 0; - double temp = 0; + double time, temp = 0; double weight = 0; CHK(path && ctx); @@ -110,14 +111,42 @@ solve_green_path(struct sdis_green_path* path, void* ctx) BA(sdis_green_path_for_each_flux_term(path, NULL, &acc)); OK(sdis_green_path_for_each_flux_term(path, accum_flux_terms, &flux)); + BA(sdis_green_path_get_elapsed_time(NULL, NULL)); + BA(sdis_green_path_get_elapsed_time(path, NULL)); + BA(sdis_green_path_get_elapsed_time(NULL, &time)); + OK(sdis_green_path_get_elapsed_time(path, &time)); + + BA(sdis_green_path_get_end_type(NULL, NULL)); + BA(sdis_green_path_get_end_type(path, NULL)); + BA(sdis_green_path_get_end_type(NULL, &end_type)); + OK(sdis_green_path_get_end_type(path, &end_type)); + + BA(sdis_green_path_get_limit_point(NULL, NULL)); BA(sdis_green_path_get_limit_point(NULL, &pt)); BA(sdis_green_path_get_limit_point(path, NULL)); - OK(sdis_green_path_get_limit_point(path, &pt)); - - switch(pt.type) { + if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) { + struct sdis_green_function* green; + struct sdis_scene* scn; + BO(sdis_green_path_get_limit_point(path, &pt)); + BA(sdis_green_path_get_green_function(NULL, NULL)); + BA(sdis_green_path_get_green_function(path, NULL)); + BA(sdis_green_path_get_green_function(NULL, &green)); + OK(sdis_green_path_get_green_function(path, &green)); + + BA(sdis_green_function_get_scene(NULL, NULL)); + BA(sdis_green_function_get_scene(NULL, &scn)); + BA(sdis_green_function_get_scene(green, NULL)); + OK(sdis_green_function_get_scene(green, &scn)); + + BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL)); + BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL)); + BA(sdis_scene_get_ambient_radiative_temperature(NULL, &temp)); + OK(sdis_scene_get_ambient_radiative_temperature(scn, &temp)); + } else { + OK(sdis_green_path_get_limit_point(path, &pt)); + switch(pt.type) { case SDIS_FRAGMENT: frag = pt.data.itfrag.fragment; - frag.time = INF; OK(sdis_interface_get_shader(pt.data.itfrag.intface, &interf)); data = sdis_interface_get_data(pt.data.itfrag.intface); temp = frag.side == SDIS_FRONT @@ -126,7 +155,6 @@ solve_green_path(struct sdis_green_path* path, void* ctx) break; case SDIS_VERTEX: vtx = pt.data.mdmvert.vertex; - vtx.time = INF; type = sdis_medium_get_type(pt.data.mdmvert.medium); data = sdis_medium_get_data(pt.data.mdmvert.medium); if(type == SDIS_FLUID) { @@ -138,6 +166,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx) } break; default: FATAL("Unreachable code.\n"); break; + } } weight = temp + power + flux; @@ -229,7 +258,7 @@ check_green_function(struct sdis_green_function* green) time_range[0] = time_range[1] = INF; - OK(sdis_green_function_solve(green, time_range, &estimator)); + OK(sdis_green_function_solve(green, &estimator)); BA(sdis_green_function_get_paths_count(NULL, &n)); BA(sdis_green_function_get_paths_count(green, NULL)); @@ -367,15 +396,14 @@ dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator) void check_green_serialization (struct sdis_green_function* green, - struct sdis_scene* scn, - const double time_range[2]) + struct sdis_scene* scn) { FILE* stream = NULL; struct sdis_estimator *e1 = NULL; struct sdis_estimator *e2 = NULL; struct sdis_green_function* green2 = NULL; - CHK(green && time_range); + CHK(green && scn); stream = tmpfile(); CHK(stream); @@ -386,8 +414,8 @@ check_green_serialization CHK(!fclose(stream)); check_green_function(green2); - OK(sdis_green_function_solve(green, time_range, &e1)); - OK(sdis_green_function_solve(green2, time_range, &e2)); + OK(sdis_green_function_solve(green, &e1)); + OK(sdis_green_function_solve(green2, &e2)); check_estimator_eq_strict(e1, e2); OK(sdis_estimator_ref_put(e1)); diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -26,6 +26,7 @@ #define OK(Cond) CHK((Cond) == RES_OK) #define BA(Cond) CHK((Cond) == RES_BAD_ARG) +#define BO(Cond) CHK((Cond) == RES_BAD_OP) /******************************************************************************* * Box geometry @@ -354,8 +355,7 @@ dump_heat_paths extern LOCAL_SYM void check_green_serialization (struct sdis_green_function* green, - struct sdis_scene* scn, - const double time_range[2]); + struct sdis_scene* scn); #endif /* TEST_SDIS_UTILS_H */ diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -19,13 +19,14 @@ #include <rsys/clock_time.h> #include <rsys/double3.h> #include <rsys/math.h> +#include <star/ssp.h> /* * The scene is composed of a solid cube/square whose temperature is unknown. * The convection coefficient with the surrounding fluid is null everywhere The * Temperature of the +/- X faces are fixed to T0, and the solid has a volumic * power of P0. This test computes the temperature of a probe position pos into - * the solid and check that it is is equal to: + * the solid and check that at t=inf it is is equal to: * * T(pos) = P0 / (2*LAMBDA) * (A^2/4 - (pos-0.5)^2) + T0 * @@ -61,6 +62,8 @@ struct solid { double cp; double delta; double vpower; + double initial_temperature; + double t0; }; static double @@ -108,9 +111,15 @@ static double solid_get_temperature (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) { - (void)data; + double t0; CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; + CHK(data != NULL); + t0 = ((const struct solid*)sdis_data_cget(data))->t0; + if(vtx->time > t0) { + return UNKNOWN_TEMPERATURE; + } else { + return ((const struct solid*)sdis_data_cget(data))->initial_temperature; + } } static double @@ -152,7 +161,7 @@ interface_get_convection_coef static void solve (struct sdis_scene* scn, - const double pos[], + struct ssp_rng* rng, struct solid* solid) { char dump[128]; @@ -165,142 +174,203 @@ solve struct sdis_mc time = SDIS_MC_NULL; size_t nreals; size_t nfails; - double x; - double ref; - const double time_range[2] = {INF, INF}; + double ref = -1; enum sdis_scene_dimension dim; - ASSERT(scn && pos && solid); - - /* Restore power value */ - solid->vpower = P0; - - x = pos[0] - 0.5; - ref = solid->vpower / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; + const int nsimuls = 4; + int isimul; + ASSERT(scn && rng && solid); 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, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - 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)); - - switch(dim) { + FOR_EACH(isimul, 0, nsimuls) { + int steady = (isimul % 2) == 0; + + /* Restore power value */ + solid->vpower = P0; + + solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9); + solve_args.position[2] = + dim == SDIS_SCENE_2D ? 0 : ssp_rng_uniform_double(rng, 0.1, 0.9); + + solve_args.nrealisations = N; + if(steady) { + solve_args.time_range[0] = solve_args.time_range[1] = INF; + } else { + solve_args.time_range[0] = 100 * (double)isimul; + solve_args.time_range[1] = 4 * solve_args.time_range[0]; + } + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + 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)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g)th Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); - - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, T.SE*3)); - - /* Check green function */ - time_current(&t0); - 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); - - OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); - OK(sdis_estimator_get_failure_count(estimator2, &nfails)); - OK(sdis_estimator_get_temperature(estimator2, &T)); - - switch(dim) { + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); + + /* Check green function */ + time_current(&t0); + OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); + time_current(&t1); + OK(sdis_green_function_solve(green, &estimator2)); + time_current(&t2); + + OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); + OK(sdis_estimator_get_failure_count(estimator2, &nfails)); + OK(sdis_estimator_get_temperature(estimator2, &T)); + + switch(dim) { case SDIS_SCENE_2D: - printf("Green temperature at (%g %g) with Power=%g = %g ~ %g +/- %g\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Green Steady temperature at (%g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Green Mean temperature at (%g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; case SDIS_SCENE_3D: - printf("Green temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE); + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Green Steady temperature at (%g, %g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Green Mean temperature at (%g, %g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } break; default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - time_sub(&t0, &t1, &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green estimation time = %s\n", dump); - time_sub(&t1, &t2, &t1); - time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); - - check_green_function(green); - check_estimator_eq(estimator, estimator2); - check_green_serialization(green, scn, time_range); - - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - printf("\n"); - - /* Check green used at a different power */ - solid->vpower = 3 * P0; - - time_current(&t0); - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_time(estimator, &time)); - - ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; - - switch (dim) { - case SDIS_SCENE_2D: - printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); - break; - case SDIS_SCENE_3D: - printf("Temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE); - break; - default: FATAL("Unreachable code.\n"); break; - } - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - printf("Elapsed time = %s\n", dump); - printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + time_sub(&t0, &t1, &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green estimation time = %s\n", dump); + time_sub(&t1, &t2, &t1); + time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); + + check_green_function(green); + check_estimator_eq(estimator, estimator2); + check_green_serialization(green, scn); + + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + printf("\n"); + + /* Check same green used at a different power level */ + solid->vpower = 3 * P0; + + time_current(&t0); + OK(sdis_solve_probe(scn, &solve_args, &estimator)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + + 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)); + + switch(dim) { + case SDIS_SCENE_2D: + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g" + " ~ %g +/- %g\n", + SPLIT2(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } + break; + case SDIS_SCENE_3D: + if(steady) { + double x = solve_args.position[0] - 0.5; + ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; + printf("Steady temperature at (%g, %g, %g) with Power=%g = %g" + " ~ %g +/- %g\n", + SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); + } else { + printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and" + " Power=%g ~ %g +/- %g\n", + SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), + solid->vpower, T.E, T.SE); + } + break; + default: FATAL("Unreachable code.\n"); break; + } + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); - CHK(nfails + nreals == N); - CHK(nfails < N / 1000); - CHK(eq_eps(T.E, ref, T.SE * 3)); + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + if(steady) CHK(eq_eps(T.E, ref, T.SE * 3)); - time_current(&t0); - OK(sdis_green_function_solve(green, time_range, &estimator2)); - time_sub(&t0, time_current(&t1), &t0); - time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); - printf("Green solve time = %s\n", dump); + time_current(&t0); + OK(sdis_green_function_solve(green, &estimator2)); + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); + printf("Green solve time = %s\n", dump); - check_green_function(green); - check_estimator_eq(estimator, estimator2); + check_green_function(green); + check_estimator_eq(estimator, estimator2); - OK(sdis_estimator_ref_put(estimator)); - OK(sdis_estimator_ref_put(estimator2)); - OK(sdis_green_function_ref_put(green)); + OK(sdis_estimator_ref_put(estimator)); + OK(sdis_estimator_ref_put(estimator2)); + OK(sdis_green_function_ref_put(green)); - printf("\n"); + printf("\n\n"); + } } /******************************************************************************* @@ -318,6 +388,7 @@ main(int argc, char** argv) struct sdis_interface* interf_T0 = NULL; struct sdis_scene* box_scn = NULL; struct sdis_scene* square_scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -325,7 +396,7 @@ main(int argc, char** argv) struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; struct solid* solid_props = NULL; - double pos[3]; + struct ssp_rng* rng = NULL; (void)argc, (void)argv; OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); @@ -350,6 +421,8 @@ main(int argc, char** argv) solid_props->rho = 25; solid_props->delta = DELTA; solid_props->vpower = P0; + solid_props->t0 = 0; + solid_props->initial_temperature = T0; OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); @@ -392,29 +465,38 @@ main(int argc, char** argv) square_interfaces[3] = interf_T0; /* Right */ /* Create the box scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = box_get_interface; + scn_args.get_position = box_get_position; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = box_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); /* Create the square scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = square_get_interface; + scn_args.get_position = square_get_position; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = square_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_adiabatic)); OK(sdis_interface_ref_put(interf_T0)); /* Solve */ - d3_splat(pos, 0.25); + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); printf(">> Box scene\n"); - solve(box_scn, pos, solid_props); + solve(box_scn, rng, solid_props); printf(">> Square scene\n"); - solve(square_scn, pos, solid_props); + solve(square_scn, rng, solid_props); OK(sdis_scene_ref_put(box_scn)); OK(sdis_scene_ref_put(square_scn)); OK(sdis_device_ref_put(dev)); + OK(ssp_rng_ref_put(rng)); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -274,6 +274,7 @@ main(int argc, char** argv) struct sdis_medium* solid1 = NULL; struct sdis_medium* solid2 = NULL; struct sdis_scene* scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -425,8 +426,13 @@ main(int argc, char** argv) interfaces[17] = interf_solid2_adiabatic; /* Create the scene */ - OK(sdis_scene_create(dev, ntriangles, get_indices, get_interface, - nvertices, get_position, interfaces, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ntriangles; + scn_args.nvertices = nvertices; + scn_args.context = interfaces; + OK(sdis_scene_create(dev, &scn_args, &scn)); #if 0 dump_mesh(stdout, vertices, nvertices, indices, ntriangles); @@ -441,8 +447,7 @@ main(int argc, char** argv) data = sdis_medium_get_data(solid1); solid_param = sdis_data_get(data); solid_param->lambda = 0.1; - OK(sdis_scene_create(dev, ntriangles, get_indices, get_interface, - nvertices, get_position, interfaces, &scn)); + OK(sdis_scene_create(dev, &scn_args, &scn)); printf("\n>>> Check 2\n"); check(scn, refs2, sizeof(refs2)/sizeof(struct reference)); diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -296,6 +296,7 @@ main(int argc, char** argv) struct sdis_medium* solid1 = NULL; struct sdis_medium* solid2 = NULL; struct sdis_scene* scn = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -446,8 +447,13 @@ main(int argc, char** argv) interfaces[7] = interf_solid1_solid2; /* Create the scene */ - OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, - nvertices, get_position, interfaces, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = nsegments; + scn_args.nvertices = nvertices; + scn_args.context = interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); printf(">>> Check 1\n"); check(scn, refs1, sizeof(refs1)/sizeof(struct reference)); @@ -457,8 +463,7 @@ main(int argc, char** argv) data = sdis_medium_get_data(solid1); solid_param = sdis_data_get(data); solid_param->lambda = 0.1; - OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, - nvertices, get_position, interfaces, &scn) ); + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); printf("\n>>> Check 2\n"); check(scn, refs2, sizeof(refs2)/sizeof(struct reference)); @@ -472,8 +477,7 @@ main(int argc, char** argv) solid_param = sdis_data_get(data); solid_param->lambda = 10; solid_param->P = SDIS_VOLUMIC_POWER_NONE; - OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, - nvertices, get_position, interfaces, &scn)); + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); printf("\n>>> Check 3\n"); check(scn, refs3, sizeof(refs3)/sizeof(struct reference)); 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_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; struct sdis_mc T = SDIS_MC_NULL; double Tref; @@ -414,8 +415,13 @@ main(int argc, char** argv) #endif /* Create the scene */ - OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, - nvertices, get_position, interfaces, &scn)); + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = nsegments; + scn_args.nvertices = nvertices; + scn_args.context = interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &scn)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_solid_adiabatic)); diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c @@ -217,6 +217,7 @@ main(int argc, char** argv) struct sdis_scene* scn_2d = NULL; struct sdis_scene* scn_3d = NULL; struct sdis_estimator* estimator = NULL; + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; @@ -326,8 +327,13 @@ main(int argc, char** argv) interfaces[3] = interf_adiabatic; /* Right */ /* Create the 2D scene */ - OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, get_interface, - square_nvertices, get_position_2d, interfaces, &scn_2d)); + scn_args.get_indices = square_get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position_2d; + scn_args.nprimitives = square_nsegments; + scn_args.nvertices = square_nvertices; + scn_args.context = interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &scn_2d)); /* Map the interfaces to their box triangles */ interfaces[0] = interfaces[1] = interf_adiabatic; /* Front */ @@ -338,8 +344,13 @@ main(int argc, char** argv) interfaces[10]= interfaces[11]= interf_solid_fluid2; /* Bottom */ /* Create the 3D scene */ - OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, get_interface, - box_nvertices, get_position_3d, interfaces, &scn_3d)); + scn_args.get_indices = box_get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position_3d; + scn_args.nprimitives = box_ntriangles; + scn_args.nvertices = box_nvertices; + scn_args.context = interfaces; + OK(sdis_scene_create(dev, &scn_args, &scn_3d)); /* Release the interfaces */ OK(sdis_interface_ref_put(interf_adiabatic));