stardis-solver

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

commit 71401c573c36f3519de415fa788d94342d6f8dcb
parent 01ffc837ac8f1ae211ab49b3c9b6f75ea7176cdc
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  9 Oct 2018 12:18:51 +0200

Test the sdis_solve_boundary function

Diffstat:
Mcmake/CMakeLists.txt | 2+-
Msrc/sdis_solve_Xd.h | 21+++++++++++++++------
Asrc/test_sdis_solve_boundary.c | 373+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dsrc/test_sdis_solve_probe_boundary.c | 305------------------------------------------------------------------------------
4 files changed, 389 insertions(+), 312 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -155,7 +155,7 @@ if(NOT NO_TEST) new_test(test_sdis_solve_probe_2d) new_test(test_sdis_solve_probe2_2d) new_test(test_sdis_solve_probe3_2d) - new_test(test_sdis_solve_probe_boundary) + new_test(test_sdis_solve_boundary) new_test(test_sdis_volumic_power) # Additionnal tests diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -1592,10 +1592,10 @@ XD(solve_boundary) struct ssp_rng** rngs = NULL; size_t i; size_t N = 0; /* #realisations that do not fail */ + size_t view_nprims; double weight=0, sqr_weight=0; int64_t irealisation; ATOMIC res = RES_OK; - ASSERT(scene_is_2d(scn)); if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives || !sides || !nprimitives || time < 0 || fp_to_meter < 0 || Tref < 0 @@ -1604,6 +1604,14 @@ XD(solve_boundary) goto error; } + SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims)); + FOR_EACH(i, 0, nprimitives) { + if(primitives[i] >= view_nprims) { + res = RES_BAD_ARG; + goto error; + } + } + /* Create the Star-XD shape of the boundary */ #if DIM == 2 res = sXd(shape_create_line_segments)(scn->dev->sXd_dev, &shape); @@ -1668,10 +1676,6 @@ XD(solve_boundary) if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ - /* Map from boundary scene to sdis scene */ - iprim = primitives[prim.prim_id]; - side = sides[prim.prim_id]; - /* Sample a position onto the boundary */ #if DIM == 2 res_local = s2d_scene_view_sample @@ -1679,7 +1683,7 @@ XD(solve_boundary) ssp_rng_canonical_float(rng), ssp_rng_canonical_float(rng), &prim, st); - uv[1] = (double)st[1]; + uv[0] = (double)st[0]; #else res_local = s3d_scene_view_sample (view, @@ -1691,6 +1695,11 @@ XD(solve_boundary) #endif if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } + /* Map from boundary scene to sdis scene */ + ASSERT(prim.prim_id < nprimitives); + iprim = primitives[prim.prim_id]; + side = sides[prim.prim_id]; + /* Invoke the boundary realisation */ res_local = XD(boundary_realisation) (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w); diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -0,0 +1,373 @@ +/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sdis.h" +#include "test_sdis_utils.h" + +#include <rsys/math.h> + +/* + * The scene is composed of a solid cube/square whose temperature is unknown. + * The convection coefficient with the surrounding fluid is null exepted for + * the +X face whose value is 'H'. The Temperature of the -X face is fixed to + * Tb. This test computes the temperature on the +X face and check that it is + * equal to: + * + * T = (H*Tf + LAMBDA/A * Tb) / (H+LAMBDA/A) + * + * with Tf the temperature of the surrounding fluid, lambda the conductivity of + * the cube and A the size of the cube/square, i.e. 1. + * + * 3D 2D + * + * ///// (1,1,1) ///// (1,1) + * +-------+ +-------+ + * /' /| _\ | | _\ + * +-------+ | / / Tf Tb | / / Tf + * Tb +.....|.+ \__/ | | \__/ + * |, |/ +-------+ + * +-------+ (0,0) ///// + * (0,0,0) ///// + */ + +#define UNKNOWN_TEMPERATURE -1 +#define N 10000 /* #realisations */ + +#define Tf 310.0 +#define Tb 300.0 +#define H 0.5 +#define LAMBDA 0.1 + +/******************************************************************************* + * Media + ******************************************************************************/ +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return Tf; +} + +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return 2.0; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return LAMBDA; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return 25.0; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return 1.0/20.0; +} + +static double +solid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return UNKNOWN_TEMPERATURE; +} + +/******************************************************************************* + * Interfaces + ******************************************************************************/ +struct interf { + double temperature; + double hc; +}; + +static double +interface_get_temperature + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf = sdis_data_cget(data); + CHK(frag && data); + return interf->temperature; +} + +static double +interface_get_convection_coef + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf = sdis_data_cget(data); + CHK(frag && data); + return interf->hc; +} + +/******************************************************************************* + * Helper function + ******************************************************************************/ +static void +check_estimator + (const struct sdis_estimator* estimator, + const size_t nrealisations, /* #realisations */ + const double ref) +{ + struct sdis_mc T = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + CHK(estimator && nrealisations); + + CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK); + CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); + CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); + printf("%g ~ %g +/- %g\n", ref, T.E, T.SE); + printf("#failures = %lu/%lu\n", + (unsigned long)nfails, (unsigned long)nrealisations); + CHK(nfails + nreals == nrealisations); + CHK(nfails < N/1000); + CHK(eq_eps(T.E, ref, 3*T.SE)); +} + +/******************************************************************************* + * Test + ******************************************************************************/ +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct sdis_data* data = NULL; + struct sdis_device* dev = NULL; + struct sdis_medium* fluid = NULL; + struct sdis_medium* solid = NULL; + struct sdis_interface* interf_adiabatic = NULL; + struct sdis_interface* interf_Tb = NULL; + struct sdis_interface* interf_H = NULL; + struct sdis_scene* box_scn = NULL; + struct sdis_scene* square_scn = NULL; + struct sdis_estimator* estimator = 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_interface* box_interfaces[12 /*#triangles*/]; + struct sdis_interface* square_interfaces[4/*#segments*/]; + struct interf* interf_props = NULL; + double uv[2]; + double pos[3]; + double ref; + size_t prims[4]; + enum sdis_side sides[4]; + size_t iprim; + (void)argc, (void)argv; + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + CHK(sdis_device_create + (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK); + + /* Create the fluid medium */ + fluid_shader.temperature = fluid_get_temperature; + CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); + + /* 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; + solid_shader.delta_solid = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); + + /* Setup the interface shader */ + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.front.temperature = interface_get_temperature; + interf_shader.front.emissivity = NULL; + interf_shader.front.specular_fraction = NULL; + interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; + + /* Create the adiabatic interface */ + CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); + interf_props = sdis_data_get(data); + interf_props->hc = 0; + interf_props->temperature = UNKNOWN_TEMPERATURE; + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Create the Tb interface */ + CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); + interf_props = sdis_data_get(data); + interf_props->hc = 0; + interf_props->temperature = Tb; + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_Tb) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Create the H interface */ + CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); + interf_props = sdis_data_get(data); + interf_props->hc = H; + interf_props->temperature = UNKNOWN_TEMPERATURE; + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_H) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Release the media */ + CHK(sdis_medium_ref_put(solid) == RES_OK); + CHK(sdis_medium_ref_put(fluid) == RES_OK); + + /* Map the interfaces to their box triangles */ + box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */ + box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */ + box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ + box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */ + box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */ + box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */ + + /* Map the interfaces to their square segments */ + square_interfaces[0] = interf_adiabatic; /* Bottom */ + square_interfaces[1] = interf_Tb; /* Lef */ + square_interfaces[2] = interf_adiabatic; /* Top */ + square_interfaces[3] = interf_H; /* Right */ + + /* Create the box scene */ + CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices, + box_get_interface, box_nvertices, box_get_position, box_interfaces, + &box_scn) == RES_OK); + + /* Create the square scene */ + CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, + square_get_interface, square_nvertices, square_get_position, + square_interfaces, &square_scn) == RES_OK); + + /* Release the interfaces */ + CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK); + CHK(sdis_interface_ref_put(interf_Tb) == RES_OK); + CHK(sdis_interface_ref_put(interf_H) == RES_OK); + + ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); + + #define SOLVE sdis_solve_probe_boundary + #define F SDIS_FRONT + uv[0] = 0.3; + uv[1] = 0.3; + iprim = 6; + + CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG); + + CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); + CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK); + printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos)); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + + uv[0] = 0.5; + iprim = 3; + CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); + CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK); + printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos)); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + #undef F + #undef SOLVE + + sides[0] = SDIS_FRONT; + sides[1] = SDIS_FRONT; + sides[2] = SDIS_FRONT; + sides[3] = SDIS_FRONT; + + #define SOLVE sdis_solve_boundary + prims[0] = 6; + prims[1] = 7; + CHK(SOLVE(NULL, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, 0, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, NULL, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, prims, NULL, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, prims, sides, 0, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, prims, sides, 2, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, NULL) == RES_BAD_ARG); + + /* Average temperature on the right side of the box */ + CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_OK); + printf("Average temperature of the right side of the box = "); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + + /* Average temperature on the right side of the square */ + prims[0] = 3; + sides[0] = SDIS_FRONT; + CHK(SOLVE(square_scn, N, prims, sides, 1, INF, 1.0, 0, 0, &estimator) == RES_OK); + printf("Average temperature of the right side of the square = "); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + + /* Check out of bound prims */ + prims[0] = 12; + CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + prims[0] = 4; + CHK(SOLVE(square_scn, N, prims, sides, 1, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + + ref = (ref + Tb) / 2; + + /* Average temperature on the left/right side of the box */ + prims[0] = 2; + prims[1] = 3; + prims[2] = 6; + prims[3] = 7; + CHK(SOLVE(box_scn, N, prims, sides, 4, INF, 1.0, 0, 0, &estimator) == RES_OK); + printf("Average temperature of the right/left side of the box = "); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + + /* Average temperature on the left/right side of the square */ + prims[0] = 1; + prims[1] = 3; + CHK(SOLVE(square_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_OK); + printf("Average temperature of the right/left side of the square = "); + check_estimator(estimator, N, ref); + CHK(sdis_estimator_ref_put(estimator) == RES_OK); + #undef sdis_solve_boundary + + CHK(sdis_scene_ref_put(box_scn) == RES_OK); + CHK(sdis_scene_ref_put(square_scn) == RES_OK); + CHK(sdis_device_ref_put(dev) == RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c @@ -1,305 +0,0 @@ -/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. */ - -#include "sdis.h" -#include "test_sdis_utils.h" - -#include <rsys/math.h> - -/* - * The scene is composed of a solid cube/square whose temperature is unknown. - * The convection coefficient with the surrounding fluid is null exepted for - * the +X face whose value is 'H'. The Temperature of the -X face is fixed to - * Tb. This test computes the temperature on the +X face and check that it is - * equal to: - * - * T = (H*Tf + LAMBDA/A * Tb) / (H+LAMBDA/A) - * - * with Tf the temperature of the surrounding fluid, lambda the conductivity of - * the cube and A the size of the cube/square, i.e. 1. - * - * 3D 2D - * - * ///// (1,1,1) ///// (1,1) - * +-------+ +-------+ - * /' /| _\ | | _\ - * +-------+ | / / Tf Tb | / / Tf - * Tb +.....|.+ \__/ | | \__/ - * |, |/ +-------+ - * +-------+ (0,0) ///// - * (0,0,0) ///// - */ - -#define UNKNOWN_TEMPERATURE -1 -#define N 10000 /* #realisations */ - -#define Tf 310 -#define Tb 300 -#define H 0.5 -#define LAMBDA 0.1 - -/******************************************************************************* - * Media - ******************************************************************************/ -static double -fluid_get_temperature - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return Tf; -} - -static double -solid_get_calorific_capacity - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return 2.0; -} - -static double -solid_get_thermal_conductivity - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return LAMBDA; -} - -static double -solid_get_volumic_mass - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return 25.0; -} - -static double -solid_get_delta - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return 1.0/20.0; -} - -static double -solid_get_temperature - (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) -{ - (void)data; - CHK(vtx != NULL); - return UNKNOWN_TEMPERATURE; -} - -/******************************************************************************* - * Interfaces - ******************************************************************************/ -struct interf { - double temperature; - double hc; -}; - -static double -interface_get_temperature - (const struct sdis_interface_fragment* frag, struct sdis_data* data) -{ - const struct interf* interf = sdis_data_cget(data); - CHK(frag && data); - return interf->temperature; -} - -static double -interface_get_convection_coef - (const struct sdis_interface_fragment* frag, struct sdis_data* data) -{ - const struct interf* interf = sdis_data_cget(data); - CHK(frag && data); - return interf->hc; -} - -/******************************************************************************* - * Test - ******************************************************************************/ -int -main(int argc, char** argv) -{ - struct mem_allocator allocator; - struct sdis_mc T = SDIS_MC_NULL; - struct sdis_data* data = NULL; - struct sdis_device* dev = NULL; - struct sdis_medium* fluid = NULL; - struct sdis_medium* solid = NULL; - struct sdis_interface* interf_adiabatic = NULL; - struct sdis_interface* interf_Tb = NULL; - struct sdis_interface* interf_H = NULL; - struct sdis_scene* box_scn = NULL; - struct sdis_scene* square_scn = NULL; - struct sdis_estimator* estimator = 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_interface* box_interfaces[12 /*#triangles*/]; - struct sdis_interface* square_interfaces[4/*#segments*/]; - struct interf* interf_props = NULL; - double uv[2]; - double pos[3]; - double ref; - size_t iprim; - size_t nreals; - size_t nfails; - (void)argc, (void)argv; - - CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); - CHK(sdis_device_create - (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK); - - /* Create the fluid medium */ - fluid_shader.temperature = fluid_get_temperature; - CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); - - /* 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; - solid_shader.delta_solid = solid_get_delta; - solid_shader.temperature = solid_get_temperature; - CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); - - /* Setup the interface shader */ - interf_shader.convection_coef = interface_get_convection_coef; - interf_shader.front.temperature = interface_get_temperature; - interf_shader.front.emissivity = NULL; - interf_shader.front.specular_fraction = NULL; - interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; - - /* Create the adiabatic interface */ - CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); - interf_props = sdis_data_get(data); - interf_props->hc = 0; - interf_props->temperature = UNKNOWN_TEMPERATURE; - CHK(sdis_interface_create - (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK); - CHK(sdis_data_ref_put(data) == RES_OK); - - /* Create the Tb interface */ - CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); - interf_props = sdis_data_get(data); - interf_props->hc = 0; - interf_props->temperature = Tb; - CHK(sdis_interface_create - (dev, solid, fluid, &interf_shader, data, &interf_Tb) == RES_OK); - CHK(sdis_data_ref_put(data) == RES_OK); - - /* Create the H interface */ - CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); - interf_props = sdis_data_get(data); - interf_props->hc = H; - interf_props->temperature = UNKNOWN_TEMPERATURE; - CHK(sdis_interface_create - (dev, solid, fluid, &interf_shader, data, &interf_H) == RES_OK); - CHK(sdis_data_ref_put(data) == RES_OK); - - /* Release the media */ - CHK(sdis_medium_ref_put(solid) == RES_OK); - CHK(sdis_medium_ref_put(fluid) == RES_OK); - - /* Map the interfaces to their box triangles */ - box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */ - box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */ - box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ - box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */ - box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */ - box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */ - - /* Map the interfaces to their square segments */ - square_interfaces[0] = interf_adiabatic; /* Bottom */ - square_interfaces[1] = interf_Tb; /* Lef */ - square_interfaces[2] = interf_adiabatic; /* Top */ - square_interfaces[3] = interf_H; /* Right */ - - /* Create the box scene */ - CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices, - box_get_interface, box_nvertices, box_get_position, box_interfaces, - &box_scn) == RES_OK); - - /* Create the square scene */ - CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, - square_get_interface, square_nvertices, square_get_position, - square_interfaces, &square_scn) == RES_OK); - - /* Release the interfaces */ - CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK); - CHK(sdis_interface_ref_put(interf_Tb) == RES_OK); - CHK(sdis_interface_ref_put(interf_H) == RES_OK); - - uv[0] = 0.3; - uv[1] = 0.3; - iprim = 6; - - #define SOLVE sdis_solve_probe_boundary - #define F SDIS_FRONT - CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG); - CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); - - ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); - - CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); - CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); - CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK); - CHK(sdis_estimator_ref_put(estimator) == RES_OK); - CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK); - printf("Boundary temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n", - SPLIT3(pos), ref, T.E, T.SE); - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, 3*T.SE)); - - uv[0] = 0.5; - iprim = 3; - CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); - CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); - CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); - CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK); - CHK(sdis_estimator_ref_put(estimator) == RES_OK); - CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK); - printf("Boundary temperature of the square at (%g %g) = %g ~ %g +/- %g\n", - SPLIT2(pos), ref, T.E, T.SE); - printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); - CHK(nfails + nreals == N); - CHK(nfails < N/1000); - CHK(eq_eps(T.E, ref, 3*T.SE)); - #undef SOLVE - - CHK(sdis_scene_ref_put(box_scn) == RES_OK); - CHK(sdis_scene_ref_put(square_scn) == RES_OK); - CHK(sdis_device_ref_put(dev) == RES_OK); - - check_memory_allocator(&allocator); - mem_shutdown_proxy_allocator(&allocator); - CHK(mem_allocated_size() == 0); - return 0; -} -