commit 4c159e942201e107b30ab330feb2769dc5f3c24d parent 9d73d988cf6bd9f7204de36f6542d3afeeaeb85a Author: Vincent Forest <vincent.forest@meso-star.com> Date: Thu, 24 Jun 2021 10:08:29 +0200 Merge branch 'release_0.12' Diffstat:
78 files changed, 2288 insertions(+), 170 deletions(-)
diff --git a/README.md b/README.md @@ -111,10 +111,17 @@ variable the install directories of its dependencies. ## Release notes +### Version 0.12 + +Add the support of thermal contact resistance between two solids: the new +`thermal_contact_resistance` functor on the data structure `struct +sdis_interface_shader` defines the thermal resistance contact in K.m^2.W^-1 at +a given time and at a specific position onto the interface. + ### 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 + then be used to quickly evaluate the system at the same time but 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. @@ -350,7 +357,7 @@ First version and implementation of the Stardis-Solver API. ## License -Copyright (C) 2016-2020 |Meso|Star> (<contact@meso-star.com>). Stardis-Solver +Copyright (C) 2016-2021 |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 @@ -1,4 +1,4 @@ -# Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +# Copyright (C) 2016-2021 |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 @@ -13,7 +13,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. -cmake_minimum_required(VERSION 3.0) +cmake_minimum_required(VERSION 3.1) project(stardis C) enable_testing() @@ -29,12 +29,12 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS # Check dependencies ############################################################################### find_package(RCMake 0.4 REQUIRED) -find_package(Star2D 0.4 REQUIRED) -find_package(Star3D 0.7 REQUIRED) +find_package(Star2D 0.5 REQUIRED) +find_package(Star3D 0.8 REQUIRED) find_package(StarSP 0.8 REQUIRED) find_package(StarEnc2D 0.5 REQUIRED) find_package(StarEnc3D 0.5 REQUIRED) -find_package(RSys 0.10 REQUIRED) +find_package(RSys 0.12 REQUIRED) find_package(OpenMP 2.0 REQUIRED) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR}) @@ -56,7 +56,7 @@ rcmake_append_runtime_dirs(_runtime_dirs # Configure and define targets ############################################################################### set(VERSION_MAJOR 0) -set(VERSION_MINOR 11) +set(VERSION_MINOR 12) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) @@ -178,6 +178,8 @@ if(NOT NO_TEST) new_test(test_sdis_compute_power) new_test(test_sdis_conducto_radiative) new_test(test_sdis_conducto_radiative_2d) + new_test(test_sdis_contact_resistance) + new_test(test_sdis_contact_resistance_2) new_test(test_sdis_convection) new_test(test_sdis_convection_non_uniform) new_test(test_sdis_data) @@ -199,6 +201,7 @@ if(NOT NO_TEST) new_test(test_sdis_solve_medium) new_test(test_sdis_solve_medium_2d) new_test(test_sdis_transcient) + new_test(test_sdis_unstationary_atm) new_test(test_sdis_volumic_power) new_test(test_sdis_volumic_power4) diff --git a/src/sdis.h b/src/sdis.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -221,11 +221,16 @@ struct sdis_interface_shader { * [0 convection_coef_upper_bound] */ double convection_coef_upper_bound; + /* May be NULL for solid/fluid or if the thermal contact resistance is 0 onto + * the whole interface. */ + sdis_interface_getter_T thermal_contact_resistance; /* In K.m^2.W^-1 */ + struct sdis_interface_side_shader front; struct sdis_interface_side_shader back; }; #define SDIS_INTERFACE_SHADER_NULL__ \ - {NULL, 0, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__} + {NULL, 0, NULL, SDIS_INTERFACE_SIDE_SHADER_NULL__, \ + SDIS_INTERFACE_SIDE_SHADER_NULL__} static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL = SDIS_INTERFACE_SHADER_NULL__; diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_Xd_end.h b/src/sdis_Xd_end.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_camera.c b/src/sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_camera.h b/src/sdis_camera.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_data.c b/src/sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_device.c b/src/sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -1191,6 +1191,7 @@ sdis_green_path_get_limit_point pt->type = SDIS_VERTEX; break; case SDIS_GREEN_PATH_END_RADIATIVE: + case SDIS_GREEN_PATH_END_ERROR: res = RES_BAD_OP; goto error; break; diff --git a/src/sdis_green.h b/src/sdis_green.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_heat_path.c b/src/sdis_heat_path.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -478,6 +478,7 @@ XD(solid_solid_boundary_path) double tmp; double r; double power; + double tcr; float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM]; float dir_front[DIM], dir_back[DIM]; float* dir; @@ -501,6 +502,9 @@ XD(solid_solid_boundary_path) ASSERT(solid_front->type == SDIS_SOLID); ASSERT(solid_back->type == SDIS_SOLID); + /* Retrieve the thermal contact resistance */ + tcr = interface_get_thermal_contact_resistance(interf, frag); + /* Fetch the properties of the media */ lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx); lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx); @@ -560,19 +564,44 @@ XD(solid_solid_boundary_path) goto error; } - /* Define the reinjection side. Note that the proba should be : - * Lf/Df' / (Lf/Df' + Lb/Db') - * - * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted - * delta of the <front|back> side, i.e. : - * D<f|b>' = reinject_dst_<front|back> / sqrt(DIM) - * - * Anyway, one can avoid to compute the adjusted delta by directly using the - * adjusted reinjection distance since the resulting proba is strictly the - * same; sqrt(DIM) can be simplified. */ r = ssp_rng_canonical(rng); - proba = (lambda_front/reinject_dst_front) - / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back); + if(tcr == 0) { /* No thermal contact resistance */ + /* Define the reinjection side. Note that the proba should be : Lf/Df' / + * (Lf/Df' + Lb/Db') + * + * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted + * delta of the <front|back> side, i.e. : D<f|b>' = + * reinject_dst_<front|back> / sqrt(DIM) + * + * Anyway, one can avoid to compute the adjusted delta by directly using the + * adjusted reinjection distance since the resulting proba is strictly the + * same; sqrt(DIM) can be simplified. */ + proba = (lambda_front/reinject_dst_front) + / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back); + } else { + const double df = reinject_dst_front/sqrt(DIM); + const double db = reinject_dst_back/sqrt(DIM); + const double tmp_front = lambda_front/df; + const double tmp_back = lambda_back/db; + const double tmp_r = tcr*tmp_front*tmp_back; + switch(rwalk->hit_side) { + case SDIS_BACK: + /* When coming from the BACK side, the probability to be reinjected on + * the FRONT side depends on the thermal contact resistance: it + * decreases when the TCR increases (and tends to 0 when TCR -> +inf) */ + proba = (tmp_front) / (tmp_front + tmp_back + tmp_r); + break; + case SDIS_FRONT: + /* Same thing when coming from the FRONT side: the probability of + * reinjection on the FRONT side depends on the thermal contact + * resistance: it increases when the TCR increases (and tends to 1 when + * the TCR -> +inf) */ + proba = (tmp_front + tmp_r) / (tmp_front + tmp_back + tmp_r); + break; + default: FATAL("Unreachable code.\n"); break; + } + } + if(r < proba) { /* Reinject in front */ dir = dir_front; hit = &hit0; diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -171,15 +171,6 @@ XD(convective_path) /* Cannot be in the fluid without starting there. */ ASSERT(path_started_in_fluid); - if(ctx->green_path) { - log_err(scn->dev, - "%s: the upper bound of the convection cannot of an enclosure cannot be " - "null when registering the green function; initial condition is not " - "supported.\n", FUNC_NAME); - res = RES_BAD_ARG; - goto error; - } - rwalk->vtx.time = fluid_get_t0(rwalk->mdm); tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx); if(tmp >= 0) { diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -56,7 +56,6 @@ check_interface_shader " shader's pointer function for this attribute should be NULL.\n", caller_name); } - if(shader->convection_coef_upper_bound < 0) { log_warn(dev, "%s: Invalid upper bound for convection coefficient (%g).\n", @@ -64,6 +63,14 @@ check_interface_shader if(type[0] == SDIS_FLUID || type[1] == SDIS_FLUID) return 0; } + if((type[0] != SDIS_SOLID || type[1] != SDIS_SOLID) + && shader->thermal_contact_resistance) { + log_warn(dev, + "%s: only solid/solid interface can have a thermal contact resistance. The " + " shader's pointer function for this attribute should be NULL.\n", + caller_name); + } + FOR_EACH(i, 0, 2) { switch(type[i]) { case SDIS_SOLID: diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -89,6 +89,16 @@ interface_get_convection_coef } static INLINE double +interface_get_thermal_contact_resistance + (const struct sdis_interface* interf, + const struct sdis_interface_fragment* frag) +{ + ASSERT(interf && frag); + return interf->shader.thermal_contact_resistance + ? interf->shader.thermal_contact_resistance(frag, interf->data) : 0; +} + +static INLINE double interface_get_convection_coef_upper_bound (const struct sdis_interface* interf) { diff --git a/src/sdis_log.c b/src/sdis_log.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_log.h b/src/sdis_log.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_medium.c b/src/sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_misc.c b/src/sdis_misc.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -20,6 +20,15 @@ #include <rsys/float3.h> #include <star/ssp.h> +struct bound_flux_result { + double Tradiative; + double Tboundary; + double Tfluid; +}; +#define BOUND_FLUX_RESULT_NULL__ {0,0,0} +static const struct bound_flux_result +BOUND_FLUX_RESULT_NULL = BOUND_FLUX_RESULT_NULL__; + struct accum { double sum; /* Sum of MC weights */ double sum2; /* Sum of square MC weights */ diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -65,7 +65,7 @@ XD(time_rewind) /* Fetch initial temperature */ temperature = solid_get_temperature(mdm, &rwalk->vtx); if(temperature < 0) { - log_err(mdm->dev, "%s: the path reaches the limit condition by the " + log_err(mdm->dev, "%s: the path reaches the limit condition but the " "temperature remains unknown.\n", FUNC_NAME); res = RES_BAD_ARG; goto error; diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -23,8 +23,10 @@ /* Forward declarations */ struct green_path_handle; +struct sdis_heat_path; struct sdis_scene; struct ssp_rng; +struct bound_flux_result; enum flux_flag { FLUX_FLAG_CONVECTIVE = BIT(FLUX_CONVECTIVE), @@ -95,7 +97,7 @@ boundary_flux_realisation_2d const double time, const enum sdis_side solid_side, const int flux_mask, /* Combination of enum flux_flag */ - double weight[FLUX_NAMES_COUNT__]); + struct bound_flux_result* result); extern LOCAL_SYM res_T boundary_flux_realisation_3d @@ -106,7 +108,7 @@ boundary_flux_realisation_3d const double time, const enum sdis_side solid_side, const int flux_mask, /* Combination of enum flux_flag */ - double weight[FLUX_NAMES_COUNT__]); + struct bound_flux_result* result); /******************************************************************************* * Realisation along a given ray at a given time. Available only in 3D. diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -158,8 +158,7 @@ XD(probe_realisation) res = register_heat_vertex(heat_path, &rwalk.vtx, 0, type); if(res != RES_OK) goto error; - /* No initial condition with green */ - if(!green_path && t0 >= rwalk.vtx.time) { + if(t0 >= rwalk.vtx.time) { double tmp; /* Check the initial condition. */ rwalk.vtx.time = t0; @@ -284,7 +283,7 @@ XD(boundary_flux_realisation) const double time, const enum sdis_side solid_side, const int flux_mask, - double weight[3]) + struct bound_flux_result* result) { struct rwalk_context ctx = RWALK_CONTEXT_NULL; struct XD(rwalk) rwalk; @@ -308,7 +307,7 @@ XD(boundary_flux_realisation) 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 && weight && time >= 0 ); + ASSERT(uv && result && time >= 0 ); #if SDIS_XD_DIMENSION == 2 #define SET_PARAM(Dest, Src) (Dest).u = (Src); @@ -349,7 +348,7 @@ XD(boundary_flux_realisation) T.func = XD(boundary_path); res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; - weight[0] = T.value; + result->Tboundary = T.value; /* Fetch the fluid medium */ interf = scene_get_interface(scn, (unsigned)iprim); @@ -361,7 +360,8 @@ XD(boundary_flux_realisation) T.func = XD(radiative_path); res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; - weight[1] = T.value; + ASSERT(T.value >= 0); + result->Tradiative = T.value; } /* Compute fluid temperature */ @@ -370,7 +370,7 @@ XD(boundary_flux_realisation) T.func = XD(convective_path); res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T); if(res != RES_OK) return res; - weight[2] = T.value; + result->Tfluid = T.value; } #undef SET_PARAM diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -524,8 +524,7 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash) } #endif - res = hash_sha256(scn->dev->allocator, data, len, hash); - if(res != RES_OK) goto error; + hash_sha256(data, len, hash); exit: #ifdef COMPILER_GCC diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -24,6 +24,9 @@ #include "sdis_scene_c.h" #include <star/ssp.h> +#include <star/senc2d.h> +#include <star/senc3d.h> + #include <rsys/cstr.h> #include <rsys/float22.h> #include <rsys/float33.h> @@ -442,17 +445,21 @@ XD(hit_filter_function) (const struct sXd(hit)* hit, const float org[DIM], const float dir[DIM], + const float range[2], void* ray_data, void* global_data) { const struct hit_filter_data* filter_data = ray_data; const struct sXd(hit)* hit_from = &filter_data->XD(hit); - (void)org, (void)dir, (void)global_data; + (void)org, (void)dir, (void)global_data, (void)range; if(!ray_data || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */ if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1; + /* No displacement => assume self intersection */ + if(hit->distance <= 0) return 1; + if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) { float pos[DIM]; fX(add)(pos, org, fX(mulf)(pos, dir, hit->distance)); @@ -547,7 +554,7 @@ XD(run_analyze) { struct geometry geom; struct sencXd(device)* senc = NULL; - struct sencXd(scene)* senc3d_scn = NULL; + struct sencXd(scene)* senc_scn = NULL; unsigned count; res_T res = RES_OK; ASSERT(scn && nprims && indices && interf && nverts && position && out_scn); @@ -564,21 +571,37 @@ XD(run_analyze) res = sencXd(scene_create)(senc, SENCXD_(CONVENTION_NORMAL_BACK) | SENCXD_(CONVENTION_NORMAL_OUTSIDE), (unsigned)nprims, XD(geometry_indices), geometry_media, - (unsigned)nverts, XD(geometry_position), &geom, &senc3d_scn); + (unsigned)nverts, XD(geometry_position), &geom, &senc_scn); if(res != RES_OK) goto error; /* With il-formed scenes, scene creation can success without being able * to extract enclosures; in this case just fail */ - res = sencXd(scene_get_enclosure_count(senc3d_scn, &count)); - if(res != RES_OK) goto error; + res = sencXd(scene_get_enclosure_count(senc_scn, &count)); + if(res != RES_OK) { + count = 0; +#if DIM==2 + senc2d_scene_get_overlapping_segments_count(senc_scn, &count); + if(count > 0) + log_err(scn->dev, + "%s: the scene includes overlapping segments.\n", + FUNC_NAME); +#else + senc3d_scene_get_overlapping_triangles_count(senc_scn, &count); + if(count > 0) + log_err(scn->dev, + "%s: the scene includes overlapping triangles.\n", + FUNC_NAME); +#endif + goto error; + } exit: if(senc) SENCXD(device_ref_put(senc)); - if(out_scn) *out_scn = senc3d_scn; + if(out_scn) *out_scn = senc_scn; return res; error: - if(senc3d_scn) { - SENCXD(scene_ref_put(senc3d_scn)); - senc3d_scn = NULL; + if(senc_scn) { + SENCXD(scene_ref_put(senc_scn)); + senc_scn = NULL; } goto exit; } @@ -589,17 +612,17 @@ error: static res_T XD(setup_properties) (struct sdis_scene* scn, - struct sencXd(scene)* senc3d_scn, + struct sencXd(scene)* senc_scn, void (*interf)(const size_t itri, struct sdis_interface**, void*), void* ctx) { unsigned iprim, nprims; res_T res = RES_OK; - ASSERT(scn && senc3d_scn && interf); + ASSERT(scn && senc_scn && interf); clear_properties(scn); - SENCXD(scene_get_primitives_count(senc3d_scn, &nprims)); + SENCXD(scene_get_primitives_count(senc_scn, &nprims)); FOR_EACH(iprim, 0, nprims) { struct prim_prop* prim_prop; struct sdis_interface* itface; @@ -610,7 +633,7 @@ XD(setup_properties) size_t ninterfaces; /* Fetch the enclosures that the segment/triangle splits */ - SENCXD(scene_get_primitive_enclosures(senc3d_scn, iprim, enclosures)); + SENCXD(scene_get_primitive_enclosures(senc_scn, iprim, enclosures)); /* Fetch the interface of the primitive */ interf(iprim, &itface, ctx); @@ -666,7 +689,7 @@ error: /* Build the Star-XD scene view of the whole scene */ static res_T -XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_scn) +XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(scene)* senc_scn) { struct sXd(device)* sXd_dev = NULL; struct sXd(shape)* sXd_shape = NULL; @@ -674,9 +697,9 @@ XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_sc struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL; unsigned nprims, nverts; res_T res = RES_OK; - ASSERT(scn && senc3d_scn); + ASSERT(scn && senc_scn); - SENCXD(scene_get_vertices_count(senc3d_scn, &nverts)); + SENCXD(scene_get_vertices_count(senc_scn, &nverts)); /* Setup the vertex data */ vdata.usage = SXD_POSITION; @@ -686,18 +709,18 @@ XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_sc /* Create the Star-XD geometry of the whole scene */ #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0 sXd_dev = scn->dev->sXd(dev); - SENCXD(scene_get_primitives_count(senc3d_scn, &nprims)); + SENCXD(scene_get_primitives_count(senc_scn, &nprims)); #if DIM == 2 CALL(sXd(shape_create_line_segments)(sXd_dev, &sXd_shape)); CALL(sXd(line_segments_set_hit_filter_function)(sXd_shape, XD(hit_filter_function), NULL)); CALL(sXd(line_segments_setup_indexed_vertices)(sXd_shape, nprims, - XD(scene_indices), nverts, &vdata, 1, senc3d_scn)); + XD(scene_indices), nverts, &vdata, 1, senc_scn)); #else CALL(sXd(shape_create_mesh)(sXd_dev, &sXd_shape)); CALL(sXd(mesh_set_hit_filter_function)(sXd_shape, XD(hit_filter_function), NULL)); CALL(sXd(mesh_setup_indexed_vertices)(sXd_shape, nprims, XD(scene_indices), - nverts, &vdata, 1, senc3d_scn)); + nverts, &vdata, 1, senc_scn)); #endif CALL(sXd(scene_create)(sXd_dev, &sXd_scn)); CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape)); @@ -814,19 +837,19 @@ error: /* Build the Star-XD scene view and define its associated data of the finite * fluid enclosures */ static res_T -XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc3d_scn) +XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc_scn) { struct sencXd(enclosure)* enc = NULL; unsigned ienc, nencs; int inner_multi = 0; res_T res = RES_OK; - ASSERT(scn && senc3d_scn); + ASSERT(scn && senc_scn); - SENCXD(scene_get_enclosure_count(senc3d_scn, &nencs)); + SENCXD(scene_get_enclosure_count(senc_scn, &nencs)); FOR_EACH(ienc, 0, nencs) { struct sencXd(enclosure_header) header; - SENCXD(scene_get_enclosure(senc3d_scn, ienc, &enc)); + SENCXD(scene_get_enclosure(senc_scn, ienc, &enc)); SENCXD(enclosure_get_header(enc, &header)); if(header.is_infinite) { @@ -868,7 +891,7 @@ XD(scene_create) const struct sdis_scene_create_args* args, struct sdis_scene** out_scn) { - struct sencXd(scene)* senc3d_scn = NULL; + struct sencXd(scene)* senc_scn = NULL; struct sdis_scene* scn = NULL; res_T res = RES_OK; @@ -905,34 +928,34 @@ XD(scene_create) args->nvertices, args->get_position, args->context, - &senc3d_scn); + &senc_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, args->get_interface, args->context); + res = XD(setup_properties)(scn, senc_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); goto error; } - res = XD(setup_scene_geometry)(scn, senc3d_scn); + res = XD(setup_scene_geometry)(scn, senc_scn); if(res != RES_OK) { log_err(dev, "%s: could not setup the scene geometry.\n", FUNC_NAME); goto error; } - res = XD(setup_enclosures)(scn, senc3d_scn); + res = XD(setup_enclosures)(scn, senc_scn); if(res != RES_OK) { log_err(dev, "%s: could not setup the enclosures.\n", FUNC_NAME); goto error; } - scn->sencXd(scn) = senc3d_scn; + scn->sencXd(scn) = senc_scn; exit: if(out_scn) *out_scn = scn; return res; error: - if(senc3d_scn) SENCXD(scene_ref_put(senc3d_scn)); + if(senc_scn) SENCXD(scene_ref_put(senc_scn)); if(scn) { SDIS(scene_ref_put(scn)); scn = NULL; @@ -1110,13 +1133,22 @@ XD(scene_get_medium) res = RES_BAD_OP; goto error; } - + +#if DIM == 2 + if(iprim > 10 && iprim > (size_t)((double)nprims * 0.05)) { + log_warn(scn->dev, + "%s: performance issue. Up to %lu primitives were tested to define the " + "current medium at {%g, %g}.\n", + FUNC_NAME, (unsigned long)iprim, SPLIT2(P)); + } +#else if(iprim > 10 && iprim > (size_t)((double)nprims * 0.05)) { log_warn(scn->dev, "%s: performance issue. Up to %lu primitives were tested to define the " "current medium at {%g, %g, %g}.\n", FUNC_NAME, (unsigned long)iprim, SPLIT3(P)); } +#endif exit: *out_medium = medium; diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -572,7 +572,7 @@ XD(solve_boundary_flux) const struct sdis_interface* interf; const struct sdis_medium *fmd, *bmd; enum sdis_side solid_side, fluid_side; - double T_brf[3] = { 0, 0, 0 }; + struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__; const double Tref = scn->reference_temperature; double epsilon, hc, hr, imposed_flux, imposed_temp; size_t iprim; @@ -656,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, flux_mask, T_brf); + solid_side, flux_mask, &result); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); @@ -666,16 +666,14 @@ XD(solve_boundary_flux) continue; } else if(res_simul == RES_OK) { /* Update accumulators */ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; - const double Tboundary = T_brf[0]; - const double Tradiative = T_brf[1]; - const double Tfluid = T_brf[2]; - const double w_conv = hc * (Tboundary - Tfluid); - const double w_rad = hr * (Tboundary - Tradiative); + const double w_conv = hc * (result.Tboundary - result.Tfluid); + const double w_rad = (result.Tradiative < 0) ? + 0 : hr * (result.Tboundary - result.Tradiative); const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0; const double w_total = w_conv + w_rad + w_imp; /* Temperature */ - acc_temp->sum += Tboundary; - acc_temp->sum2 += Tboundary*Tboundary; + acc_temp->sum += result.Tboundary; + acc_temp->sum2 += result.Tboundary*result.Tboundary; ++acc_temp->count; /* Time */ acc_time->sum += usec; diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -413,8 +413,8 @@ XD(solve_probe_boundary_flux) log_err(scn->dev, "%s: Attempt to compute a flux at a %s-%s interface.\n", FUNC_NAME, - (fmd->type == SDIS_FLUID ? "fluid" : "solid"), - (bmd->type == SDIS_FLUID ? "fluid" : "solid")); + (!fmd ? "undefined" : (fmd->type == SDIS_FLUID ? "fluid" : "solid")), + (!bmd ? "undefined" : (bmd->type == SDIS_FLUID ? "fluid" : "solid"))); res = RES_BAD_ARG; goto error; } @@ -479,7 +479,7 @@ XD(solve_probe_boundary_flux) struct accum* acc_fimp = &acc_fi[ithread]; double time, epsilon, hc, hr, imposed_flux, imposed_temp; int flux_mask = 0; - double T_brf[3] = { 0, 0, 0 }; + struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__; const double Tref = scn->reference_temperature; size_t n; int pcent; @@ -514,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, flux_mask, T_brf); + time, solid_side, flux_mask, &result); /* Stop time registration */ time_sub(&t0, time_current(&t1), &t0); @@ -524,16 +524,14 @@ XD(solve_probe_boundary_flux) continue; } else if(res_simul == RES_OK) { /* Update accumulators */ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; - const double Tboundary = T_brf[0]; - const double Tradiative = T_brf[1]; - const double Tfluid = T_brf[2]; - const double w_conv = hc * (Tboundary - Tfluid); - const double w_rad = hr * (Tboundary - Tradiative); + const double w_conv = hc * (result.Tboundary - result.Tfluid); + const double w_rad = (result.Tradiative < 0) ? + 0 : hr * (result.Tboundary - result.Tradiative); const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0; const double w_total = w_conv + w_rad + w_imp; /* Temperature */ - acc_temp->sum += Tboundary; - acc_temp->sum2 += Tboundary*Tboundary; + acc_temp->sum += result.Tboundary; + acc_temp->sum2 += result.Tboundary*result.Tboundary; ++acc_temp->count; /* Time */ acc_time->sum += usec; diff --git a/src/test_sdis_accum_buffer.c b/src/test_sdis_accum_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_camera.c b/src/test_sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_compute_power.c b/src/test_sdis_compute_power.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_contact_resistance.c b/src/test_sdis_contact_resistance.c @@ -0,0 +1,444 @@ +/* Copyright (C) 2016-2021 |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 "test_sdis_contact_resistance.h" + +#include <rsys/clock_time.h> +#include <rsys/mem_allocator.h> +#include <rsys/double3.h> +#include <rsys/math.h> +#include <star/ssp.h> + +/* + * The scene is composed of a solid cube/square of size L whose temperature is + * unknown. The cube is made of 2 solids that meet at x=e in ]0 L[ and the + * thermal contact resistance is R, thus T(X0-) differs from T(X0+). + * The faces are adiabatic exept at x=0 where T(0)=T0 and at x=L where T(L)=TL. + * At steady state: + * + * T(X0-) = (T0 * (R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0)) + * + TL * (R * LAMBDA2 / (L - X0))) + * / ((1 + R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0)) - 1) + * + * T(X0+) = T(X0-) * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0 + * + * T(x) is linear between T(0) and T(X0-) if x in [0 X0[ + * T(x) is linear between T(X0+) and T(L) if x in ]X0 L] + * + * 3D 2D + * + * /////////(L,L,L) /////////(L,L) + * +-------+ +-------+ + * /' / /| | ! | + * +-------+ TL T0 r TL + * | | ! | | | ! | + * T0 +.r...|.+ +-------+ + * |, ! |/ (0,0)///x=X0/// + * +-------+ + * (0,0,0)///x=X0/// + */ + +#define UNKNOWN_TEMPERATURE -1 +#define N 10000 /* #realisations */ + +#define T0 0.0 +#define LAMBDA1 0.1 + +#define TL 100.0 +#define LAMBDA2 0.2 + +#define DELTA1 X0/25.0 +#define DELTA2 (L-X0)/25.0 + +/******************************************************************************* + * Media + ******************************************************************************/ +struct solid { + double lambda; + double rho; + double cp; + double delta; +}; + +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx); + return UNKNOWN_TEMPERATURE; +} + +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->cp; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->lambda; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->rho; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->delta; +} + +static double +solid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return UNKNOWN_TEMPERATURE; +} + +/******************************************************************************* + * Interfaces + ******************************************************************************/ +struct interf { + double temperature; + double resistance; +}; + +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) +{ + CHK(frag && data); + return 0; +} + +static double +interface_get_contact_resistance + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf = sdis_data_cget(data); + CHK(frag && data); + return interf->resistance; +} + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +solve + (struct sdis_scene* scn, + struct interf* interf_props, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + double ref_L, ref_R; + enum sdis_scene_dimension dim; + const int nsimuls = 8; + int isimul; + ASSERT(scn && interf_props && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + FOR_EACH(isimul, 0, nsimuls) { + double x, ref; + double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); + + interf_props->resistance = r; + + ref_L = ( + T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) + + TL * (r * LAMBDA2 / (L - X0)) + ) + / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); + + ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; + + if(isimul % 2) { /* In solid 1 */ + x = ssp_rng_uniform_double(rng, 0.05 * X0, 0.95 * X0); + ref = T0 * (1 - x / X0) + ref_L * x / X0; + } else { /* In solid 2 */ + x = X0 + ssp_rng_uniform_double(rng, 0.05 * (L - X0), 0.95 * (L - X0)); + ref = ref_R * (1 - (x - X0) / (L - X0)) + TL * (x - X0) / (L - X0); + } + + solve_args.position[0] = x; + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L); + solve_args.position[2] = (dim == SDIS_SCENE_2D) + ? 0 : ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L); + + solve_args.nrealisations = N; + solve_args.time_range[0] = 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) { + case SDIS_SCENE_2D: + printf("Steady temperature at (%g, %g) with R=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), r, ref, T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Steady temperature at (%g, %g, %g) with R=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), r, 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); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + CHK(eq_eps(T.E, ref, T.SE * 3)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +/******************************************************************************* + * 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* solid1 = NULL; + struct sdis_medium* solid2 = NULL; + struct sdis_interface* interf_adiabatic1 = NULL; + struct sdis_interface* interf_adiabatic2 = NULL; + struct sdis_interface* interf_T0 = NULL; + struct sdis_interface* interf_TL = NULL; + struct sdis_interface* interf_R = 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* model3d_interfaces[22 /*#triangles*/]; + struct sdis_interface* model2d_interfaces[7/*#segments*/]; + struct interf* interf_props = NULL; + struct solid* solid_props = NULL; + struct ssp_rng* rng = NULL; + (void)argc, (void)argv; + + OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); + OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev)); + + fluid_shader.temperature = fluid_get_temperature; + OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + + /* Setup the solid shader */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + + /* Create the solid medium #1 */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = LAMBDA1; + solid_props->cp = 2; + solid_props->rho = 25; + solid_props->delta = DELTA1; + OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); + OK(sdis_data_ref_put(data)); + + /* Create the solid medium #2 */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = LAMBDA2; + solid_props->cp = 2; + solid_props->rho = 25; + solid_props->delta = DELTA2; + OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); + OK(sdis_data_ref_put(data)); + + /* Setup the interface shader */ + interf_shader.front.temperature = interface_get_temperature; + interf_shader.back.temperature = interface_get_temperature; + interf_shader.convection_coef = interface_get_convection_coef; + + /* Create the adiabatic interfaces */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + OK(sdis_interface_create + (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1)); + OK(sdis_interface_create + (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2)); + OK(sdis_data_ref_put(data)); + + /* Create the T0 interface */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = T0; + OK(sdis_interface_create + (dev, solid1, fluid, &interf_shader, data, &interf_T0)); + OK(sdis_data_ref_put(data)); + + /* Create the TL interface */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = TL; + OK(sdis_interface_create + (dev, solid2, fluid, &interf_shader, data, &interf_TL)); + OK(sdis_data_ref_put(data)); + + /* Create the solid1-solid2 interface */ + interf_shader.convection_coef = NULL; + interf_shader.thermal_contact_resistance = interface_get_contact_resistance; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + OK(sdis_interface_create + (dev, solid1, solid2, &interf_shader, data, &interf_R)); + OK(sdis_data_ref_put(data)); + + /* Release the media */ + OK(sdis_medium_ref_put(solid1)); + OK(sdis_medium_ref_put(solid2)); + OK(sdis_medium_ref_put(fluid)); + + /* Map the interfaces to their box triangles */ + /* Front */ + model3d_interfaces[0] = interf_adiabatic1; + model3d_interfaces[1] = interf_adiabatic1; + model3d_interfaces[2] = interf_adiabatic2; + model3d_interfaces[3] = interf_adiabatic2; + /* Left */ + model3d_interfaces[4] = interf_T0; + model3d_interfaces[5] = interf_T0; + /* Back */ + model3d_interfaces[6] = interf_adiabatic1; + model3d_interfaces[7] = interf_adiabatic1; + model3d_interfaces[8] = interf_adiabatic2; + model3d_interfaces[9] = interf_adiabatic2; + /* Right */ + model3d_interfaces[10] = interf_TL; + model3d_interfaces[11] = interf_TL; + /* Top */ + model3d_interfaces[12] = interf_adiabatic1; + model3d_interfaces[13] = interf_adiabatic1; + model3d_interfaces[14] = interf_adiabatic2; + model3d_interfaces[15] = interf_adiabatic2; + /* Bottom */ + model3d_interfaces[16] = interf_adiabatic1; + model3d_interfaces[17] = interf_adiabatic1; + model3d_interfaces[18] = interf_adiabatic2; + model3d_interfaces[19] = interf_adiabatic2; + /* Inside */ + model3d_interfaces[20] = interf_R; + model3d_interfaces[21] = interf_R; + + /* Map the interfaces to their square segments */ + /* Bottom */ + model2d_interfaces[0] = interf_adiabatic2; + model2d_interfaces[1] = interf_adiabatic1; + /* Left */ + model2d_interfaces[2] = interf_T0; + /* Top */ + model2d_interfaces[3] = interf_adiabatic1; + model2d_interfaces[4] = interf_adiabatic2; + /* Right */ + model2d_interfaces[5] = interf_TL; + /* Contact */ + model2d_interfaces[6] = interf_R; + + /* Create the box scene */ + scn_args.get_indices = model3d_get_indices; + scn_args.get_interface = model3d_get_interface; + scn_args.get_position = model3d_get_position; + scn_args.nprimitives = model3d_ntriangles; + scn_args.nvertices = model3d_nvertices; + scn_args.context = model3d_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); + + /* Create the square scene */ + scn_args.get_indices = model2d_get_indices; + scn_args.get_interface = model2d_get_interface; + scn_args.get_position = model2d_get_position; + scn_args.nprimitives = model2d_nsegments; + scn_args.nvertices = model2d_nvertices; + scn_args.context = model2d_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); + + /* Release the interfaces */ + OK(sdis_interface_ref_put(interf_adiabatic1)); + OK(sdis_interface_ref_put(interf_adiabatic2)); + OK(sdis_interface_ref_put(interf_T0)); + OK(sdis_interface_ref_put(interf_TL)); + OK(sdis_interface_ref_put(interf_R)); + + /* Solve */ + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); + printf(">> Box scene\n"); + solve(box_scn, interf_props, rng); + printf("\n>> Square scene\n"); + solve(square_scn, interf_props, rng); + + 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); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sdis_contact_resistance.h b/src/test_sdis_contact_resistance.h @@ -0,0 +1,161 @@ +/* Copyright (C) 2016-2021 |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/>. */ + + +/* + * The scene is composed of a solid cube/square of size L. The cube/square + * is made of 2 solids that meet at x=e in ]0 L[. + * + * 3D 2D + * + * /////////(L,L,L) /////////(L,L) + * +-------+ +-------+ + * /' / /| | ! | + * +-------+ TL T0 r TL + * | | ! | | | ! | + * T0 +.r...|.+ +-------+ + * |, ! |/ (0,0)///x=X0/// + * +-------+ + * (0,0,0)///x=X0/// + */ + +#define L 4.0 +#define X0 3.0 + + /******************************************************************************* + * Box geometry + ******************************************************************************/ +static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/] += { + 0, 0, 0, + X0, 0, 0, + L, 0, 0, + 0, L, 0, + X0, L, 0, + L, L, 0, + 0, 0, L, + X0, 0, L, + L, 0, L, + 0, L, L, + X0, L, L, + L, L, L +}; +static const size_t model3d_nvertices = sizeof(model3d_vertices) / sizeof(double[3]); + +/* The following array lists the indices toward the 3D vertices of each + * triangle. + * ,3---,4---,5 ,3----4----5 ,4 + * ,' | ,' | ,'/| ,'/| \ | \ | ,'/| + * 9----10---11 / | 9' / | \ | \ | 10 / | Y + * |', |', | / ,2 | / ,0---,1---,2 | / ,1 | + * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X + * 6----7----8' 6----7'---8' 7 / + * Front, right Back, left and Internal Z + * and Top faces bottom faces face */ +static const size_t model3d_indices[22/*#triangles*/ * 3/*#indices per triangle*/] += { + 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */ + 0, 6, 3, 3, 6, 9, /* -X */ + 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */ + 5, 11, 8, 8, 2, 5, /* +X */ + 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */ + 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */ + 4, 10, 7, 7, 1, 4 /* Inside */ +}; +static const size_t model3d_ntriangles = sizeof(model3d_indices) / sizeof(size_t[3]); + +static INLINE void +model3d_get_indices(const size_t itri, size_t ids[3], void* context) +{ + (void)context; + CHK(ids); + CHK(itri < model3d_ntriangles); + ids[0] = model3d_indices[itri * 3 + 0]; + ids[1] = model3d_indices[itri * 3 + 1]; + ids[2] = model3d_indices[itri * 3 + 2]; +} + +static INLINE void +model3d_get_position(const size_t ivert, double pos[3], void* context) +{ + (void)context; + CHK(pos); + CHK(ivert < model3d_nvertices); + pos[0] = model3d_vertices[ivert * 3 + 0]; + pos[1] = model3d_vertices[ivert * 3 + 1]; + pos[2] = model3d_vertices[ivert * 3 + 2]; +} + +static INLINE void +model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context) +{ + struct sdis_interface** interfaces = context; + CHK(context && bound); + CHK(itri < model3d_ntriangles); + *bound = interfaces[itri]; +} + +/******************************************************************************* + * Square geometry + ******************************************************************************/ +static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = { + L, 0, + X0, 0, + 0, 0, + 0, L, + X0, L, + L, L +}; +static const size_t model2d_nvertices = sizeof(model2d_vertices) / sizeof(double[2]); + +static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = { + 0, 1, 1, 2, /* Bottom */ + 2, 3, /* Left */ + 3, 4, 4, 5, /* Top */ + 5, 0, /* Right */ + 4, 1 /* Inside */ +}; +static const size_t model2d_nsegments = sizeof(model2d_indices) / sizeof(size_t[2]); + + +static INLINE void +model2d_get_indices(const size_t iseg, size_t ids[2], void* context) +{ + (void)context; + CHK(ids); + CHK(iseg < model2d_nsegments); + ids[0] = model2d_indices[iseg * 2 + 0]; + ids[1] = model2d_indices[iseg * 2 + 1]; +} + +static INLINE void +model2d_get_position(const size_t ivert, double pos[2], void* context) +{ + (void)context; + CHK(pos); + CHK(ivert < model2d_nvertices); + pos[0] = model2d_vertices[ivert * 2 + 0]; + pos[1] = model2d_vertices[ivert * 2 + 1]; +} + +static INLINE void +model2d_get_interface +(const size_t iseg, struct sdis_interface** bound, void* context) +{ + struct sdis_interface** interfaces = context; + CHK(context && bound); + CHK(iseg < model2d_nsegments); + *bound = interfaces[iseg]; +} diff --git a/src/test_sdis_contact_resistance_2.c b/src/test_sdis_contact_resistance_2.c @@ -0,0 +1,540 @@ +/* Copyright (C) 2016-2021 |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 "test_sdis_contact_resistance.h" + +#include <rsys/clock_time.h> +#include <rsys/mem_allocator.h> +#include <rsys/double3.h> +#include <rsys/math.h> +#include <star/ssp.h> + +/* + * The scene is composed of a solid cube/square of size L whose temperature is + * unknown. The cube is made of 2 solids that meet at x=e in ]0 L[ and the + * thermal contact resistance is R, thus T(X0-) differs from T(X0+). + * The faces are adiabatic exept at x=0 where T(0)=T0 and at x=L where T(L)=TL. + * At steady state: + * + * Flux(x0) = (T(x0+) - T(x0-)) / R + * + * 3D 2D + * + * /////////(L,L,L) /////////(L,L) + * +-------+ +-------+ + * /' / /| | ! | + * +-------+ TL T0 r TL + * | | ! | | | ! | + * T0 +.r...|.+ +-------+ + * |, ! |/ (0,0)///x=X0/// + * +-------+ + * (0,0,0)///x=X0/// + */ + +#define UNKNOWN_TEMPERATURE -1 +#define N 10000 /* #realisations */ + +#define T0 0.0 +#define LAMBDA1 0.1 + +#define TL 100.0 +#define LAMBDA2 0.2 + +#define DELTA1 X0/25.0 +#define DELTA2 (L-X0)/25.0 + +/******************************************************************************* + * Media + ******************************************************************************/ +struct solid { + double lambda; + double rho; + double cp; + double delta; +}; + +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx); + return UNKNOWN_TEMPERATURE; +} + +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->cp; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->lambda; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->rho; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return ((struct solid*)sdis_data_cget(data))->delta; +} + +static double +solid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + CHK(vtx && data); + return UNKNOWN_TEMPERATURE; +} + +/******************************************************************************* + * Interfaces + ******************************************************************************/ +struct interf { + double temperature; + double resistance; +}; + +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) +{ + CHK(frag && data); + return 0; +} + +static double +interface_get_contact_resistance + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf = sdis_data_cget(data); + CHK(frag && data); + return interf->resistance; +} + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +solve_probe + (struct sdis_scene* scn, + struct interf* interf_props, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_boundary_args solve_args + = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + double ref_L, ref_R; + enum sdis_scene_dimension dim; + int nsimuls; + int isimul; + ASSERT(scn && interf_props && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4; + FOR_EACH(isimul, 0, nsimuls) { + double ref; + double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); + + interf_props->resistance = r; + + ref_L = ( + T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) + + TL * (r * LAMBDA2 / (L - X0)) + ) + / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); + + ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; + + if(dim == SDIS_SCENE_2D) + /* last segment */ + solve_args.iprim = model2d_nsegments - 1; + else + /* last 2 triangles */ + solve_args.iprim = model3d_ntriangles - ((isimul % 2) ? 2 : 1); + + solve_args.uv[0] = ssp_rng_canonical(rng); + solve_args.uv[1] = (dim == SDIS_SCENE_2D) + ? 0 : ssp_rng_uniform_double(rng, 0, 1 - solve_args.uv[0]); + + if(isimul < nsimuls / 2) { + solve_args.side = SDIS_FRONT; + ref = ref_L; + } else { + solve_args.side = SDIS_BACK; + ref = ref_R; + } + + solve_args.nrealisations = N; + solve_args.time_range[0] = solve_args.time_range[1] = INF; + + time_current(&t0); + OK(sdis_solve_probe_boundary(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("Steady temperature at (%lu/%s/%g) with R=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, + (solve_args.side == SDIS_FRONT ? "front" : "back"), + solve_args.uv[0], + r, ref, T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Steady temperature at (%lu/%s/%g,%g) with R=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, + (solve_args.side == SDIS_FRONT ? "front" : "back"), + SPLIT2(solve_args.uv), + r, 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); + + CHK(nfails + nreals == N); + CHK(nfails <= N/1000); + CHK(eq_eps(T.E, ref, T.SE * 3)); + + OK(sdis_estimator_ref_put(estimator)); + } +} +static void +solve + (struct sdis_scene* scn, + struct interf* interf_props, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_boundary_args solve_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + const enum sdis_side all_front[] = { SDIS_FRONT, SDIS_FRONT }; + const enum sdis_side all_back[] = { SDIS_BACK, SDIS_BACK }; + size_t plist[2]; + size_t nreals; + size_t nfails; + double ref_L, ref_R; + enum sdis_scene_dimension dim; + int nsimuls; + int isimul; + ASSERT(scn && interf_props && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4; + FOR_EACH(isimul, 0, nsimuls) { + double ref; + double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); + + interf_props->resistance = r; + + ref_L = ( + T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) + + TL * (r * LAMBDA2 / (L - X0)) + ) + / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); + + ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; + + if(dim == SDIS_SCENE_2D) { + /* last segment */ + solve_args.nprimitives = 1; + plist[0] = model2d_nsegments - 1; + } else { + /* last 2 triangles */ + solve_args.nprimitives = 2; + plist[0] = model3d_ntriangles - 2; + plist[1] = model3d_ntriangles - 1; + } + solve_args.primitives = plist; + + if(isimul < nsimuls / 2) { + solve_args.sides = all_front; + ref = ref_L; + } else { + solve_args.sides = all_back; + ref = ref_R; + } + + solve_args.nrealisations = N; + solve_args.time_range[0] = solve_args.time_range[1] = INF; + + time_current(&t0); + OK(sdis_solve_boundary(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)); + + printf("Steady temperature at the %s side with R=%g = %g ~ %g +/- %g\n", + (solve_args.sides[0] == SDIS_FRONT ? "front" : "back"), + r, ref, T.E, T.SE); + printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); + printf("Elapsed time = %s\n", dump); + printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); + + CHK(nfails + nreals == N); + CHK(nfails <= N / 1000); + CHK(eq_eps(T.E, ref, T.SE * 3)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +/******************************************************************************* + * 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* solid1 = NULL; + struct sdis_medium* solid2 = NULL; + struct sdis_interface* interf_adiabatic1 = NULL; + struct sdis_interface* interf_adiabatic2 = NULL; + struct sdis_interface* interf_T0 = NULL; + struct sdis_interface* interf_TL = NULL; + struct sdis_interface* interf_R = 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* model3d_interfaces[22 /*#triangles*/]; + struct sdis_interface* model2d_interfaces[7/*#segments*/]; + struct interf* interf_props = NULL; + struct solid* solid_props = NULL; + struct ssp_rng* rng = NULL; + (void)argc, (void)argv; + + OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); + OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev)); + + fluid_shader.temperature = fluid_get_temperature; + OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + + /* Setup the solid shader */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + + /* Create the solid medium #1 */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = LAMBDA1; + solid_props->cp = 2; + solid_props->rho = 25; + solid_props->delta = DELTA1; + OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); + OK(sdis_data_ref_put(data)); + + /* Create the solid medium #2 */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = LAMBDA2; + solid_props->cp = 2; + solid_props->rho = 25; + solid_props->delta = DELTA2; + OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); + OK(sdis_data_ref_put(data)); + + /* Setup the interface shader */ + interf_shader.front.temperature = interface_get_temperature; + interf_shader.back.temperature = interface_get_temperature; + interf_shader.convection_coef = interface_get_convection_coef; + + /* Create the adiabatic interfaces */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + OK(sdis_interface_create + (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1)); + OK(sdis_interface_create + (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2)); + OK(sdis_data_ref_put(data)); + + /* Create the T0 interface */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = T0; + OK(sdis_interface_create + (dev, solid1, fluid, &interf_shader, data, &interf_T0)); + OK(sdis_data_ref_put(data)); + + /* Create the TL interface */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = TL; + OK(sdis_interface_create + (dev, solid2, fluid, &interf_shader, data, &interf_TL)); + OK(sdis_data_ref_put(data)); + + /* Create the solid1-solid2 interface */ + interf_shader.convection_coef = NULL; + interf_shader.thermal_contact_resistance = interface_get_contact_resistance; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + OK(sdis_interface_create + (dev, solid1, solid2, &interf_shader, data, &interf_R)); + OK(sdis_data_ref_put(data)); + + /* Release the media */ + OK(sdis_medium_ref_put(solid1)); + OK(sdis_medium_ref_put(solid2)); + OK(sdis_medium_ref_put(fluid)); + + /* Map the interfaces to their box triangles */ + /* Front */ + model3d_interfaces[0] = interf_adiabatic1; + model3d_interfaces[1] = interf_adiabatic1; + model3d_interfaces[2] = interf_adiabatic2; + model3d_interfaces[3] = interf_adiabatic2; + /* Left */ + model3d_interfaces[4] = interf_T0; + model3d_interfaces[5] = interf_T0; + /* Back */ + model3d_interfaces[6] = interf_adiabatic1; + model3d_interfaces[7] = interf_adiabatic1; + model3d_interfaces[8] = interf_adiabatic2; + model3d_interfaces[9] = interf_adiabatic2; + /* Right */ + model3d_interfaces[10] = interf_TL; + model3d_interfaces[11] = interf_TL; + /* Top */ + model3d_interfaces[12] = interf_adiabatic1; + model3d_interfaces[13] = interf_adiabatic1; + model3d_interfaces[14] = interf_adiabatic2; + model3d_interfaces[15] = interf_adiabatic2; + /* Bottom */ + model3d_interfaces[16] = interf_adiabatic1; + model3d_interfaces[17] = interf_adiabatic1; + model3d_interfaces[18] = interf_adiabatic2; + model3d_interfaces[19] = interf_adiabatic2; + /* Inside */ + model3d_interfaces[20] = interf_R; + model3d_interfaces[21] = interf_R; + + /* Map the interfaces to their square segments */ + /* Bottom */ + model2d_interfaces[0] = interf_adiabatic2; + model2d_interfaces[1] = interf_adiabatic1; + /* Left */ + model2d_interfaces[2] = interf_T0; + /* Top */ + model2d_interfaces[3] = interf_adiabatic1; + model2d_interfaces[4] = interf_adiabatic2; + /* Right */ + model2d_interfaces[5] = interf_TL; + /* Contact */ + model2d_interfaces[6] = interf_R; + + /* Create the box scene */ + scn_args.get_indices = model3d_get_indices; + scn_args.get_interface = model3d_get_interface; + scn_args.get_position = model3d_get_position; + scn_args.nprimitives = model3d_ntriangles; + scn_args.nvertices = model3d_nvertices; + scn_args.context = model3d_interfaces; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); + + /* Create the square scene */ + scn_args.get_indices = model2d_get_indices; + scn_args.get_interface = model2d_get_interface; + scn_args.get_position = model2d_get_position; + scn_args.nprimitives = model2d_nsegments; + scn_args.nvertices = model2d_nvertices; + scn_args.context = model2d_interfaces; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); + + /* Release the interfaces */ + OK(sdis_interface_ref_put(interf_adiabatic1)); + OK(sdis_interface_ref_put(interf_adiabatic2)); + OK(sdis_interface_ref_put(interf_T0)); + OK(sdis_interface_ref_put(interf_TL)); + OK(sdis_interface_ref_put(interf_R)); + + /* Solve */ + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); + printf(">> Box scene\n"); + solve_probe(box_scn, interf_props, rng); + solve(box_scn, interf_props, rng); + printf("\n>> Square scene\n"); + solve_probe(square_scn, interf_props, rng); + solve(square_scn, interf_props, rng); + + 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); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_data.c b/src/test_sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_device.c b/src/test_sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_unstationary_atm.c b/src/test_sdis_unstationary_atm.c @@ -0,0 +1,883 @@ +/* Copyright (C) 2016-2021 |Meso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sdis.h" +#include "test_sdis_utils.h" + +#include <rsys/clock_time.h> +#include <rsys/mem_allocator.h> +#include <rsys/double3.h> +#include <rsys/math.h> +#include <star/ssp.h> + +/* + * + * 3D 2D + * + * /////////////// /////////////// + * +-----------+-+ +-----------+-+ + * /' / /| | | | + * +-----------+-+ | HA _\ <---- TR | _\ | | HA _\ <---- TR + * | | _\ | | | / / <---- TR TG| HG / / HC | | / / <---- TR + * | |HG / / HC| | | TA \__/ <---- TR | \__/ | | TA \__/ <---- TR + * TG| | \__/ | | | | | | + * | +.........|.|.+ +-----------+-+ + * |, |/|/ /////H///////E/ + * +-----------+-+ + * //////H//////E/ + */ + +#define XH 3 +#define XHpE 3.2 +#define XE (XHpE - XH) + +#define T0_SOLID 300 +#define T0_FLUID 300 + +#define UNKNOWN_TEMPERATURE -1 + +#define N 10000 /* #realisations */ + +#define TG 310 +#define HG 400 + +#define HC 400 + +#define TA 290 +#define HA 400 +#define TR 260 + +/* hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon + * Tref = (hr / (4 * 5.6696e-8 * epsilon)) ^ 1/3, hr = 6 */ +#define TREF 297.974852286 + +#define RHO_F 1.3 +#define CP_F 1005 +#define RHO_S 2400 +#define CP_S 800 +#define LAMBDA 0.6 + +#define X_PROBE (XH + 0.2 * XE) + +#define DELTA (XE/30.0) + + /******************************************************************************* + * Box geometry + ******************************************************************************/ +static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/] += { + 0, 0, 0, + XH, 0, 0, + XHpE, 0, 0, + 0, XHpE, 0, + XH, XHpE, 0, + XHpE, XHpE, 0, + 0, 0, XHpE, + XH, 0, XHpE, + XHpE, 0, XHpE, + 0, XHpE, XHpE, + XH, XHpE, XHpE, + XHpE, XHpE, XHpE +}; +static const size_t model3d_nvertices = sizeof(model3d_vertices) / sizeof(double[3]); + +/* The following array lists the indices toward the 3D vertices of each + * triangle. + * ,3---,4---,5 ,3----4----5 ,4 + * ,' | ,' | ,'/| ,'/| \ | \ | ,'/| + * 9----10---11 / | 9' / | \ | \ | 10 / | Y + * |', |', | / ,2 | / ,0---,1---,2 | / ,1 | + * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X + * 6----7----8' 6----7'---8' 7 / + * Front, right Back, left and Internal Z + * and Top faces bottom faces face */ +static const size_t model3d_indices[22/*#triangles*/ * 3/*#indices per triangle*/] += { + 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */ + 0, 6, 3, 3, 6, 9, /* -X */ + 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */ + 5, 11, 8, 8, 2, 5, /* +X */ + 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */ + 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */ + 4, 10, 7, 7, 1, 4 /* Inside */ +}; +static const size_t model3d_ntriangles = sizeof(model3d_indices) / sizeof(size_t[3]); + +static INLINE void +model3d_get_indices(const size_t itri, size_t ids[3], void* context) +{ + (void)context; + CHK(ids); + CHK(itri < model3d_ntriangles); + ids[0] = model3d_indices[itri * 3 + 0]; + ids[1] = model3d_indices[itri * 3 + 1]; + ids[2] = model3d_indices[itri * 3 + 2]; +} + +static INLINE void +model3d_get_position(const size_t ivert, double pos[3], void* context) +{ + (void)context; + CHK(pos); + CHK(ivert < model3d_nvertices); + pos[0] = model3d_vertices[ivert * 3 + 0]; + pos[1] = model3d_vertices[ivert * 3 + 1]; + pos[2] = model3d_vertices[ivert * 3 + 2]; +} + +static INLINE void +model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context) +{ + struct sdis_interface** interfaces = context; + CHK(context && bound); + CHK(itri < model3d_ntriangles); + *bound = interfaces[itri]; +} + +/******************************************************************************* + * Square geometry + ******************************************************************************/ +static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = { + XHpE, 0, + XH, 0, + 0, 0, + 0, XHpE, + XH, XHpE, + XHpE, XHpE +}; +static const size_t model2d_nvertices = sizeof(model2d_vertices) / sizeof(double[2]); + +static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = { + 0, 1, 1, 2, /* Bottom */ + 2, 3, /* Left */ + 3, 4, 4, 5, /* Top */ + 5, 0, /* Right */ + 4, 1 /* Inside */ +}; +static const size_t model2d_nsegments = sizeof(model2d_indices) / sizeof(size_t[2]); + + +static INLINE void +model2d_get_indices(const size_t iseg, size_t ids[2], void* context) +{ + (void)context; + CHK(ids); + CHK(iseg < model2d_nsegments); + ids[0] = model2d_indices[iseg * 2 + 0]; + ids[1] = model2d_indices[iseg * 2 + 1]; +} + +static INLINE void +model2d_get_position(const size_t ivert, double pos[2], void* context) +{ + (void)context; + CHK(pos); + CHK(ivert < model2d_nvertices); + pos[0] = model2d_vertices[ivert * 2 + 0]; + pos[1] = model2d_vertices[ivert * 2 + 1]; +} + +static INLINE void +model2d_get_interface +(const size_t iseg, struct sdis_interface** bound, void* context) +{ + struct sdis_interface** interfaces = context; + CHK(context && bound); + CHK(iseg < model2d_nsegments); + *bound = interfaces[iseg]; +} + +/******************************************************************************* + * Media + ******************************************************************************/ +struct solid { + double lambda; + double rho; + double cp; + double delta; + double temperature; + double t0; +}; + +static double +solid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct solid* solid; + CHK(vtx && data); + solid = ((struct solid*)sdis_data_cget(data)); + return solid->cp; +} + +static double +solid_get_thermal_conductivity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct solid* solid; + CHK(vtx && data); + solid = ((struct solid*)sdis_data_cget(data)); + return solid->lambda; +} + +static double +solid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct solid* solid; + CHK(vtx && data); + solid = ((struct solid*)sdis_data_cget(data)); + return solid->rho; +} + +static double +solid_get_delta + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct solid* solid; + CHK(vtx && data); + solid = ((struct solid*)sdis_data_cget(data)); + return solid->delta; +} + +static double +solid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct solid* solid; + CHK(vtx && data); + solid = ((struct solid*)sdis_data_cget(data)); + if(vtx->time <= solid->t0) + return solid->temperature; + return UNKNOWN_TEMPERATURE; +} + +struct fluid { + double rho; + double cp; + double t0; + double temperature; +}; + +static double +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct fluid* fluid; + CHK(vtx && data); + fluid = ((struct fluid*)sdis_data_cget(data)); + if(vtx->time <= fluid->t0) + return fluid->temperature; + return UNKNOWN_TEMPERATURE; +} + +static double +fluid_get_volumic_mass + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct fluid* fluid; + CHK(vtx && data); + fluid = ((struct fluid*)sdis_data_cget(data)); + return fluid->rho; +} + +static double +fluid_get_calorific_capacity + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + struct fluid* fluid; + CHK(vtx && data); + fluid = ((struct fluid*)sdis_data_cget(data)); + return fluid->cp; +} + +/******************************************************************************* + * Interfaces + ******************************************************************************/ +struct interf { + double temperature; + double emissivity; + double h; +}; + +static double +interface_get_temperature + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf; + CHK(frag && data); + interf = sdis_data_cget(data); + return interf->temperature; +} + +static double +interface_get_convection_coef + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf; + CHK(frag && data); + interf = sdis_data_cget(data); + return interf->h; +} + +static double +interface_get_emissivity + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf; + CHK(frag && data); + interf = sdis_data_cget(data); + return interf->emissivity; +} + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +solve_tbound1 + (struct sdis_scene* scn, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_boundary_args solve_args + = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + enum sdis_scene_dimension dim; + double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; + double ref[sizeof(t) / sizeof(*t)] + = { 290.046375, 289.903935, 289.840490, 289.802690, 289.777215, + 289.759034, 289.745710, 289.735826, 289.728448, 289.722921 }; + const int nsimuls = sizeof(t) / sizeof(*t); + int isimul; + ASSERT(scn && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + solve_args.side = SDIS_FRONT; + FOR_EACH(isimul, 0, nsimuls) { + solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; + if(dim == SDIS_SCENE_2D) { + solve_args.iprim = model2d_nsegments - 1; + solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1); + } else { + double u = ssp_rng_uniform_double(rng, 0, 1); + double v = ssp_rng_uniform_double(rng, 0, 1); + double w = ssp_rng_uniform_double(rng, 0, 1); + double x = 1 / (u + v + w); + solve_args.iprim = (isimul % 2) ? 10 : 11; /* +X face */ + solve_args.uv[0] = u * x; + solve_args.uv[1] = v * x; + } + + time_current(&t0); + OK(sdis_solve_probe_boundary(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("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul], + ref[isimul], T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul], + ref[isimul], 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(t[isimul] == 0 || eq_eps(T.E, ref[isimul], T.SE * 4)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +static void +solve_tbound2 + (struct sdis_scene* scn, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_boundary_args solve_args + = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + enum sdis_scene_dimension dim; + double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; + double ref[sizeof(t) / sizeof(*t)] + = { 309.08032, 309.34626, 309.46525, 309.53625, 309.58408, + 309.618121, 309.642928, 309.661167, 309.674614, 309.684524 }; + const int nsimuls = sizeof(t) / sizeof(*t); + int isimul; + ASSERT(scn && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + solve_args.side = SDIS_FRONT; + FOR_EACH(isimul, 0, nsimuls) { + solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; + if(dim == SDIS_SCENE_2D) { + solve_args.iprim = model2d_nsegments - 1; + solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1); + } else { + double u = ssp_rng_uniform_double(rng, 0, 1); + double v = ssp_rng_uniform_double(rng, 0, 1); + double w = ssp_rng_uniform_double(rng, 0, 1); + double x = 1 / (u + v + w); + solve_args.iprim = (isimul % 2) ? 20 : 21; /* Internal face */ + solve_args.uv[0] = u * x; + solve_args.uv[1] = v * x; + } + + time_current(&t0); + OK(sdis_solve_probe_boundary(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("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul], + ref[isimul], T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n", + (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul], + ref[isimul], 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[isimul], T.SE * 4)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +static void +solve_tsolid + (struct sdis_scene* scn, + struct ssp_rng* rng) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + enum sdis_scene_dimension dim; + double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; + double ref[sizeof(t) / sizeof(*t)] + = { 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030, + 304.75041, 305.01595, 305.21193, 305.35641, 305.46271 }; + const int nsimuls = sizeof(t) / sizeof(*t); + int isimul; + ASSERT(scn && rng); + + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + FOR_EACH(isimul, 0, nsimuls) { + solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; + solve_args.position[0] = X_PROBE; + solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE); + + if(dim == SDIS_SCENE_3D) + solve_args.position[2] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE); + + 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("Unstationary temperature at (%g,%g) at t=%g = %g ~ %g +/- %g\n", + SPLIT2(solve_args.position), t[isimul], ref[isimul], T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Unstationary temperature at (%g,%g,%g) at t=%g = %g ~ %g +/- %g\n", + SPLIT3(solve_args.position), t[isimul], ref[isimul], 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[isimul], T.SE * 4)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +static void +solve_tfluid + (struct sdis_scene* scn) +{ + char dump[128]; + struct time t0, t1; + struct sdis_estimator* estimator; + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_mc T = SDIS_MC_NULL; + struct sdis_mc time = SDIS_MC_NULL; + size_t nreals; + size_t nfails; + enum sdis_scene_dimension dim; + double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; + double ref[sizeof(t) / sizeof(*t)] + = { 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899, + 309.82141, 309.83055, 309.83728, 309.84224 }; + const int nsimuls = sizeof(t) / sizeof(*t); + int isimul; + ASSERT(scn); + + OK(sdis_scene_get_dimension(scn, &dim)); + + solve_args.nrealisations = N; + solve_args.position[0] = XH * 0.5; + solve_args.position[1] = XH * 0.5; + solve_args.position[2] = XH * 0.5; + FOR_EACH(isimul, 0, nsimuls) { + solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; + + 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("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n", + t[isimul], ref[isimul], T.E, T.SE); + break; + case SDIS_SCENE_3D: + printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n", + t[isimul], ref[isimul], 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[isimul], T.SE * 4)); + + OK(sdis_estimator_ref_put(estimator)); + } +} + +/******************************************************************************* + * 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* fluid_A = NULL; + struct sdis_medium* solid = NULL; + struct sdis_medium* dummy_solid = NULL; + struct sdis_interface* interf_adiabatic_1 = NULL; + struct sdis_interface* interf_adiabatic_2 = NULL; + struct sdis_interface* interf_TG = NULL; + struct sdis_interface* interf_P = NULL; + struct sdis_interface* interf_TA = 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* model3d_interfaces[22 /*#triangles*/]; + struct sdis_interface* model2d_interfaces[7/*#segments*/]; + struct interf* interf_props = NULL; + struct solid* solid_props = NULL; + struct fluid* fluid_props = NULL; + struct ssp_rng* rng = NULL; + (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 solid shader */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + + /* Create the solid media */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = LAMBDA; + solid_props->cp = CP_S; + solid_props->rho = RHO_S; + solid_props->delta = DELTA; + solid_props->t0 = 0; + solid_props->temperature = T0_SOLID; + OK(sdis_solid_create(dev, &solid_shader, data, &solid)); + OK(sdis_data_ref_put(data)); + + /* Create a dummy solid media to be used outside the model */ + OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); + solid_props = sdis_data_get(data); + solid_props->lambda = 0; + solid_props->cp = 1; + solid_props->rho = 1; + solid_props->delta = 1; + solid_props->t0 = INF; + solid_props->temperature = UNKNOWN_TEMPERATURE; + OK(sdis_solid_create(dev, &solid_shader, data, &dummy_solid)); + OK(sdis_data_ref_put(data)); + + /* Setup the fluid shader */ + fluid_shader.calorific_capacity = fluid_get_calorific_capacity; + fluid_shader.volumic_mass = fluid_get_volumic_mass; + fluid_shader.temperature = fluid_get_temperature; + + /* Create the internal fluid media */ + OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data)); + fluid_props = sdis_data_get(data); + fluid_props->cp = CP_F; + fluid_props->rho = RHO_F; + fluid_props->t0 = 0; + fluid_props->temperature = T0_FLUID; + OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); + OK(sdis_data_ref_put(data)); + + /* Create the 'A' fluid media */ + OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data)); + fluid_props = sdis_data_get(data); + fluid_props->cp = 1; + fluid_props->rho = 1; + fluid_props->t0 = INF; + fluid_props->temperature = TA; + OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid_A)); + OK(sdis_data_ref_put(data)); + + /* Setup the interface shader */ + interf_shader.front.temperature = interface_get_temperature; + interf_shader.back.temperature = interface_get_temperature; + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.convection_coef_upper_bound = 0; + + /* Create the adiabatic interfaces */ + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + interf_props->h = 0; + OK(sdis_interface_create + (dev, fluid, dummy_solid, &interf_shader, data, &interf_adiabatic_1)); + OK(sdis_data_ref_put(data)); + + interf_shader.convection_coef = NULL; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + interf_props->h = 0; + OK(sdis_interface_create + (dev, solid, dummy_solid, &interf_shader, data, &interf_adiabatic_2)); + OK(sdis_data_ref_put(data)); + + /* Create the P interface */ + interf_shader.convection_coef_upper_bound = HC; + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.front.emissivity = interface_get_emissivity; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + interf_props->h = HC; + interf_props->emissivity = 1; + OK(sdis_interface_create + (dev, fluid, solid, &interf_shader, data, &interf_P)); + OK(sdis_data_ref_put(data)); + + /* Create the TG interface */ + interf_shader.convection_coef_upper_bound = HG; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = TG; + interf_props->h = HG; + interf_props->emissivity = 1; + OK(sdis_interface_create + (dev, fluid, dummy_solid, &interf_shader, data, &interf_TG)); + OK(sdis_data_ref_put(data)); + + /* Create the TA interface */ + interf_shader.convection_coef_upper_bound = HA; + interf_shader.front.emissivity = NULL; + interf_shader.back.emissivity = interface_get_emissivity; + OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); + interf_props = sdis_data_get(data); + interf_props->temperature = UNKNOWN_TEMPERATURE; + interf_props->h = HA; + interf_props->emissivity = 1; + OK(sdis_interface_create + (dev, solid, fluid_A, &interf_shader, data, &interf_TA)); + OK(sdis_data_ref_put(data)); + + /* Release the media */ + OK(sdis_medium_ref_put(solid)); + OK(sdis_medium_ref_put(dummy_solid)); + OK(sdis_medium_ref_put(fluid)); + OK(sdis_medium_ref_put(fluid_A)); + + /* Front */ + model3d_interfaces[0] = interf_adiabatic_1; + model3d_interfaces[1] = interf_adiabatic_1; + model3d_interfaces[2] = interf_adiabatic_2; + model3d_interfaces[3] = interf_adiabatic_2; + /* Left */ + model3d_interfaces[4] = interf_TG; + model3d_interfaces[5] = interf_TG; + /* Back */ + model3d_interfaces[6] = interf_adiabatic_1; + model3d_interfaces[7] = interf_adiabatic_1; + model3d_interfaces[8] = interf_adiabatic_2; + model3d_interfaces[9] = interf_adiabatic_2; + /* Right */ + model3d_interfaces[10] = interf_TA; + model3d_interfaces[11] = interf_TA; + /* Top */ + model3d_interfaces[12] = interf_adiabatic_1; + model3d_interfaces[13] = interf_adiabatic_1; + model3d_interfaces[14] = interf_adiabatic_2; + model3d_interfaces[15] = interf_adiabatic_2; + /* Bottom */ + model3d_interfaces[16] = interf_adiabatic_1; + model3d_interfaces[17] = interf_adiabatic_1; + model3d_interfaces[18] = interf_adiabatic_2; + model3d_interfaces[19] = interf_adiabatic_2; + /* Inside */ + model3d_interfaces[20] = interf_P; + model3d_interfaces[21] = interf_P; + + /* Bottom */ + model2d_interfaces[0] = interf_adiabatic_2; + model2d_interfaces[1] = interf_adiabatic_1; + /* Left */ + model2d_interfaces[2] = interf_TG; + /* Top */ + model2d_interfaces[3] = interf_adiabatic_1; + model2d_interfaces[4] = interf_adiabatic_2; + /* Right */ + model2d_interfaces[5] = interf_TA; + /* Contact */ + model2d_interfaces[6] = interf_P; + + /* Create the box scene */ + scn_args.get_indices = model3d_get_indices; + scn_args.get_interface = model3d_get_interface; + scn_args.get_position = model3d_get_position; + scn_args.nprimitives = model3d_ntriangles; + scn_args.nvertices = model3d_nvertices; + scn_args.context = model3d_interfaces; + scn_args.trad = TR; + scn_args.tref = TREF; + OK(sdis_scene_create(dev, &scn_args, &box_scn)); + + /* Create the square scene */ + scn_args.get_indices = model2d_get_indices; + scn_args.get_interface = model2d_get_interface; + scn_args.get_position = model2d_get_position; + scn_args.nprimitives = model2d_nsegments; + scn_args.nvertices = model2d_nvertices; + scn_args.context = model2d_interfaces; + scn_args.trad = TR; + scn_args.tref = TREF; + OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); + + /* Release the interfaces */ + OK(sdis_interface_ref_put(interf_adiabatic_1)); + OK(sdis_interface_ref_put(interf_adiabatic_2)); + OK(sdis_interface_ref_put(interf_TG)); + OK(sdis_interface_ref_put(interf_P)); + OK(sdis_interface_ref_put(interf_TA)); + + /* Solve */ + OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng)); + printf(">> Box scene\n"); + solve_tfluid(box_scn); + solve_tbound1(box_scn, rng); + solve_tbound2(box_scn, rng); + solve_tsolid(box_scn, rng); + printf("\n>> Square scene\n"); + solve_tfluid(square_scn); + solve_tbound1(box_scn, rng); + solve_tbound2(box_scn, rng); + solve_tsolid(square_scn, rng); + + 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); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -191,6 +191,7 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = { static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = { dummy_interface_getter, 0, + dummy_interface_getter, DUMMY_INTERFACE_SIDE_SHADER__, DUMMY_INTERFACE_SIDE_SHADER__ }; diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -233,7 +233,6 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) struct sdis_mc T = SDIS_MC_NULL; size_t nreals; size_t nfails; - double pos[3] = {0,0}; size_t i; solve_args.time_range[0] = INF; @@ -252,7 +251,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) OK(sdis_estimator_get_failure_count(estimator, &nfails)); Tc = T.E - 273.15; /* Convert in Celcius */ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n", - SPLIT3(pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); + SPLIT3(refs[i].pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); printf("#realisations: %lu; #failures: %lu\n", (unsigned long)nreals, (unsigned long)nfails); /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/ diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 @@ -253,7 +253,6 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; size_t nreals; size_t nfails; - double pos[2] = {0,0}; size_t i; solve_args.nrealisations = N; @@ -271,7 +270,7 @@ check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) OK(sdis_estimator_get_failure_count(estimator, &nfails)); Tc = T.E - 273.15; /* Convert in Celcius */ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", - SPLIT2(pos), refs[i].temperature, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); + SPLIT2(refs[i].pos), refs[i].temperature, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); printf("#realisations: %lu; #failures: %lu\n", (unsigned long)nreals, (unsigned long)nfails); /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/ diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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 diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2021 |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