commit 2182cfceaff407632f10e6dfe171a375f4e90ea4
parent 7201428a5b2f4ec5a81019ce6fe33a210799bec8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 19 Feb 2018 16:38:45 +0100
Add a conducto radiative test in 2D
Diffstat:
5 files changed, 397 insertions(+), 7 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
+compile_commands.json
.gitignore
CMakeCache.txt
CMakeFiles
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -119,10 +119,11 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe)
new_test(test_sdis_solve_probe2)
new_test(test_sdis_solve_probe3)
+ new_test(test_sdis_conducto_radiative)
new_test(test_sdis_solve_probe_2d)
new_test(test_sdis_solve_probe2_2d)
new_test(test_sdis_solve_probe3_2d)
- new_test(test_sdis_conducto_radiative)
+ new_test(test_sdis_conducto_radiative_2d)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
if(CMAKE_COMPILER_IS_GNUCC)
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -312,7 +312,7 @@ main(int argc, char** argv)
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
- /* Create the interface that force to keep in conduction */
+ /* Create the interface that forces to keep in conduction */
interf.temperature = UNKNOWN_TEMPERATURE;
interf.convection_coef = -1;
interf.emissivity = -1;
@@ -333,14 +333,14 @@ main(int argc, char** argv)
interf.specular_fraction = 1;
create_interface(dev, fluid, solid2, &interf, interfaces+2);
- /* Create a the interface with a limit condition of 300K */
+ /* Create the interface with a limit condition of T0 Kelvin */
interf.temperature = T0;
interf.convection_coef = 0;
interf.emissivity = 1;
interf.specular_fraction = 1;
create_interface(dev, fluid, solid2, &interf, interfaces+3);
- /* Create a the interface with a limit condition of 310K */
+ /* Create the interface with a limit condition of T1 Kelvin */
interf.temperature = T1;
interf.convection_coef = 0;
interf.emissivity = 1;
@@ -355,14 +355,14 @@ main(int argc, char** argv)
prim_interfaces[8] = prim_interfaces[9] = interfaces[0];
prim_interfaces[10] = prim_interfaces[11] = interfaces[0];
- /* Setup the per primitive interfaces of the fluid on the left of the medium */
+ /* Setup the per primitive interface of the fluid on the left of the medium */
prim_interfaces[12] = prim_interfaces[13] = interfaces[2];
prim_interfaces[14] = prim_interfaces[15] = interfaces[3];
prim_interfaces[16] = prim_interfaces[17] = interfaces[2];
prim_interfaces[18] = prim_interfaces[19] = interfaces[2];
prim_interfaces[20] = prim_interfaces[21] = interfaces[2];
- /* Setup the per primitive interfaces of the fluid on the right of the medium */
+ /* Setup the per primitive interface of the fluid on the right of the medium */
prim_interfaces[22] = prim_interfaces[23] = interfaces[2];
prim_interfaces[24] = prim_interfaces[25] = interfaces[2];
prim_interfaces[26] = prim_interfaces[27] = interfaces[4];
@@ -382,7 +382,7 @@ main(int argc, char** argv)
Ts1 = T1 - tmp;
/* Run the simulations */
- ssp_rng_create(&allocator, &ssp_rng_kiss, &rng);
+ CHK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng) == RES_OK);
FOR_EACH(isimul, 0, nsimuls) {
struct sdis_mc T = SDIS_MC_NULL;
struct sdis_estimator* estimator;
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -0,0 +1,386 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (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 <star/ssp.h>
+
+#define UNKNOWN_TEMPERATURE -1
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct geometry {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface** interfaces;
+};
+
+static const double vertices[8/*#vertices*/*2/*#coords par vertex*/] = {
+ 1.0, -1.0,
+ -1.0, -1.0,
+ -1.0, 1.0,
+ 1.0, 1.0,
+ 1.5, -1.0,
+ -1.5, -1.0,
+ -1.5, 1.0,
+ 1.5, 1.0
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(double[2]);
+
+static const size_t indices[10/*#segments*/*2/*#indices per segment*/] = {
+ 0, 1, /* Solid bottom segment */
+ 1, 2, /* Solid left segment */
+ 2, 3, /* Solid top segment */
+ 3, 0, /* Solid right segment */
+
+ 1, 5, /* Left fluid bottom segment */
+ 5, 6, /* Left fluid left segment */
+ 6, 2, /* Left fluid top segment */
+
+ 4, 0, /* Right fluid bottom segment */
+ 3, 7, /* Right fluid top segment */
+ 7, 4 /* Right fluid right segment */
+};
+static const size_t nsegments = sizeof(indices) / sizeof(size_t[2]);
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ ids[0] = geom->indices[iseg*2+0];
+ ids[1] = geom->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ pos[0] = geom->positions[ivert*2+0];
+ pos[1] = geom->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ *bound = geom->interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+};
+
+static double
+temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return -1;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL);
+ CHK(data != NULL);
+ return ((const 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 != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0/10.0;
+}
+
+static double
+solid_get_delta_boundary
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 2.1/10.0;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interface {
+ double temperature;
+ double convection_coef;
+ double emissivity;
+ double specular_fraction;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interface*)sdis_data_cget(data))->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interface*)sdis_data_cget(data))->convection_coef;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interface*)sdis_data_cget(data))->emissivity;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interface*)sdis_data_cget(data))->specular_fraction;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct interface* interf,
+ struct sdis_interface** out_interf)
+{
+ struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_data* data = NULL;
+
+ CHK(interf != NULL);
+
+ shader.temperature = interface_get_temperature;
+ shader.convection_coef = interface_get_convection_coef;
+ shader.emissivity = interface_get_emissivity;
+ shader.specular_fraction = interface_get_specular_fraction;
+
+ CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface),
+ NULL, &data) == RES_OK);
+ *((struct interface*)sdis_data_get(data)) = *interf;
+
+ CHK(sdis_interface_create(dev, front, back, &shader, data, out_interf) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+}
+
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct interface interf;
+ struct geometry geom;
+ struct ssp_rng* rng = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_interface* interfaces[5] = {NULL};
+ struct sdis_interface* prim_interfaces[10/*#segment*/];
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ const size_t nsimuls = 4;
+ size_t isimul;
+ const double emissivity = 1;/* Emissivity of the side +/-X of the solid */
+ const double lambda = 0.1; /* Conductivity of the solid */
+ const double Tref = 300; /* Reference temperature */
+ const double T0 = 300; /* Fixed temperature on the left side of the system */
+ const double T1 = 310; /* Fixed temperature on the right side of the system */
+ const double thickness = 2.0; /* Thickness of the solid along X */
+ double Ts0, Ts1, hr, tmp;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = temperature_unknown;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = lambda;
+ 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.delta_boundary = solid_get_delta_boundary;
+ solid_shader.temperature = temperature_unknown;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the surrounding solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
+ NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = 0;
+ solid_shader.calorific_capacity = solid_get_thermal_conductivity;
+ 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.delta_boundary = solid_get_delta_boundary;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the interface that forces to keep in conduction */
+ interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.convection_coef = -1;
+ interf.emissivity = -1;
+ interf.specular_fraction = -1;
+ create_interface(dev, solid, solid2, &interf, interfaces+0);
+
+ /* Create the interface that emits radiative heat from the solid */
+ interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.convection_coef = 0;
+ interf.emissivity = emissivity;
+ interf.specular_fraction = -1;
+ create_interface(dev, solid, fluid, &interf, interfaces+1);
+
+ /* Create the interface that forces the radiative heat to bounce */
+ interf.temperature = UNKNOWN_TEMPERATURE;
+ interf.convection_coef = 0;
+ interf.emissivity = 0;
+ interf.specular_fraction = 1;
+ create_interface(dev, fluid, solid2, &interf, interfaces+2);
+
+ /* Create the interface with a limit condition of T0 Kelvin */
+ interf.temperature = T0;
+ interf.convection_coef = 0;
+ interf.emissivity = 1;
+ interf.specular_fraction = 1;
+ create_interface(dev, fluid, solid2, &interf, interfaces+3);
+
+ /* Create the interface with a limit condition of T1 Kelvin */
+ interf.temperature = T1;
+ interf.convection_coef = 0;
+ interf.emissivity = 1;
+ interf.specular_fraction = 1;
+ create_interface(dev, fluid, solid2, &interf, interfaces+4);
+
+ /* Setup the per primitive interface of the solid medium */
+ prim_interfaces[0] = interfaces[0];
+ prim_interfaces[1] = interfaces[1];
+ prim_interfaces[2] = interfaces[0];
+ prim_interfaces[3] = interfaces[1];
+
+ /* Setup the per primitive interface of the fluid on the left of the medium */
+ prim_interfaces[4] = interfaces[2];
+ prim_interfaces[5] = interfaces[3];
+ prim_interfaces[6] = interfaces[2];
+
+ /* Setup the per primitive interface of the fluid on the right of the medium */
+ prim_interfaces[7] = interfaces[2];
+ prim_interfaces[8] = interfaces[2];
+ prim_interfaces[9] = interfaces[4];
+
+ /* Create the scene */
+ geom.positions = vertices;
+ geom.indices = indices;
+ geom.interfaces = prim_interfaces;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, nvertices,
+ get_position, &geom, &scn) == RES_OK);
+
+ hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
+ tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0);
+ Ts0 = T0 + tmp;
+ Ts1 = T1 - tmp;
+
+ /* Run the simulations */
+ CHK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng) == RES_OK);
+ FOR_EACH(isimul, 0, nsimuls) {
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator;
+ double pos[2];
+ double ref, u;
+ size_t nreals = 0;
+ size_t nfails = 0;
+
+ pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
+ pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
+
+ CHK(sdis_solve_probe(scn, 10000, pos, INF, 1, -1, Tref, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+
+ u = (pos[0] + 1) / thickness;
+ ref = u * Ts1 + (1-u) * Ts0;
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+
+ CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+
+ /* Release memory */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[0]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[1]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[2]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[3]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[4]) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+ CHK(ssp_rng_ref_put(rng) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+
+ return 0;
+}
+
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -16,6 +16,8 @@
#ifndef TEST_SDIS_UTILS_H
#define TEST_SDIS_UTILS_H
+#include "sdis.h"
+
#include <rsys/mem_allocator.h>
#include <stdio.h>
