commit 1e1e73f8d36d51aa4d06c568f0127d4c5cb75143
parent 0333e1bf05eb29a0474a38584f04131cc10fedb8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 8 Sep 2025 15:34:45 +0200
Testing the imposed flux
Two methods were used to test the imposed flux: as a boundary condition
or as an imposed flux on a solid-fluid connection.
Diffstat:
2 files changed, 327 insertions(+), 0 deletions(-)
diff --git a/Makefile b/Makefile
@@ -126,6 +126,7 @@ TEST_SRC=\
src/sadist_conducto_radiative.c\
src/sadist_custom_conductive_path.c\
src/sadist_external_flux.c\
+ src/sadist_flux_with_h.c\
src/sadist_probe_boundary.c\
src/sadist_unsteady.c
TEST_OBJ=$(TEST_SRC:.c=.o)
@@ -163,6 +164,7 @@ $(TEST_OBJ) $(TEST_DEP): config.mk
sadist_conducto_radiative \
sadist_custom_conductive_path \
sadist_external_flux \
+sadist_flux_with_h \
sadist_probe_boundary \
sadist_unsteady \
: config.mk
diff --git a/src/sadist_flux_with_h.c b/src/sadist_flux_with_h.c
@@ -0,0 +1,325 @@
+/* Copyright (C) 2024 |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/>. */
+
+#define _POSIX_C_SOURCE 200112L /* popen */
+
+#include "sadist.h"
+
+#include <star/s3dut.h>
+#include <star/sstl.h>
+
+#include <rsys/cstr.h>
+#include <rsys/float3.h>
+#include <rsys/math.h>
+
+#include <errno.h>
+#include <stdio.h>
+#include <string.h> /* strerror */
+#include <wait.h> /* WIFEXITED, WEXITSTATUS */
+
+/*
+ * The configuration is a rectangle with conductivity lambda and unknown
+ * temperature. Its upper and lower boundaries are adiabatic, while its left and
+ * right boundaries have fixed fluxes phi1 and phi2, respectively. The left
+ * boundary also has convective exchange with the surrounding fluid, whose
+ * temperature is Text. In steady state, the temperature at a given position in
+ * the solid rectangle is:
+ *
+ * T(x) = phi2/lambda*x + (Text + (phi1 + phi2)/h)
+ *
+ * with h the convective coefficient on the left boundary
+ *
+ * Text /////// (1,1)
+ * +---------+
+ * | |
+ * h _\ |--> <--|
+ * / / phi1-> <-phi2
+ * \__/ |--> <--|
+ * | |
+ * +---------+
+ * Y (0,0) ///////
+ * |__ X
+ * /
+ * Z
+ *
+ * This test is performed with different system descriptions in order to verify
+ * how Stardis processes them. Interfaces with imposed flows are defined either
+ * as boundary conditions or as fluid-solid connections with a surrounding fluid
+ * whose temperature is known. The objective is to verify that Stardis correctly
+ * processes these different descriptions.
+ */
+
+#define FILENAME1 "scene.txt"
+#define FILENAME2 "scene2.txt"
+
+/* Geometry */
+#define CUBOID_LEN_X 0.2 /* [m] */
+#define CUBOID_LEN_Y 0.5 /* [m] */
+#define CUBOID_LEN_Z 1.0 /* [m] */
+
+/* Physical properties */
+#define LAMBDA 25.0 /* [W.m^-1.K^-1] */
+#define RHO 7500.0 /* [kg.m^-3] */
+#define CP 500.0 /* [J.K^-1.kg^-1] */
+#define Text 373.15 /* [K] */
+#define PHI1 1000.0 /* [W.m^-2] */
+#define PHI2 5000.0 /* [W.m^-2] */
+#define H 10.0 /* [W.K^-1.m^-2] */
+
+/* Probe position */
+#define PX 0.1
+#define PY 0.25
+#define PZ 0.0
+
+/* The commands to be executed */
+#define COMMAND "stardis -a wos -V3 -p"STR(PX)","STR(PY)","STR(PZ)" -n1000 -M"
+
+enum facet_flag {
+ LX = BIT(0), /* Lower X */
+ UX = BIT(1), /* Upper X */
+ LY = BIT(2), /* Lower Y */
+ UY = BIT(3), /* Upper Y */
+ LZ = BIT(4), /* Lower Z */
+ UZ = BIT(5), /* Upper Z */
+ X = LX|UX, /* Lower and upper X */
+ Y = LY|UY, /* Lower and upper Y */
+ Z = LZ|UZ /* Lower and upper Z */
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static enum facet_flag
+get_facet_flag(const struct sstl_facet* facet)
+{
+ ASSERT(facet);
+
+ if(facet->vertices[0][0] == facet->vertices[1][0]
+ && facet->vertices[0][0] == facet->vertices[2][0]) {
+ return facet->vertices[0][0] > 0 ? UX : LX;
+ }
+
+ if(facet->vertices[0][1] == facet->vertices[1][1]
+ && facet->vertices[0][1] == facet->vertices[2][1]) {
+ return facet->vertices[0][1] > 0 ? UY : LY;
+ }
+
+ if(facet->vertices[0][2] == facet->vertices[1][2]
+ && facet->vertices[0][2] == facet->vertices[2][2]) {
+ return facet->vertices[0][2] > 0 ? UY : LY;
+ }
+
+ FATAL("Unreachable code\n");
+ return 0;
+}
+
+static res_T
+write_faces
+ (FILE* fp,
+ const struct s3dut_mesh* mesh,
+ const unsigned int facet_mask)
+{
+ struct sstl_writer_create_args args = SSTL_WRITER_CREATE_ARGS_DEFAULT__;
+ struct sstl_writer* writer = NULL;
+ struct s3dut_mesh_data data;
+ size_t i = 0;
+ res_T res = RES_OK;
+ ASSERT(fp && mesh);
+
+ S3DUT(mesh_get_data(mesh, &data));
+
+ args.filename = "Cuboid facets";
+ args.stream = fp;
+ args.type = SSTL_ASCII;
+ args.verbose = 1;
+ if((res = sstl_writer_create(&args, &writer)) != RES_OK) goto error;
+
+ FOR_EACH(i, 0, data.nprimitives) {
+ struct sstl_facet facet = SSTL_FACET_NULL;
+ const float offset[3] = {(float)PX*0.5f, (float)PY*0.5f, (float)PZ*0.5f};
+ enum facet_flag flag = 0;
+
+ f3_set_d3(facet.vertices[0], data.positions + data.indices[i*3+0] * 3);
+ f3_set_d3(facet.vertices[1], data.positions + data.indices[i*3+1] * 3);
+ f3_set_d3(facet.vertices[2], data.positions + data.indices[i*3+2] * 3);
+ f3_add(facet.vertices[0], facet.vertices[0], offset);
+ f3_add(facet.vertices[1], facet.vertices[1], offset);
+ f3_add(facet.vertices[2], facet.vertices[2], offset);
+ flag = get_facet_flag(&facet);
+
+ if(!(flag & facet_mask)) continue; /* Discard the facet */
+
+ if((res = sstl_write_facet(writer, &facet)) != RES_OK) goto error;
+ }
+
+exit:
+ if(writer) SSTL(writer_ref_put(writer));
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_scene(FILE* fp)
+{
+ ASSERT(fp);
+ fprintf(fp, "SOLID cuboid %f %f %f AUTO 300 UNKNOWN 0 BACK cuboid.stl\n",
+ LAMBDA, RHO, CP);
+ fprintf(fp, "H_BOUNDARY_FOR_SOLID Adiabatic 0 0 0 0 0 cuboid_yYzZ.stl\n");
+ fprintf(fp, "HF_BOUNDARY_FOR_SOLID HFlux 0 0 0 %f %f %f cuboid_x.stl\n",
+ H, Text, PHI1);
+ fprintf(fp, "F_BOUNDARY_FOR_SOLID Flux %f cuboid_X.stl\n", PHI2);
+ return RES_OK;
+}
+
+static res_T
+setup_scene2(FILE* fp)
+{
+ fprintf(fp, "SOLID cuboid %f %f %f AUTO 300 UNKNOWN 0 BACK cuboid.stl\n",
+ LAMBDA, RHO, CP);
+ fprintf(fp, "FLUID Extern 1 1 %f %f FRONT cuboid.stl\n", Text, Text);
+ fprintf(fp, "SOLID_FLUID_CONNECTION Adiabatic 0 0 0 0 cuboid_yYzZ.stl\n");
+ fprintf(fp, "F_SOLID_FLUID_CONNECTION HFlux 0 0 0 %f %f cuboid_x.stl\n", H, PHI1);
+ fprintf(fp, "F_SOLID_FLUID_CONNECTION Flux 0 0 0 0 %f cuboid_X.stl\n", PHI2);
+ return RES_OK;
+}
+
+static int
+init(void)
+{
+ struct s3dut_mesh* mesh = NULL;
+ FILE* fp_cuboid = NULL;
+ FILE* fp_cuboid_x = NULL;
+ FILE* fp_cuboid_X = NULL;
+ FILE* fp_cuboid_yYzZ = NULL;
+ FILE* fp_scene1 = NULL;
+ FILE* fp_scene2 = NULL;
+ res_T res = RES_OK;
+ int err = 0;
+
+ #define FOPEN(Name) { \
+ if(!(fp_##Name = fopen(STR(Name)".stl", "w"))) { \
+ perror(STR(Name)".stl"); \
+ goto error; \
+ } \
+ } (void)0
+ FOPEN(cuboid);
+ FOPEN(cuboid_x);
+ FOPEN(cuboid_X);
+ FOPEN(cuboid_yYzZ);
+ #undef FOPEN
+
+ if(!(fp_scene1 = fopen(FILENAME1, "w"))) { perror(FILENAME1); goto error; }
+ if(!(fp_scene2 = fopen(FILENAME2, "w"))) { perror(FILENAME2); goto error; }
+
+ res = s3dut_create_cuboid(NULL, CUBOID_LEN_X, CUBOID_LEN_Y, CUBOID_LEN_Z, &mesh);
+ if(res != RES_OK) {
+ fprintf(stderr, "Error creating the cuboid -- %s\n", res_to_cstr(res));
+ goto error;
+ }
+
+ if((res = write_faces(fp_cuboid, mesh, X|Y|Z)) != RES_OK) goto error;
+ if((res = write_faces(fp_cuboid_x, mesh, LX)) != RES_OK) goto error;
+ if((res = write_faces(fp_cuboid_X, mesh, UX)) != RES_OK) goto error;
+ if((res = write_faces(fp_cuboid_yYzZ, mesh, Y|Z)) != RES_OK) goto error;
+ if((res = setup_scene(fp_scene1)) != RES_OK) goto error;
+ if((res = setup_scene2(fp_scene2)) != RES_OK) goto error;
+
+exit:
+ if(mesh) S3DUT(mesh_ref_put(mesh));
+ if(fp_cuboid) fclose(fp_cuboid);
+ if(fp_cuboid_x) fclose(fp_cuboid_x);
+ if(fp_cuboid_X) fclose(fp_cuboid_X);
+ if(fp_scene1) fclose(fp_scene1);
+ if(fp_scene2) fclose(fp_scene2);
+ return err;
+error:
+ err = 1;
+ goto exit;
+}
+
+static double
+analytical_solution(void)
+{
+ return PHI2 / LAMBDA * PX + (Text + (PHI1 + PHI2)/H);
+}
+
+static int
+run(const char* command)
+{
+ FILE* output = NULL;
+ double E = 0;
+ double SE = 0;
+ double ref = 0;
+ int n = 0;
+ int err = 0;
+ int status = 0;
+
+ printf("%s\n", command);
+
+ if(!(output = popen(command, "r"))) {
+ fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno));
+ err = errno;
+ goto error;
+ }
+
+ if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) {
+ fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno));
+ err = errno;
+ goto error;
+ }
+
+ /* Check command exit status */
+ if((status=pclose(output)), output=NULL, status) {
+ if(status == -1) err = errno;
+ else if(WIFEXITED(status)) err = WEXITSTATUS(status);
+ else if(WIFSIGNALED(status)) err = WTERMSIG(status);
+ else if(WIFSTOPPED(status)) err = WSTOPSIG(status);
+ else FATAL("Unreachable code.\n");
+ goto error;
+ }
+
+ ref = analytical_solution();
+ printf("T = %g ~ %g +/- %g\n", ref, E, SE);
+ if(!eq_eps(ref, E, SE*3)) {
+ err = 1;
+ goto error;
+ }
+
+exit:
+ if(output) pclose(output);
+ return err;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ int err = 0;
+ (void)argc, (void)argv;
+
+ if((err = init())) goto error;
+ if((err = run(COMMAND STR(FILENAME1)))) goto error;
+ if((err = run(COMMAND STR(FILENAME2)))) goto error;
+
+exit:
+ return err;
+error:
+ goto exit;
+}
