commit ca87a91855b78e95e096fdd4600123bfcb205438
parent 551b5ce63b058bdb4e2acebe09bd1fa94bbda6d6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 27 Mar 2024 16:07:31 +0100
Add a test for unsteady conduction
It is designed around an infinite unstable temperature profile, meaning
that at any point, at any time, we can analytically calculate the
temperature. In Stardis, we calculate the temperature in a super shape whose
boundary temperature is set to the aforementioned unstable temperature
profile. The initial temperature of the super shape is never reached
(initial time to -infinity) for two reasons: firstly, because the
non-uniform initial conditions are still not supported; and secondly, to
only test the unstable diffusion and its time travel along the sampled
trajectory, and not its premature termination.
Note that the temperature profile returns negative values, which is no
longer a problem since negative temperatures no longer mean that the
temperature is unknown (see commit 75ae9b8). This temperature profile
can be interpreted as a variation profile with respect to a given
temperature.
Diffstat:
2 files changed, 321 insertions(+), 0 deletions(-)
diff --git a/Makefile b/Makefile
@@ -200,6 +200,7 @@ TEST_SRC =\
src/test_sdis_unstationary_atm.c\
src/test_sdis_unsteady.c\
src/test_sdis_unsteady_1d.c\
+ src/test_sdis_unsteady_analytic_profile.c\
src/test_sdis_volumic_power.c\
src/test_sdis_volumic_power4.c
TEST_SRC_LONG =\
@@ -382,18 +383,21 @@ test_sdis_volumic_power4 \
src/test_sdis_draw_external_flux.d \
src/test_sdis_solid_random_walk_robustness.d \
src/test_sdis_solve_probe3.d \
+src/test_sdis_unsteady_analytic_profile.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis_draw_external_flux.o \
src/test_sdis_solid_random_walk_robustness.o \
src/test_sdis_solve_probe3.o \
+src/test_sdis_unsteady_analytic_profile.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
test_sdis_draw_external_flux \
test_sdis_solid_random_walk_robustness \
test_sdis_solve_probe3 \
+test_sdis_unsteady_analytic_profile \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS) $(S3DUT_LIBS)
diff --git a/src/test_sdis_unsteady_analytic_profile.c b/src/test_sdis_unsteady_analytic_profile.c
@@ -0,0 +1,317 @@
+/* Copyright (C) 2016-2023 |Méso|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 <star/s3dut.h>
+
+#include <rsys/mem_allocator.h>
+
+/*
+ * The system is an unsteady-state temperature profile, meaning that at any
+ * point, at any:e time, we can analytically calculate the temperature. We
+ * immerse in this temperature field a supershape representing a solid in which
+ * we want to evaluate the temperature by Monte Carlo at a given position and
+ * observation time. On the Monte Carlo side, the temperature of the supershape
+ * is unknown. Only the boundary temperature is fixed at the temperature of the
+ * unsteady trilinear profile mentioned above. Monte Carlo would then have to
+ * find the temperature defined by the unsteady profile.
+ *
+ * T(z) /\ <-- T(x,y,z,t)
+ * | T(y) ___/ \___
+ * |/ \ . T=? /
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ *
+ */
+
+#define LAMBDA 0.1 /* [W/(m.K)] */
+#define RHO 25.0 /* [kg/m^3] */
+#define CP 2.0 /* [J/K/kg)] */
+#define DELTA 1.0/20.0 /* [m/fp_to_meter] */
+
+#define NREALISATIONS 100000
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+temperature(const double pos[3], const double time)
+{
+ const double kx = PI/4.0;
+ const double ky = PI/4.0;
+ const double kz = PI/4.0;
+ const double alpha = LAMBDA / (RHO*CP); /* Diffusivity */
+
+ const double A = 0; /* No termal source */
+ const double B1 = 10;
+ const double B2 = 1000;
+
+ double x, y, z, t;
+ double a, b, c;
+ double temp;
+
+ ASSERT(pos);
+
+ x = pos[0];
+ y = pos[1];
+ z = pos[2];
+ t = time;
+
+ a = B1*(x*x*x*z-3*x*y*y*z);
+ b = B2*sin(kx*x)*sin(ky*y)*sin(kz*z)*exp(-alpha*(kx*kx + ky*ky + kz*kz)*t);
+ c = A * x*x*x*x * y*y*y * z*z;
+
+ temp = (a + b - c) / LAMBDA;
+ return temp;
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static struct s3dut_mesh*
+create_super_shape(void)
+{
+ struct s3dut_mesh* mesh = NULL;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 1;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &mesh));
+
+ return mesh;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ struct s3dut_mesh_data mesh_data;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {{0}, NULL};
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < context->mesh_data.nprimitives);
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh_data.indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh_data.indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh_data.indices[itri*3+1];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < context->mesh_data.nprimitives);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < context->mesh_data.nvertices);
+ pos[0] = context->mesh_data.positions[ivert*3+0];
+ pos[1] = context->mesh_data.positions[ivert*3+1];
+ pos[2] = context->mesh_data.positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct s3dut_mesh* mesh,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ OK(s3dut_mesh_get_data(mesh, &context.mesh_data));
+ context.interf = interf;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = context.mesh_data.nprimitives;
+ scn_args.nvertices = context.mesh_data.nvertices;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of the super shape
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, CP) /* [J/K/kg] */
+SOLID_PROP(thermal_conductivity, LAMBDA) /* [W/m/K] */
+SOLID_PROP(volumic_mass, RHO) /* [kg/m^3] */
+SOLID_PROP(delta, DELTA) /* [m/fp_to_meter] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
+#undef SOLID_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.t0 = -INF;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is fixed to the unsteady-state profile
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data;
+ return temperature(frag->P, frag->time);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe
+ (struct sdis_scene* scn,
+ const enum sdis_diffusion_algorithm diff_algo,
+ const double pos[3],
+ const double time) /* [s] */
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ args.nrealisations = NREALISATIONS;
+ args.position[0] = pos[0];
+ args.position[1] = pos[1];
+ args.position[2] = pos[2];
+ args.time_range[0] = time;
+ args.time_range[1] = time;
+ args.diff_algo = diff_algo;
+
+ OK(sdis_solve_probe(scn, &args, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = temperature(pos, time);
+ printf("T(%g, %g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), time, ref, T.E, T.SE);
+ CHK(eq_eps(ref, T.E, 3*T.SE));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ struct s3dut_mesh* super_shape = NULL;
+ const double pos[3] = {0.2,0.3,0.4}; /* [m/fp_to_meter] */
+ const double time = 5; /* [s] */
+
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ super_shape = create_super_shape();
+
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scn = create_scene(sdis, super_shape, interf);
+
+ check_probe(scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time);
+ check_probe(scn, SDIS_DIFFUSION_WOS, pos, time);
+
+ OK(s3dut_mesh_ref_put(super_shape));
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scn));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
