commit 82746bff05b3bc03507591dd1b18b630394dcfc9
parent aee9f6c7e4b84062e127cf93125839450eb59cb5
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 13 Apr 2018 11:35:53 +0200
Add a simple transient convection test
Diffstat:
2 files changed, 302 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -125,6 +125,7 @@ if(NOT NO_TEST)
new_test(test_sdis_camera)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
+ new_test(test_sdis_convection)
new_test(test_sdis_data)
new_test(test_sdis_device)
new_test(test_sdis_interface)
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -0,0 +1,301 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/double3.h>
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of an unit fluid cube/square whose temperature is
+ * unknown. The convection coefficient with the surrounding solid is H
+ * everywhere the temperature of the -/+X, -/+Y and -/+Z faces are fixed to T0
+ * and T1, T2, T3, T4 and T5, respectively. This test computes the temperature
+ * of the fluid Tf at an observation time t. This temperature is equal to:
+ *
+ * Tf(t) = Tf(0) * e^(-nu*t) + Tinf*(1-e^(-nu*t))
+ *
+ * nu = (Sum_{i=0..5}(H*Si)) / (RHO*CP*V)
+ * Tinf = (Sum_{i=0..5}(H*Si*Ti) / (Sum_{i=0..5}(H*Si));
+ *
+ * with Si surface of the faces (i.e. one), V the volume of the cube (i.e.
+ * one), RHO the volumic mass of the fluid and CP its calorific capacity.
+ *
+ * 3D 2D
+ *
+ * (1,1,1) (1,1)
+ * +---------+ +-----T3----+
+ * /' T3 /|T4 | |
+ * +---------+ | | H _\ |
+ * | ' H _\ |T1 T0 / / T1
+ * |T0 / / | | | \__/ |
+ * | +..\__/.|.+ | |
+ * T5|, T2 |/ +-----------+
+ * +---------+ (0,0)
+ * (0,0,0)
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 100000 /* #realisations */
+
+#define Tf_0 280.0
+
+#define T0 300.0
+#define T1 310.0
+#define T2 320.0
+#define T3 330.0
+#define T4 340.0
+#define T5 350.0
+
+#define H 10.0
+#define RHO 25.0
+#define CP 2.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return RHO;
+}
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return CP;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+};
+
+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 H;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct sdis_interface_shader* interf_shader,
+ const double temperature)
+{
+ struct sdis_data* data;
+ struct sdis_interface* interf;
+ struct interf* interf_props;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = temperature;
+ CHK(sdis_interface_create
+ (dev, front, back, interf_shader, data, &interf) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+ return interf;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_T1 = NULL;
+ struct sdis_interface* interf_T2 = NULL;
+ struct sdis_interface* interf_T3 = NULL;
+ struct sdis_interface* interf_T4 = NULL;
+ struct sdis_interface* interf_T5 = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface* box_interfaces[12/*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ double pos[3];
+ double ref;
+ double Tinf;
+ double nu;
+ double time;
+ size_t nreals;
+ size_t nfails;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.temperature = interface_get_temperature;
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.emissivity = interface_get_emissivity;
+ interf_shader.specular_fraction = interface_get_specular_fraction;
+
+ /* Create the interfaces */
+ interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0);
+ interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1);
+ interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2);
+ interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3);
+ interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4);
+ interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_T5; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_T2; /* Bottom */
+ square_interfaces[1] = interf_T0; /* Left */
+ square_interfaces[2] = interf_T3; /* Top */
+ square_interfaces[3] = interf_T1; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_T0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T3) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T4) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T5) == RES_OK);
+
+ d3_splat(pos, 0.25);
+
+ nu = (6*H)/(RHO*CP);
+ time = 1.0/nu;
+ Tinf = (H*(T0+T1+T2+T3+T4+T5))/(6*H);
+ ref = Tf_0 * exp(-nu*time) + Tinf*(1-exp(-nu*time));
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE*2));
+
+ nu = (4*H)/(RHO*CP);
+ time = 1.0/nu;
+ Tinf = (H*(T0+T1+T2+T3))/(4*H);
+ ref = Tf_0 * exp(-nu*time) + Tinf*(1-exp(-nu*time));
+
+ /* Solve in 2D */
+ CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE*2.0));
+
+ CHK(sdis_scene_ref_put(box_scn) == RES_OK);
+ CHK(sdis_scene_ref_put(square_scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
