commit 2a15100483e3a6cfcaf687d2c8a8d7e50526f0a0
parent 25bee04e88d63598844f9736e4bc173cbb9edb49
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 17 Jul 2020 15:38:34 +0200
Test the sdis_compute_mean_power function
Diffstat:
2 files changed, 345 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -172,6 +172,7 @@ if(NOT NO_TEST)
endfunction()
new_test(test_sdis_camera)
+ new_test(test_sdis_compute_mean_power)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
new_test(test_sdis_convection)
@@ -208,6 +209,7 @@ if(NOT NO_TEST)
add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
endif()
+ target_link_libraries(test_sdis_compute_mean_power Star3DUT)
target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT)
target_link_libraries(test_sdis_solve_medium Star3DUT)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
diff --git a/src/test_sdis_compute_mean_power.c b/src/test_sdis_compute_mean_power.c
@@ -0,0 +1,343 @@
+/* Copyright (C) 2016-2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/stretchy_array.h>
+#include <star/s3dut.h>
+
+#define UNKOWN_TEMPERATURE -1
+#define N 100000ul /* #realisations */
+#define POWER0 10
+#define POWER1 5
+
+static INLINE void
+check_intersection
+ (const double val0,
+ const double eps0,
+ const double val1,
+ const double eps1)
+{
+ double interval0[2], interval1[2];
+ double intersection[2];
+ interval0[0] = val0 - eps0;
+ interval0[1] = val0 + eps0;
+ interval1[0] = val1 - eps1;
+ interval1[1] = val1 + eps1;
+ intersection[0] = MMAX(interval0[0], interval1[0]);
+ intersection[1] = MMIN(interval0[1], interval1[1]);
+ CHK(intersection[0] <= intersection[1]);
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ struct s3dut_mesh_data msh0;
+ struct s3dut_mesh_data msh1;
+ struct sdis_interface* interf0;
+ struct sdis_interface* interf1;
+};
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ const struct context* ctx = context;
+ /* Note that we swap the indices to ensure that the triangle normals point
+ * inward the super shape */
+ if(itri < ctx->msh0.nprimitives) {
+ ids[0] = ctx->msh0.indices[itri*3+0];
+ ids[2] = ctx->msh0.indices[itri*3+1];
+ ids[1] = ctx->msh0.indices[itri*3+2];
+ } else {
+ const size_t itri2 = itri - ctx->msh0.nprimitives;
+ ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices;
+ ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices;
+ ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices;
+ }
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ const struct context* ctx = context;
+ if(ivert < ctx->msh0.nvertices) {
+ pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0;
+ pos[1] = ctx->msh0.positions[ivert*3+1];
+ pos[2] = ctx->msh0.positions[ivert*3+2];
+ } else {
+ const size_t ivert2 = ivert - ctx->msh0.nvertices;
+ pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0;
+ pos[1] = ctx->msh1.positions[ivert2*3+1];
+ pos[2] = ctx->msh1.positions[ivert2*3+2];
+ }
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL); (void)data;
+ return 1.0;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx == NULL); (void)data;
+ return 300;
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0 / 20.0;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return vtx->P[0] < 0 ? POWER0 : POWER1;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct context ctx;
+ struct s3dut_mesh* sphere = NULL;
+ struct s3dut_mesh* cylinder = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid0 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_interface* interf0 = NULL;
+ struct sdis_interface* interf1 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_mc mpow = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_compute_mean_power_args args = SDIS_COMPUTE_MEAN_POWER_ARGS_DEFAULT;
+ size_t nverts = 0;
+ size_t ntris = 0;
+ double ref = 0;
+ (void)argc, (void) argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the fluid */
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* Create the solids */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+
+ /* Create the interface0 */
+ OK(sdis_interface_create(dev, solid0, fluid, &interf_shader, NULL, &interf0));
+ OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, NULL, &interf1));
+ ctx.interf0 = interf0;
+ ctx.interf1 = interf1;
+
+ /* Create the geometry */
+ OK(s3dut_create_sphere(&allocator, 1, 512, 256, &sphere));
+ OK(s3dut_create_cylinder(&allocator, 1, 10, 512, 8, &cylinder));
+ OK(s3dut_mesh_get_data(sphere, &ctx.msh0));
+ OK(s3dut_mesh_get_data(cylinder, &ctx.msh1));
+
+ /* Create the scene */
+ ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives;
+ nverts = ctx.msh0.nvertices + ctx.msh1.nvertices;
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ /* Test sdis_compute_mean_power function */
+ args.nrealisations = N;
+ args.medium = solid0;
+ args.time_range[0] = INF;
+ args.time_range[1] = INF;
+ BA(sdis_compute_mean_power(NULL, &args, &estimator));
+ BA(sdis_compute_mean_power(scn, NULL, &estimator));
+ BA(sdis_compute_mean_power(scn, &args, NULL));
+ args.nrealisations = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.nrealisations = N;
+ args.medium = NULL;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.medium = solid0;
+ args.fp_to_meter = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.fp_to_meter = 1;
+ args.time_range[0] = args.time_range[1] = -1;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[0] = 1;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[1] = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[0] = args.time_range[1] = INF;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+
+ BA(sdis_estimator_get_mean_power(NULL, &mpow));
+ BA(sdis_estimator_get_mean_power(estimator, NULL));
+ OK(sdis_estimator_get_mean_power(estimator, &mpow));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ /* Check results for solid 0 */
+ ref = 4.0/3.0 * PI * POWER0;
+ printf("Mean power of the solid0 = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Check results for solid 1 */
+ args.medium = solid1;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_power(estimator, &mpow));
+ ref = PI * 10 * POWER1;
+ printf("Mean power of the solid1 = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Check for a not null time range */
+ args.time_range[0] = 0;
+ args.time_range[1] = 10;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_power(estimator, &mpow));
+ ref = PI * 10 * POWER1 / 10;
+ printf("Mean power of the solid1 in [0, 10] s = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Reset the scene with only one solid medium */
+ OK(sdis_scene_ref_put(scn));
+ ctx.interf0 = interf0;
+ ctx.interf1 = interf0;
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ /* Check invalid medium */
+ args.time_range[0] = args.time_range[1] = 1;
+ args.medium = solid1;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+
+ /* Check non constant volumic power */
+ args.medium = solid0;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_power(estimator, &mpow));
+ ref = 4.0/3.0*PI*POWER0 + PI*10*POWER1;
+ printf("Mean power of the sphere+cylinder = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-1, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+#if 0
+ {
+ double* vertices = NULL;
+ size_t* indices = NULL;
+ size_t i;
+ CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices)));
+ CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices)));
+ FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx);
+ FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx);
+ dump_mesh(stdout, vertices, nverts, indices, ntris);
+ MEM_RM(&allocator, vertices);
+ MEM_RM(&allocator, indices);
+ }
+#endif
+
+ /* Clean up memory */
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid0));
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_interface_ref_put(interf0));
+ OK(sdis_interface_ref_put(interf1));
+ OK(sdis_scene_ref_put(scn));
+ OK(s3dut_mesh_ref_put(sphere));
+ OK(s3dut_mesh_ref_put(cylinder));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
