commit 373fa3f5c56e5e94f0436708275b4e45c61bcea7
parent 00a9259698486399b75ce42878b5feef6b7a2c6f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 14 Feb 2024 17:18:13 +0100
In-depth rewriting of the compute_probe_on_interface function
This commit is a rereading of the function with a restructuring and
rewriting of its code to make it easier to read and reduce its technical
debt. It's a commit that prepares the resolution of several probes in
parallel: its code should be shared with this upcoming feature.
Temperature calculation at the boundary is now based on this new
function. The calculation at the boundary is therefore no longer
broken... except by bugs. Indeed, although the code compiles without
warning or error, nothing has yet been tested. This rewrite must still
be very buggy! The rest of the message gives a brief summary of the
revision.
The function has been moved to a separate file. From now on, nothing is
shared with the calculation of a probe point in a medium. For example,
the check_probe_conform_to_type function did too many things, processing
both a probe in a medium and on a boundary. It was therefore very
long, complicated and difficult to read. In this commit, it's still
hard to read, but it now only handles probes in the medium, which
simplifies it a little. What it used to do for a boundary probe is now
moved to the move_to_boundary helper function, which is a heavily
reworked version of the previous implementation. Among other things, it
has been completely restructured, relying on sub-functions that make the
source code easier to read. In addition, we now rely on the solver to
project the position of the probe onto the boundary, so as to share the
same structure and avoid having to build a new one from scratch. The
filter function has also been simplified; the old filter performed
several calculations and data accesses that could be performed once and
for all after traversing the structure.
The same work was carried out throughout the function review. For
example, new functions were written to read/write the state of the RNG;
the use of macros was a source of errors while adding nothing in
particular except complexity. The estimation of a temperature or a
green function is now done in two separate functions. The side of the
probe is also defined in its own function. And so on.
Diffstat:
9 files changed, 924 insertions(+), 435 deletions(-)
diff --git a/Makefile b/Makefile
@@ -30,6 +30,7 @@ SRC =\
src/stardis-app.c\
src/stardis-args.c\
src/stardis-compute.c\
+ src/stardis-compute-probe-boundary.c\
src/stardis-description.c\
src/stardis-fluid.c\
src/stardis-fluid-prog.c\
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -319,8 +319,8 @@ stardis_init
else if(args->mode & MODE_MEDIUM_COMPUTE) {
ERR(str_set(&stardis->solve_name, args->medium_name));
}
- else if(args->mode & MODE_PROBE_COMPUTE_ON_INTERFACE && args->medium_name) {
- ERR(str_set(&stardis->solve_name, args->medium_name));
+ else if(args->mode & MODE_PROBE_COMPUTE_ON_INTERFACE) {
+ stardis->probe_boundary = args->probe_boundary;
}
else if(args->mode & SURFACE_COMPUTE_MODES) {
ERR(str_set(&stardis->solve_name, args->solve_filename));
@@ -484,7 +484,7 @@ stardis_init
/* If computation is on a volume, check medium is known */
if(args->mode & MODE_MEDIUM_COMPUTE) {
- if(!find_medium_by_name(stardis, &stardis->solve_name, NULL)) {
+ if(!find_medium_by_name(stardis, str_cget(&stardis->solve_name), NULL)) {
logger_print(stardis->logger, (is_for_compute ? LOG_ERROR : LOG_WARNING),
"Cannot find medium '%s'\n", str_cget(&stardis->solve_name));
res = RES_BAD_ARG;
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -16,6 +16,7 @@
#ifndef STARDIS_APP_H
#define STARDIS_APP_H
+#include "stardis-args.h"
#include "stardis-description.h"
#include <star/sg3d.h>
@@ -238,6 +239,8 @@ struct stardis {
int geometry_initialized;
int mpi_initialized;
int mpi_rank;
+
+ struct stardis_probe_boundary probe_boundary;
};
unsigned
diff --git a/src/stardis-compute-probe-boundary.c b/src/stardis-compute-probe-boundary.c
@@ -0,0 +1,832 @@
+/* Copyright (C) 2018-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 "stardis-app.h"
+#include "stardis-compute.h"
+#include "stardis-output.h"
+#include "stardis-prog-properties.h"
+#include "stardis-solid.h"
+#include "stardis-solid-prog.h"
+#include "stardis-ssconnect.h"
+
+#include <sdis.h>
+
+#include <star/senc3d.h>
+#include <star/ssp.h>
+
+#include <rsys/logger.h>
+#include <rsys/str.h>
+#include <rsys/clock_time.h>
+
+#include <strings.h>
+
+struct filter_ctx {
+ float distance;
+ enum sg3d_property_type side;
+};
+#define FILTER_CTX_DEFAULT__ {0.f, SG3D_INTFACE}
+static const struct filter_ctx FILTER_CTX_DEFAULT = FILTER_CTX_DEFAULT__;
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE const char*
+sdis_side_to_cstr(const enum sdis_side side)
+{
+ const char* cstr = NULL;
+ switch(side) {
+ case SDIS_FRONT: cstr = "FRONT"; break;
+ case SDIS_BACK: cstr = "BACK"; break;
+ case SDIS_SIDE_NULL__: cstr = "UNDEFINED"; break;
+ default: FATAL("Unreachable code.\n");
+ }
+ return cstr;
+}
+
+static res_T
+read_rng_state
+ (struct stardis* stardis,
+ const char* filename,
+ struct ssp_rng* rng)
+{
+ FILE* fp = NULL;
+ res_T res = RES_OK;
+ ASSERT(stardis && filename && rng);
+
+ if(stardis->mpi_initialized && stardis->mpi_rank != 0) {
+ goto exit; /* Non master process. Nothing to do */
+ }
+
+ if((fp = fopen(filename, "r")) == NULL) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot open generator's state file ('%s').\n", filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ res = ssp_rng_read(rng, fp);
+ if(res != RES_OK) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot read random generator's state ('%s').\n", filename);
+ goto error;
+ }
+
+exit:
+ if(fp) CHK(fclose(fp) == 0);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+write_rng_state
+ (struct stardis* stardis,
+ const char* filename,
+ struct ssp_rng* rng_state)
+{
+ FILE* fp = NULL;
+ res_T res = RES_OK;
+ ASSERT(stardis && filename && rng_state);
+
+ if(stardis->mpi_initialized && stardis->mpi_rank != 0) {
+ goto exit; /* Non master process. Nothing to do */
+ }
+
+ if((fp = fopen(filename, "wb")) == NULL) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot open generator's state file ('%s').\n", filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ res = ssp_rng_write(rng_state, fp);
+ if(res != RES_OK) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot write random generator's state ('%s').\n", filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+exit:
+ if(fp) CHK(fclose(fp) == 0);
+ return res;
+error:
+ goto exit;
+}
+
+/* Filter used from a point query to determine not only one of the closest
+ * point, but the better one if there are more than one. In some circumstances
+ * it is not possible to determine the medium we are in using a given hit, but
+ * it is possible using another equidistant hit :
+ *
+ *
+ * P C
+ * +.............+---trg 1---
+ * |
+ * medium 1 trg 2 medium 2
+ * |
+ *
+ * C is the closest point from P, and 2 different hits at C are possible (one
+ * on each triangle). However, only hit on trg 2 allows to find out that P is
+ * in medium 1 using sign(PC.Ntrg1) as PC.Ntrg2 = 0.
+ * The following filter function aims at selecting the hit on trg2 regardless
+ * of the order in which the 2 triangles are checked.
+ * One unexpected case cannot be decided though, but it implies that the
+ * closest triangle has 2 different media on its sides and that P lies on the
+ * triangle's plane :
+ *
+ * P C medium 1
+ * + +---trg---
+ * medium 2 */
+static int
+hit_filter
+ (const struct s3d_hit* hit,
+ const float ray_org[3],
+ const float ray_dir[3],
+ const float ray_range[2],
+ void* ray_data,
+ void* filter_data)
+{
+ struct filter_ctx* filter_ctx = ray_data;
+
+ (void)ray_org, (void)ray_range, (void)filter_data;
+ ASSERT(hit && filter_ctx);
+ ASSERT(hit->uv[0] == CLAMP(hit->uv[0], 0, 1));
+ ASSERT(hit->uv[1] == CLAMP(hit->uv[1], 0, 1));
+
+ /* That's not the closest point. Keep the previous one if it can be used to
+ * detect the medium (i.e. side != SG3D_INTFACE) */
+ if(filter_ctx->distance == hit->distance && filter_ctx->side != SG3D_INTFACE) {
+ return 1; /* Skip */
+ }
+
+ filter_ctx->distance = hit->distance;
+
+ if(filter_ctx->distance == 0) {
+ filter_ctx->side = SG3D_INTFACE;
+ } else {
+ float sign = 0;
+ float N[3] = {0,0,0}; /* Normalized normal */
+
+ /* Calculate the dot product with normalized vectors limits the numerical
+ * inaccuracies on its sign */
+ f3_normalize(N, hit->normal);
+ sign = f3_dot(ray_dir, N);
+
+ /* Star3D hit normals are left-handed */
+ if(sign < 0) filter_ctx->side = SG3D_FRONT;
+ else if(sign > 0) filter_ctx->side = SG3D_BACK;
+ else/*sign == 0*/ filter_ctx->side = SG3D_INTFACE;
+ }
+
+ return 0; /* Keep */
+}
+
+static INLINE res_T
+find_closest_point
+ (struct stardis* stardis,
+ const double pos[3],
+ struct filter_ctx* filter_ctx,
+ size_t* iprim,
+ double uv[2])
+{
+ struct sdis_scene_find_closest_point_args closest_pt_args =
+ SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
+ res_T res = RES_OK;
+ ASSERT(stardis && pos && filter_ctx && iprim && uv);
+
+ /* Find the surface point closest to the input position */
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.position[2] = pos[2];
+ closest_pt_args.radius = INF;
+ closest_pt_args.filter_3d = hit_filter;
+ closest_pt_args.filter_data = &filter_ctx;
+ ERR(sdis_scene_find_closest_point
+ (stardis->sdis_scn , &closest_pt_args, iprim, uv));
+ if(*iprim == SDIS_PRIMITIVE_NONE) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+check_move_to_solid_boundary
+ (const struct stardis* stardis,
+ const double pos[3], /* Original position */
+ const double time, /* [s] */
+ const struct description* desc, /* Solid medium in which pos lies */
+ const size_t iprim, /* Triangle index to which to move */
+ const double uv[2], /* Triangle coordinates to which to move */
+ const double distance) /* Move distance */
+{
+ struct stardis_vertex vtx = STARDIS_VERTEX_NULL;
+ const char* prefix = "";
+ const char* solid_name = "";
+ double delta = 0;
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(stardis && pos && time > 0 && desc && uv && distance >= 0);
+
+ /* Retrieve the delta and define the prefix of the solid for log messages */
+ switch(desc->type) {
+ /* Regular solid, i.e. solid with constant properties */
+ case DESC_MAT_SOLID:
+ delta = desc->d.solid->delta;
+ prefix = "";
+ break;
+
+ /* Solid with programmed properties */
+ case DESC_MAT_SOLID_PROG:
+ vtx.P[0] = pos[0];
+ vtx.P[1] = pos[1];
+ vtx.P[2] = pos[2];
+ vtx.time = time;
+ delta = desc->d.solid_prog->delta(&vtx, desc->d.solid_prog->prog_data);
+ prefix = "programmed";
+ break;
+
+ default: FATAL("Unreachable code.\n");
+ }
+
+ solid_name = str_cget(get_description_name(desc));
+ logger_print(stardis->logger, LOG_OUTPUT, "Probe was in %ssolid '%s'.\n",
+ prefix, solid_name);
+
+ /* The position is closed from the triangle */
+ if(distance < 0.5*delta) {
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe was %g delta from closest boundary.\n",
+ distance/delta);
+
+ /* Notice that the position is a little far from the triangle */
+ } else if(distance < 2*delta) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe was %g delta from closest boundary. "
+ "Consider using -p instead of -P.\n",
+ distance/delta);
+
+ /* The position is too far from the triangle */
+ } else {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Probe moved to (%g, %g, %g), primitive %lu, uv = (%g, %g). "
+ "Move is %g delta long. Use -p instead of -P.\n",
+ SPLIT3(pos), (unsigned long)iprim, SPLIT2(uv), distance/delta);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/* This function checks nothing. It only records the status. It is named as the
+ * one used to check the projection on the solid limit to make it symmetrical,
+ * and thus simplify the reading of sources */
+static res_T
+check_move_to_fluid_boundary
+ (struct stardis* stardis,
+ const struct description* desc, /* Fluid medium in which pos lies */
+ const double distance) /* Move distance */
+{
+ const char* prefix = "";
+ const char* fluid_name = "";
+
+ ASSERT(stardis && desc && distance >= 0);
+
+ switch(desc->type) {
+ case DESC_MAT_FLUID: prefix = ""; break;
+ case DESC_MAT_FLUID_PROG: prefix = "programmed"; break;
+ default: FATAL("Unreachable code.\n");
+ }
+
+ fluid_name = str_cget(get_description_name(desc));
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe was in %sfluid '%s'.\n", prefix, fluid_name);
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe distance from closest boundary was %g.\n", distance);
+
+ return RES_OK;
+}
+
+static res_T
+move_to_boundary
+ (struct stardis* stardis,
+ const double pos[3],
+ const double time, /* [s] */
+ size_t* out_iprim,
+ double uv[2])
+{
+ /* Position on boundary */
+ struct filter_ctx filter_ctx = FILTER_CTX_DEFAULT;
+ double proj_pos[3] = {0,0,0};
+ size_t iprim = 0;
+
+ /* Properties */
+ const struct description* desc_list = NULL;
+ const struct description* desc = NULL;
+ unsigned desc_ids[SG3D_PROP_TYPES_COUNT__];
+
+ /* Miscellaneous */
+ size_t nvertices_close = 0;
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(stardis && pos && time >= 0 && out_iprim && uv);
+
+ ERR(find_closest_point(stardis, pos, &filter_ctx, &iprim, uv));
+ SG3D(geometry_get_unique_triangle_properties
+ (stardis->geometry.sg3d, (unsigned)iprim, desc_ids));
+
+ desc_list = darray_descriptions_cdata_get(&stardis->descriptions);
+
+ /* Undefined medium */
+ if(filter_ctx.side == SG3D_INTFACE
+ || desc_ids[filter_ctx.side] == SG3D_UNSPECIFIED_PROPERTY) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Could not determine the medium probe is in.\n");
+ }
+
+ if(filter_ctx.side != SG3D_INTFACE) {
+
+ /* Probe is outside the system */
+ if(desc_ids[filter_ctx.side] == SG3D_UNSPECIFIED_PROPERTY) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe was outside the system.\n");
+
+ /* Probe is in a medium */
+ } else {
+ desc = desc_list + desc_ids[filter_ctx.side];
+
+ switch(desc->type) {
+ case DESC_MAT_SOLID:
+ case DESC_MAT_SOLID_PROG:
+ ERR(check_move_to_solid_boundary
+ (stardis, pos, time, desc, iprim, uv, filter_ctx.distance));
+ break;
+ case DESC_MAT_FLUID:
+ case DESC_MAT_FLUID_PROG:
+ ERR(check_move_to_fluid_boundary
+ (stardis, desc, filter_ctx.distance));
+ break;
+ default: FATAL("Unreachable code.\n");
+ }
+ }
+ }
+
+ SDIS(scene_get_boundary_position(stardis->sdis_scn, iprim, uv, proj_pos));
+
+ /* Count the number of vertices that are close to the boundary position
+ * and issue a warning if necessary */
+ nvertices_close += CLAMP(uv[0], 0.0005, 0.9995) != uv[0];
+ nvertices_close += CLAMP(uv[1], 0.0005, 0.9995) != uv[1];
+ if(nvertices_close) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe %s close to / on %s. "
+ "If computation fails, try moving it slightly.\n",
+ filter_ctx.distance == 0 ? "is" : "moved",
+ nvertices_close == 1 ? "an edge" : "a vertex");
+ }
+
+ /* Probe is on a boundary */
+ if(filter_ctx.distance == 0) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Probe is on primitive %lu, uv = (%g, %g), not moved.\n",
+ (unsigned long)iprim, SPLIT2(uv));
+
+ /* Probe was projected on a boundary */
+ } else {
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe moved to (%g, %g, %g), primitive %lu, uv = (%g, %g).\n",
+ SPLIT3(proj_pos), (unsigned long)iprim, SPLIT2(uv));
+ }
+
+exit:
+ *out_iprim = iprim;
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_probe_side
+ (struct stardis* stardis,
+ const unsigned desc_ids[SG3D_PROP_TYPES_COUNT__],
+ const char* side_str,
+ const size_t iprim,
+ enum sdis_side *out_side)
+{
+ const struct description* desc_list = NULL;
+ const struct description* desc_front = NULL;
+ const struct description* desc_back = NULL;
+ size_t ndescs = 0;
+ enum sdis_side side = SDIS_SIDE_NULL__;
+ res_T res = RES_OK;
+ (void)ndescs; /* Avoid "Unused variable" warnings in release */
+
+ /* Check pre-conditions */
+ ASSERT(stardis && side_str && desc_ids && out_side);
+
+ /* Fetch the properties */
+ desc_list = darray_descriptions_cdata_get(&stardis->descriptions);
+ ndescs = darray_descriptions_size_get(&stardis->descriptions);
+ desc_front = desc_list + desc_ids[SG3D_FRONT];
+ desc_back = desc_list + desc_ids[SG3D_BACK];
+
+ /* No side specified */
+ if(!side_str || !strlen(side_str)) {
+ side = SDIS_SIDE_NULL__;
+
+ /* Set probe to front side */
+ } else if(!strcasecmp(side_str, "FRONT")) {
+ ASSERT(desc_ids[SG3D_FRONT] < ndescs && DESC_IS_MEDIUM(desc_front));
+ side = SDIS_FRONT;
+
+ /* Set probe to back side */
+ } else if(!strcasecmp(side_str, "BACK")) {
+ ASSERT(desc_ids[SG3D_BACK] < ndescs && DESC_IS_MEDIUM(desc_back));
+ side = SDIS_BACK;
+
+ /* Set the probe to the side that points to the submitted medium name */
+ } else {
+ unsigned med_id_probe = 0; /* Medium defined on the probe */
+ unsigned med_id_front = 0; /* Medium on front side */
+ unsigned med_id_back = 0; /* Medium on back side */
+ ASSERT(DESC_IS_MEDIUM(desc_front) && DESC_IS_MEDIUM(desc_back));
+
+ if(!find_medium_by_name(stardis, side_str, &med_id_probe)) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot locate side from medium name '%s' (unknown medium)\n",
+ side_str);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ description_get_medium_id(desc_front, &med_id_front);
+ description_get_medium_id(desc_back, &med_id_back);
+
+ /* Invalid probe medium wrt the boundary on which it is located */
+ if(med_id_probe != med_id_front
+ && med_id_probe != med_id_back) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Medium '%s' is not used at this interface (prim id=%lu)\n",
+ side_str, (unsigned long)iprim);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* The same medium is used on both sides: cannot differentiate */
+ if(med_id_front == med_id_back) {
+ unsigned encs[2]; /* Identifier of the enclosures */
+
+ ERR(senc3d_scene_get_triangle_enclosures
+ (stardis->senc3d_scn, (unsigned)iprim, encs));
+ logger_print(stardis->logger, LOG_ERROR,
+ "Medium '%s' is used on both sides of this interface (prim id=%lu).\n",
+ side_str, (unsigned long)iprim);
+ logger_print(stardis->logger, LOG_ERROR,
+ "Side must be defined using either FRONT or BACK.\n");
+ logger_print(stardis->logger, LOG_ERROR,
+ "FRONT side is related to enclosure %u, BACK side to enclosure %u.\n",
+ encs[SENC3D_FRONT], encs[SENC3D_BACK]);
+
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ side = med_id_probe == med_id_front ? SDIS_FRONT : SDIS_BACK;
+ }
+
+exit:
+ *out_side = side;
+ return res;
+error:
+ side = SDIS_SIDE_NULL__;
+ goto exit;
+}
+
+/* This function checks the conformity between the potential thermal contact
+ * resistance at the probe location and the specified probe side. */
+static res_T
+setup_thermal_contact_resistance
+ (struct stardis* stardis,
+ const unsigned desc_ids[SG3D_PROP_TYPES_COUNT__],
+ const enum sdis_side probe_side)
+{
+ struct str log_msg;
+ const struct description* desc_list = NULL;
+ const struct description* desc_front = NULL;
+ const struct description* desc_back = NULL;
+ const struct description* desc_intface = NULL;
+ size_t ndescs = 0;
+ double tcr = 0;
+ res_T res = RES_OK;
+ (void)ndescs; /* Avoid "Unused variable" warnings in release */
+
+ /* Check pre-conditions */
+ ASSERT(stardis && desc_ids);
+
+ str_init(stardis->allocator, &log_msg);
+
+ /* Fetch the properties */
+ desc_list = darray_descriptions_cdata_get(&stardis->descriptions);
+ ndescs = darray_descriptions_size_get(&stardis->descriptions);
+ desc_front = desc_list + desc_ids[SG3D_FRONT];
+ desc_back = desc_list + desc_ids[SG3D_BACK];
+ desc_intface = desc_list + desc_ids[SG3D_INTFACE];
+
+ /* Get the thermal contact resistance between solid/solid connection if any */
+ if(desc_ids[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY
+ && desc_intface->type == DESC_SOLID_SOLID_CONNECT) {
+ ASSERT(desc_ids[SG3D_INTFACE] < ndescs);
+ tcr = desc_intface->d.ss_connect->tcr;
+ }
+
+ /* Warn if side defined and no resistance defined */
+ if(tcr == 0 && probe_side != SDIS_SIDE_NULL__) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Specifying a compute side at an interface with no contact resistance "
+ "is meaningless.\n");
+ }
+
+ #define GET_DESC_NAME(Desc) str_cget(get_description_name(Desc))
+
+ /* A thermal contact resistance cannot be defined if probe side is NULL */
+ if(tcr != 0 && probe_side == SDIS_SIDE_NULL__) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot let probe computation side unspecified on an interface with a "
+ "non-nul thermal resistance.\n");
+
+ /* Format the log string */
+ if(desc_ids[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) {
+ ASSERT(desc_ids[SG3D_FRONT] < ndescs);
+ ERR(str_append_printf(&log_msg, " FRONT: '%s'", GET_DESC_NAME(desc_front)));
+ }
+ if(desc_ids[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY
+ && desc_ids[SG3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) {
+ ERR(str_append_char(&log_msg, ','));
+ }
+ if(desc_ids[SG3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) {
+ ASSERT(desc_ids[SG3D_BACK] < ndescs);
+ ERR(str_append_printf(&log_msg, " BACK: '%s'", GET_DESC_NAME(desc_back)));
+ }
+
+ /* Print error message */
+ logger_print(stardis->logger, LOG_ERROR,
+ "Interface '%s',%s, resistance=%g K.m^2/W.\n",
+ GET_DESC_NAME(desc_intface), str_cget(&log_msg), tcr);
+
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Log that a calculation is done on a boundary with tcr */
+ if(tcr > 0) {
+ const char* medium_name = probe_side == SDIS_FRONT
+ ? GET_DESC_NAME(desc_front)
+ : GET_DESC_NAME(desc_back);
+
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe computation on an interface with a thermal resistance = %g K.m^2/W "
+ "on %s side (medium is '%s').\n",
+ tcr, sdis_side_to_cstr(probe_side), medium_name);
+ }
+
+ #undef GET_DESC_NAME
+
+exit:
+ str_release(&log_msg);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+solve
+ (struct stardis* stardis,
+ struct time* start,
+ struct sdis_solve_probe_boundary_args* args,
+ res_T output[2])
+{
+ struct time t0, t1;
+ struct dump_path_context ctx = DUMP_PATH_CONTEXT_NULL;
+ struct sdis_estimator* estimator = NULL;
+ const struct str* rng_in = NULL;
+ const struct str* rng_out = NULL;
+ struct ssp_rng* rng = NULL;
+ int is_master_process = 0;
+ res_T res = RES_OK;
+ ASSERT(stardis && args && output);
+
+ is_master_process = !stardis->mpi_initialized || stardis->mpi_rank == 0;
+
+ rng_in = &stardis->rndgen_state_in_filename;
+ rng_out = &stardis->rndgen_state_in_filename;
+
+ /* Read RNG state from file */
+ if(!str_is_empty(rng_in)) {
+ ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &rng));
+ ERR(read_rng_state(stardis, str_cget(rng_in), rng));
+ args->rng_state = rng;
+ }
+
+ /* Run the calculation */
+ time_current(&t0);
+ ERR(sdis_solve_probe_boundary(stardis->sdis_scn, args, &estimator));
+ time_current(&t1);
+
+ /* No more to do for non master processes */
+ if(!is_master_process) goto exit;
+
+ ERR(print_computation_time(estimator, stardis, start, &t0, &t1, NULL));
+
+ /* Write outputs and save their status */
+ ctx.stardis = stardis;
+ ctx.rank = 0;
+ output[0] = print_single_MC_result(estimator, stardis, stdout);
+ output[1] = sdis_estimator_for_each_path(estimator, dump_path, &ctx);
+
+ /* Write the resulting RNG state to a file */
+ if(!str_is_empty(rng_out)) {
+ struct ssp_rng* rng_state = NULL;
+ ERR(sdis_estimator_get_rng_state(estimator, &rng_state));
+ ERR(write_rng_state(stardis, str_cget(rng_out), rng_state));
+ }
+
+exit:
+ if(estimator) SDIS(estimator_ref_put(estimator));
+ if(rng) SSP(rng_ref_put(rng));
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+solve_green
+ (struct stardis* stardis,
+ struct time* start,
+ struct sdis_solve_probe_boundary_args* args)
+{
+ struct time t0/*calculation start*/, t1/*calculation end*/, t2/*output end*/;
+ struct sdis_green_function* green = NULL;
+ FILE* fp_green = NULL;
+ FILE* fp_path = NULL;
+ const struct str* rng_in = NULL;
+ struct ssp_rng* rng = NULL;
+ int is_master_process = 0;
+ res_T res = RES_OK;
+
+ ASSERT(stardis && args);
+
+ is_master_process = !stardis->mpi_initialized || stardis->mpi_rank == 0;
+
+ rng_in = &stardis->rndgen_state_in_filename;
+
+ /* Read RNG state from file */
+ if(!str_is_empty(rng_in)) {
+ ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &rng));
+ ERR(read_rng_state(stardis, str_cget(rng_in), rng));
+ args->rng_state = rng;
+ }
+
+ /* Try to open output files to detect errors early */
+ if(is_master_process && (stardis->mode & MODE_BIN_GREEN)) {
+ const char* green_filename = str_cget(&stardis->bin_green_filename);
+ const char* path_filename = str_cget(&stardis->end_paths_filename);
+
+ if((fp_green = fopen(green_filename, "wb")) == NULL) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot open file '%s' for binary writing.\n", green_filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ if(strlen(path_filename) != 0) {
+ if((fp_path = fopen(path_filename, "w")) == NULL) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot open file '%s' for writing.\n", path_filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+ }
+ }
+
+ /* Run the Green estimation */
+ time_current(&t0); /* Calculation starts */
+ ERR(sdis_solve_probe_boundary_green_function(stardis->sdis_scn, args, &green));
+ time_current(&t1); /* Calculation ends */
+
+ /* No more to do for non master processes */
+ if(is_master_process) goto exit;
+
+ /* Write ASCII Green */
+ if(stardis->mode & MODE_GREEN) {
+ ERR(dump_green_ascii(green, stardis, stdout));
+ }
+
+ /* Write binary Green */
+ if(stardis->mode & MODE_BIN_GREEN) {
+ ERR(dump_green_bin(green, stardis, fp_green));
+ if(fp_path) {
+ ERR(dump_paths_end(green, stardis, fp_path));
+ }
+ }
+
+ time_current(&t2); /* Output ends */
+
+ ERR(print_computation_time(NULL, stardis, start, &t0, &t1, &t2));
+
+ /* Note that the resulting RNG state is not written in an output file because
+ * the solver API does not provide a function to recover it. But in fact, the
+ * green function saves the RNG state after its estimation. Therefore, the API
+ * can be expected to provide such functionality soon.
+ *
+ * TODO write the RNG status of the Green function when it is available */
+
+exit:
+ if(fp_green) CHK(fclose(fp_green) == 0);
+ if(fp_path) CHK(fclose(fp_path) == 0);
+ if(green) SDIS(green_function_ref_put(green));
+ if(rng) SSP(rng_ref_put(rng));
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+compute_probe_on_interface(struct stardis* stardis, struct time* start)
+{
+ /* Probe */
+ struct sdis_solve_probe_boundary_args args
+ = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ const struct stardis_probe_boundary* probe = NULL;
+ enum sdis_side probe_side = SDIS_SIDE_NULL__;
+
+ /* Miscellaneous */
+ unsigned desc_ids[SG3D_PROP_TYPES_COUNT__] = {0};
+ double uv[2] = {0, 0};
+ size_t iprim = SIZE_MAX;
+ res_T output_status[2] = {RES_OK, RES_OK};
+ res_T res = RES_OK;
+
+ ASSERT(stardis && start && (stardis->mode & MODE_PROBE_COMPUTE_ON_INTERFACE));
+
+ probe = &stardis->probe_boundary;
+
+ /* Calculate the probe position on the boundary */
+ ERR(move_to_boundary(stardis, probe->position, probe->time[0], &iprim, uv));
+
+ ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d,
+ (unsigned)iprim, desc_ids));
+
+ ERR(setup_probe_side(stardis, desc_ids, probe->side, iprim, &probe_side));
+ ERR(setup_thermal_contact_resistance(stardis, desc_ids, probe_side));
+
+ /* Setup of solve input parameters */
+ args.nrealisations = stardis->samples;
+ args.iprim = iprim;
+ args.uv[0] = uv[0];
+ args.uv[1] = uv[1];
+ args.time_range[0] = probe->time[0];
+ args.time_range[1] = probe->time[1];
+ args.picard_order = stardis->picard_order;
+ args.side = probe_side;
+ args.register_paths = stardis->dump_paths;
+
+ /* Run the calculation */
+ if(stardis->mode & (MODE_GREEN | MODE_BIN_GREEN)) {
+ ERR(solve_green(stardis, start, &args));
+ } else {
+ ERR(solve(stardis, start, &args, output_status));
+ }
+
+exit:
+ res = (res != RES_OK ? res : output_status[0]);
+ res = (res != RES_OK ? res : output_status[1]);
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -192,7 +192,6 @@ hit_filter
static res_T
check_probe_conform_to_type
(const struct stardis* stardis,
- const int move2boundary,
double pos[3],
double time,
unsigned* iprim,
@@ -207,7 +206,7 @@ check_probe_conform_to_type
struct s3d_vertex_data attribs;
struct filter_ctx filter_ctx = FILTER_CTX_DEFAULT__;
float origin[3];
- unsigned vcount, tcount, j;
+ unsigned vcount, tcount;
ASSERT(stardis && pos && iprim && uv);
@@ -244,175 +243,80 @@ check_probe_conform_to_type
}
ASSERT(filter_ctx.dist == hit.distance);
- if(move2boundary) {
- /* Need to move probe to closest boundary */
- int danger = 0;
- if(filter_ctx.outside) {
- /* Not an error as probe moved to boundary */
- logger_print(stardis->logger, LOG_WARNING,
- "Probe was outside the model.\n");
- }
- if(!filter_ctx.desc) {
- /* Not an error as probe moved to boundary */
- if(!filter_ctx.probe_on_boundary && !filter_ctx.outside)
- logger_print(stardis->logger, LOG_WARNING,
- "Could not determine the medium probe is in.\n");
- } else {
- if(DESC_IS_SOLID(filter_ctx.desc)) {
- double delta;
- const char* pppp;
- if(filter_ctx.desc->type == DESC_MAT_SOLID) {
- struct solid* solid = filter_ctx.desc->d.solid;
- delta = solid->delta;
- pppp = "";
- } else {
- struct solid_prog* solid_prog = filter_ctx.desc->d.solid_prog;
- struct stardis_vertex vtx;
- ASSERT(filter_ctx.desc->type == DESC_MAT_SOLID_PROG);
- d3_set(vtx.P, pos);
- vtx.time = time;
- delta = solid_prog->delta(&vtx, solid_prog->prog_data);
- pppp = "programmed ";
- }
- ASSERT(delta < INF);
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe was in %ssolid '%s'.\n",
- pppp, str_cget(get_description_name(filter_ctx.desc)));
- if(filter_ctx.dist > 2 * delta) {
- logger_print(stardis->logger, LOG_ERROR,
- "Probe moved to (%g, %g, %g), primitive %u, uv = (%g, %g).\n",
- SPLIT3(filter_ctx.pos), hit.prim.prim_id, SPLIT2(hit.uv));
- logger_print(stardis->logger, LOG_ERROR,
- "Move is %g delta long. Use -p instead of -P.\n",
- filter_ctx.dist / delta);
- res = RES_BAD_ARG;
- goto error;
- }
- if(filter_ctx.dist > 0.5 * delta) {
- logger_print(stardis->logger, LOG_WARNING,
- "Probe was %g delta from closest boundary. "
- "Consider using -p instead of -P.\n",
- filter_ctx.dist / delta);
- } else {
- if(filter_ctx.dist != 0)
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe was %g delta from closest boundary.\n",
- filter_ctx.dist / delta);
- }
- } else {
- const char* pppp;
- ASSERT(DESC_IS_FLUID(filter_ctx.desc));
- /* TODO: check move length wrt local geometry? */
- if(filter_ctx.desc->type == DESC_MAT_FLUID) {
- pppp = "";
- } else {
- ASSERT(filter_ctx.desc->type == DESC_MAT_FLUID_PROG);
- pppp = "programmed ";
- }
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe was in %sfluid '%s'.\n",
- pppp, str_cget(get_description_name(filter_ctx.desc)));
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe distance from closest boundary was %g.\n", filter_ctx.dist);
- }
- }
- /* Check if moved to vertex/edge */
- FOR_EACH(j, 0, 2) {
- if(CLAMP(hit.uv[j], 0.0005, 0.9995) != hit.uv[j])
- danger++;
- }
- if(danger) {
- logger_print(stardis->logger, LOG_WARNING,
- "Probe %s close to / on %s. "
- "If computation fails, try moving it slightly.\n",
- (filter_ctx.dist != 0 ? "moved" : "is"), (danger == 1 ? "an edge" : "a vertex"));
- }
- if(filter_ctx.probe_on_boundary) {
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe is on primitive %u, uv = (%g, %g), not moved.\n",
- hit.prim.prim_id, SPLIT2(hit.uv));
+ /* Need to check medium */
+ if(filter_ctx.outside) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Probe is outside the model.\n");
+ logger_print(stardis->logger, LOG_ERROR,
+ "Closest geometry is primitive %u, uv = (%g, %g), pos (%g, %g, %g).\n",
+ hit.prim.prim_id, SPLIT2(hit.uv), SPLIT3(filter_ctx.pos));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(filter_ctx.probe_on_boundary) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Probe is on primitive %u, uv = (%g, %g). Use -P instead of -p.\n",
+ hit.prim.prim_id, SPLIT2(hit.uv));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!filter_ctx.desc) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Could not determine the medium probe is in. "
+ "Try moving it slightly.\n");
+ logger_print(stardis->logger, LOG_ERROR,
+ "Closest geometry is primitive %u, uv = (%g, %g), distance %g.\n",
+ hit.prim.prim_id, SPLIT2(hit.uv), filter_ctx.dist);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(DESC_IS_SOLID(filter_ctx.desc)) {
+ double delta;
+ if(filter_ctx.desc->type == DESC_MAT_SOLID) {
+ struct solid* solid = filter_ctx.desc->d.solid;
+ delta = solid->delta;
} else {
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe moved to (%g, %g, %g), primitive %u, uv = (%g, %g).\n",
- SPLIT3(filter_ctx.pos), hit.prim.prim_id, SPLIT2(hit.uv));
+ struct solid_prog* solid_prog = filter_ctx.desc->d.solid_prog;
+ struct stardis_vertex vtx;
+ ASSERT(filter_ctx.desc->type == DESC_MAT_SOLID_PROG);
+ d3_set(vtx.P, pos);
+ vtx.time = time;
+ delta = solid_prog->delta(&vtx, solid_prog->prog_data);
}
- d3_set_f3(pos, filter_ctx.pos);
- } else {
- /* Need to check medium */
- if(filter_ctx.outside) {
+ if(filter_ctx.dist < 0.25 * delta) {
logger_print(stardis->logger, LOG_ERROR,
- "Probe is outside the model.\n");
- logger_print(stardis->logger, LOG_ERROR,
- "Closest geometry is primitive %u, uv = (%g, %g), pos (%g, %g, %g).\n",
- hit.prim.prim_id, SPLIT2(hit.uv), SPLIT3(filter_ctx.pos));
- res = RES_BAD_ARG;
- goto error;
- }
- if(filter_ctx.probe_on_boundary) {
+ "Probe is %g delta from closest boundary. Use -P instead of -p.\n",
+ filter_ctx.dist / delta);
logger_print(stardis->logger, LOG_ERROR,
- "Probe is on primitive %u, uv = (%g, %g). Use -P instead of -p.\n",
+ "Closest geometry is primitive %u, uv = (%g, %g).\n",
hit.prim.prim_id, SPLIT2(hit.uv));
res = RES_BAD_ARG;
goto error;
}
- if(!filter_ctx.desc) {
- logger_print(stardis->logger, LOG_ERROR,
- "Could not determine the medium probe is in. "
- "Try moving it slightly.\n");
- logger_print(stardis->logger, LOG_ERROR,
- "Closest geometry is primitive %u, uv = (%g, %g), distance %g.\n",
- hit.prim.prim_id, SPLIT2(hit.uv), filter_ctx.dist);
- res = RES_BAD_ARG;
- goto error;
+ if(filter_ctx.dist < 0.5 * delta) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe is %g delta from closest boundary. "
+ "Consider using -P instead of -p.\n",
+ filter_ctx.dist / delta);
+ } else {
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Probe is %g delta from closest boundary.\n",
+ filter_ctx.dist / delta);
}
- if(DESC_IS_SOLID(filter_ctx.desc)) {
- double delta;
- if(filter_ctx.desc->type == DESC_MAT_SOLID) {
- struct solid* solid = filter_ctx.desc->d.solid;
- delta = solid->delta;
- } else {
- struct solid_prog* solid_prog = filter_ctx.desc->d.solid_prog;
- struct stardis_vertex vtx;
- ASSERT(filter_ctx.desc->type == DESC_MAT_SOLID_PROG);
- d3_set(vtx.P, pos);
- vtx.time = time;
- delta = solid_prog->delta(&vtx, solid_prog->prog_data);
- }
- if(filter_ctx.dist < 0.25 * delta) {
- logger_print(stardis->logger, LOG_ERROR,
- "Probe is %g delta from closest boundary. Use -P instead of -p.\n",
- filter_ctx.dist / delta);
- logger_print(stardis->logger, LOG_ERROR,
- "Closest geometry is primitive %u, uv = (%g, %g).\n",
- hit.prim.prim_id, SPLIT2(hit.uv));
- res = RES_BAD_ARG;
- goto error;
- }
- if(filter_ctx.dist < 0.5 * delta) {
- logger_print(stardis->logger, LOG_WARNING,
- "Probe is %g delta from closest boundary. "
- "Consider using -P instead of -p.\n",
- filter_ctx.dist / delta);
- } else {
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe is %g delta from closest boundary.\n",
- filter_ctx.dist / delta);
- }
+ } else {
+ ASSERT(DESC_IS_FLUID(filter_ctx.desc));
+ /* In fluid; TODO: check distance wrt local geometry (use 4V/S?) */
+ if(filter_ctx.desc->type == DESC_MAT_FLUID) {
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe is in fluid '%s': computing fluid temperature, "
+ "not using a specific position.\n",
+ str_cget(&filter_ctx.desc->d.fluid->name));
} else {
- ASSERT(DESC_IS_FLUID(filter_ctx.desc));
- /* In fluid; TODO: check distance wrt local geometry (use 4V/S?) */
- if(filter_ctx.desc->type == DESC_MAT_FLUID) {
- logger_print(stardis->logger, LOG_WARNING,
- "Probe is in fluid '%s': computing fluid temperature, "
- "not using a specific position.\n",
- str_cget(&filter_ctx.desc->d.fluid->name));
- } else {
- ASSERT(filter_ctx.desc->type == DESC_MAT_FLUID_PROG);
- logger_print(stardis->logger, LOG_WARNING,
- "Probe is in fluid_prog '%s': computing fluid temperature, "
- "not using a specific position.\n",
- str_cget(&filter_ctx.desc->d.fluid_prog->name));
- }
+ ASSERT(filter_ctx.desc->type == DESC_MAT_FLUID_PROG);
+ logger_print(stardis->logger, LOG_WARNING,
+ "Probe is in fluid_prog '%s': computing fluid temperature, "
+ "not using a specific position.\n",
+ str_cget(&filter_ctx.desc->d.fluid_prog->name));
}
}
@@ -490,7 +394,7 @@ compute_probe(struct stardis* stardis, struct time* start)
ASSERT(stardis && start && (stardis->mode & MODE_PROBE_COMPUTE));
- ERR(check_probe_conform_to_type(stardis, 0, stardis->probe,
+ ERR(check_probe_conform_to_type(stardis, stardis->probe,
stardis->time_range[0], &iprim, uv));
args.nrealisations = stardis->samples;
@@ -585,268 +489,6 @@ error:
}
static res_T
-compute_probe_on_interface(struct stardis* stardis, struct time* start)
-{
- res_T res = RES_OK;
- double uv[2] = { 0,0 };
- unsigned iprim = UINT_MAX;
- struct sdis_estimator* estimator = NULL;
- struct sdis_green_function* green = NULL;
- struct dump_path_context dump_ctx;
- struct sdis_solve_probe_boundary_args args
- = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
- FILE* stream_r = NULL;
- FILE* stream_g = NULL;
- FILE* stream_p = NULL;
- struct time compute_start, compute_end;
- enum sdis_side compute_side;
- unsigned prop[3];
- const struct description* descriptions
- = darray_descriptions_cdata_get(&stardis->descriptions);
- double tcr = 0;
- const char* solve_name;
- const char* compute_side_name = NULL;
- const char* medium_name = NULL;
- const struct description* intface_desc = NULL;
-#ifndef NDEBUG
- const size_t dcount = darray_descriptions_size_get(&stardis->descriptions);
-#endif
-
- ASSERT(stardis && start && (stardis->mode & MODE_PROBE_COMPUTE_ON_INTERFACE));
- ERR(check_probe_conform_to_type(stardis, 1, stardis->probe,
- stardis->time_range[0], &iprim, uv));
- ASSERT(iprim != UINT_MAX);
-
- ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d,
- (unsigned)iprim, prop));
-
- if(prop[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY) {
- ASSERT(prop[SG3D_INTFACE] < dcount);
- intface_desc = descriptions + prop[SG3D_INTFACE];
- if(intface_desc->type == DESC_SOLID_SOLID_CONNECT)
- tcr = intface_desc->d.ss_connect->tcr;
- }
-
- if(str_is_empty(&stardis->solve_name)) {
- solve_name = NULL;
- compute_side = SDIS_SIDE_NULL__; /* Side is let unspecified */
- } else {
- solve_name = str_cget(&stardis->solve_name);
- if(0 == strcasecmp(solve_name, "FRONT")) {
- const struct description* d;
- ASSERT(prop[SG3D_FRONT] < dcount);
- d = descriptions + prop[SG3D_FRONT];
- ASSERT(DESC_IS_MEDIUM(d));
- medium_name = str_cget(get_description_name(d));
- compute_side = SDIS_FRONT;
- compute_side_name = "FRONT";
- }
- else if(0 == strcasecmp(solve_name, "BACK")) {
- const struct description* d;
- ASSERT(prop[SG3D_BACK] < dcount);
- d = descriptions + prop[SG3D_BACK];
- ASSERT(DESC_IS_MEDIUM(d));
- medium_name = str_cget(get_description_name(d));
- compute_side = SDIS_BACK;
- compute_side_name = "BACK";
- } else { /* solve_name must be a medium name */
- unsigned med_id, fmat_id, bmat_id;
- const struct description* fd;
- const struct description* bd;
- struct sdis_medium* medium
- = find_medium_by_name(stardis, &stardis->solve_name, &med_id);
- if(medium == NULL) { /* Not found */
- logger_print(stardis->logger, LOG_ERROR,
- "Cannot locate side from medium name '%s' (unknown medium)\n",
- solve_name);
- res = RES_BAD_ARG;
- goto error;
- }
-
- fd = descriptions + prop[SG3D_FRONT];
- bd = descriptions + prop[SG3D_BACK];
- ASSERT(DESC_IS_MEDIUM(fd));
- ASSERT(DESC_IS_MEDIUM(bd));
- description_get_medium_id(fd, &fmat_id);
- description_get_medium_id(bd, &bmat_id);
-
- if(med_id != fmat_id && med_id != bmat_id) { /* Not here */
- logger_print(stardis->logger, LOG_ERROR,
- "Medium '%s' is not used at this interface (prim id=%u)\n",
- solve_name, iprim);
- res = RES_BAD_ARG;
- goto error;
- }
- if(fmat_id == bmat_id) { /* Cannot differentiate */
- unsigned encs[2];
- ERR(senc3d_scene_get_triangle_enclosures(stardis->senc3d_scn, iprim,
- encs));
- logger_print(stardis->logger, LOG_ERROR,
- "Medium '%s' is used on both sides of this interface (prim id=%u)\n",
- solve_name,iprim);
- logger_print(stardis->logger, LOG_ERROR,
- "Side must be defined using either FRONT or BACK.\n");
- logger_print(stardis->logger, LOG_ERROR,
- "FRONT side is related to enclosure %u, BACK side to enclosure %u.\n",
- encs[SENC3D_FRONT], encs[SENC3D_BACK]);
- res = RES_BAD_ARG;
- goto error;
- }
- medium_name = solve_name;
- if(med_id == fmat_id) {
- compute_side = SDIS_FRONT;
- compute_side_name = "FRONT";
- } else {
- compute_side = SDIS_BACK;
- compute_side_name = "BACK";
- }
- }
- }
-
- if(compute_side == SDIS_SIDE_NULL__) {
- /* Cannot have defined a contact resistance here */
- if(tcr != 0) {
- const struct description* fdesc = NULL;
- const struct description* bdesc = NULL;
- if(prop[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) {
- ASSERT(prop[SG3D_FRONT] < dcount);
- fdesc = descriptions + prop[SG3D_FRONT];
- }
- if(prop[SG3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) {
- ASSERT(prop[SG3D_BACK] < dcount);
- bdesc = descriptions + prop[SG3D_BACK];
- }
- ASSERT(fdesc || bdesc);
- logger_print(stardis->logger, LOG_ERROR,
- "Cannot let probe computation side unspecified on an interface with a "
- "non-nul thermal resistance.\n");
- logger_print(stardis->logger, LOG_ERROR,
- "Interface '%s',%s%s%s%s%s%s%s, resistance=%g K.m^2/W.\n",
- str_cget(get_description_name(intface_desc)),
- (fdesc ? " FRONT:'" : ""),
- (fdesc ? str_cget(get_description_name(fdesc)) : ""),
- (fdesc ? "'" : ""),
- (fdesc && bdesc ? "," : ""),
- (bdesc ? " BACK:'" : ""),
- (bdesc ? str_cget(get_description_name(bdesc)) : ""),
- (bdesc ? "'" : ""),
- tcr);
- res = RES_BAD_ARG;
- goto error;
- }
- /* Any side is OK as temperature is the same */
- compute_side = SDIS_FRONT;
- } else if (tcr == 0) {
- /* Warn if side defined and no resistance defined */
- logger_print(stardis->logger, LOG_WARNING,
- "Specifying a compute side at an interface with no contact resistance "
- "is meaningless.\n");
- }
-
- if(tcr > 0) {
- ASSERT(compute_side_name && medium_name);
- logger_print(stardis->logger, LOG_OUTPUT,
- "Probe computation on an interface with a thermal resistance = %g K.m^2/W "
- "on %s side (medium is '%s').\n",
- tcr, compute_side_name, medium_name);
- }
-
- args.nrealisations = stardis->samples;
- args.iprim = iprim;
- d2_set(args.uv, uv);
- args.side = compute_side;
- d2_set(args.time_range, stardis->time_range);
-
- /* Input random state? */
- READ_RANDOM_STATE(&stardis->rndgen_state_in_filename);
-
- if(stardis->mode & (MODE_BIN_GREEN | MODE_GREEN)) {
- if(stardis->mpi_initialized && stardis->mpi_rank != 0) {
- ERR(sdis_solve_probe_boundary_green_function(stardis->sdis_scn, &args,
- &green));
- } else {
- /* Try to open output files to detect errors early */
- if(stardis->mode & MODE_BIN_GREEN) {
- ASSERT(!str_is_empty(&stardis->bin_green_filename));
- stream_g = fopen(str_cget(&stardis->bin_green_filename), "wb");
- if(!stream_g) {
- logger_print(stardis->logger, LOG_ERROR,
- "cannot open file '%s' for binary writing.\n",
- str_cget(&stardis->bin_green_filename));
- res = RES_IO_ERR;
- goto error;
- }
- if(str_cget(&stardis->end_paths_filename)
- && strlen(str_cget(&stardis->end_paths_filename)))
- {
- stream_p = fopen(str_cget(&stardis->end_paths_filename), "w");
- if(!stream_p) {
- logger_print(stardis->logger, LOG_ERROR,
- "cannot open file '%s' for writing.\n",
- str_cget(&stardis->end_paths_filename));
- res = RES_IO_ERR;
- goto error;
- }
- }
- }
- /* Call solve() */
- time_current(&compute_start);
- ERR(sdis_solve_probe_boundary_green_function(stardis->sdis_scn, &args,
- &green));
- time_current(&compute_end);
- if(stardis->mode & MODE_BIN_GREEN) {
- struct time output_end;
- ERR(dump_green_bin(green, stardis, stream_g));
- if(stream_p) {
- ERR(dump_paths_end(green, stardis, stream_p));
- }
- time_current(&output_end);
- ERR(print_computation_time(NULL, stardis,
- start, &compute_start, &compute_end, &output_end));
- }
- if(stardis->mode & MODE_GREEN) {
- ERR(dump_green_ascii(green, stardis, stdout));
- }
- }
- } else {
- args.register_paths = stardis->dump_paths;
- args.picard_order = stardis->picard_order;
- if(stardis->mpi_initialized && stardis->mpi_rank != 0) {
- ERR(sdis_solve_probe_boundary(stardis->sdis_scn, &args, &estimator));
- } else {
- res_T tmp_res1, tmp_res2;
- time_current(&compute_start);
- ERR(sdis_solve_probe_boundary(stardis->sdis_scn, &args, &estimator));
- time_current(&compute_end);
- ERR(print_computation_time(estimator, stardis,
- start, &compute_start, &compute_end, NULL));
- tmp_res1 = print_single_MC_result(estimator, stardis, stdout);
-
- /* Dump recorded paths according to user settings */
- dump_ctx.stardis = stardis;
- dump_ctx.rank = 0;
- tmp_res2 = sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx);
- if(tmp_res1 != RES_OK) res = tmp_res1;
- else if(tmp_res2 != RES_OK) res = tmp_res2;
- }
- }
-
- /* Output random state? */
- WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename);
-
-end:
- if(stream_r) fclose(stream_r);
- if(stream_g) fclose(stream_g);
- if(stream_p) fclose(stream_p);
- if(estimator) SDIS(estimator_ref_put(estimator));
- if(green) SDIS(green_function_ref_put(green));
- if(args.rng_state) SSP(rng_ref_put(args.rng_state));
- return res;
-error:
- goto end;
-}
-
-static res_T
auto_look_at
(struct sdis_scene* scn,
const double fov_x, /* Horizontal field of view in radian */
@@ -1040,7 +682,7 @@ compute_medium(struct stardis* stardis, struct time* start)
ASSERT(stardis && start && (stardis->mode & MODE_MEDIUM_COMPUTE));
/* Find medium */
- medium = find_medium_by_name(stardis, &stardis->solve_name, NULL);
+ medium = find_medium_by_name(stardis, str_cget(&stardis->solve_name), NULL);
if(medium == NULL) { /* Not found */
logger_print(stardis->logger, LOG_ERROR,
"Cannot solve medium '%s' (unknown medium)\n",
@@ -1418,7 +1060,7 @@ error:
struct sdis_medium*
find_medium_by_name
(struct stardis* stardis,
- const struct str* name,
+ const char* name,
unsigned* out_id)
{
struct sdis_medium* medium = NULL;
@@ -1430,7 +1072,7 @@ find_medium_by_name
unsigned id;
struct description* desc =
darray_descriptions_data_get(&stardis->descriptions) + i;
- if(!str_eq(name, get_description_name(desc)))
+ if(strcmp(name, str_cget(get_description_name(desc))) != 0)
continue;
description_get_medium_id(desc, &id);
ASSERT(darray_media_ptr_size_get(&stardis->media) > id);
diff --git a/src/stardis-compute.h b/src/stardis-compute.h
@@ -27,19 +27,24 @@ struct str;
#define DARRAY_DATA struct sdis_estimator*
#include <rsys/dynamic_array.h>
-struct sdis_medium*
+extern LOCAL_SYM struct sdis_medium*
find_medium_by_name
(struct stardis* stardis,
- const struct str* name,
+ const char* name,
unsigned* id); /* Can be NULL */
-res_T
+extern LOCAL_SYM res_T
stardis_compute
(struct stardis* stardis,
struct time* start);
-res_T
+extern LOCAL_SYM res_T
read_compute_surface
(struct stardis* stardis);
+extern LOCAL_SYM res_T
+compute_probe_on_interface
+ (struct stardis* stardis,
+ struct time* start);
+
#endif
diff --git a/src/stardis-output.c b/src/stardis-output.c
@@ -1956,7 +1956,8 @@ dump_compute_region_at_the_end_of_vtk
const struct description* descriptions;
unsigned medium_id;
ASSERT(stardis->mode & MODE_MEDIUM_COMPUTE);
- medium = find_medium_by_name(stardis, &stardis->solve_name, &medium_id);
+ medium = find_medium_by_name
+ (stardis, str_cget(&stardis->solve_name), &medium_id);
ASSERT(medium != NULL); (void)medium;
descriptions = darray_descriptions_cdata_get(&stardis->descriptions);
FOR_EACH(i, 0, tsz) {
diff --git a/src/stardis-output.h b/src/stardis-output.h
@@ -38,6 +38,9 @@ struct dump_path_context {
unsigned long rank;
struct stardis* stardis;
};
+#define DUMP_PATH_CONTEXT_NULL__ {0,NULL}
+static const struct dump_path_context DUMP_PATH_CONTEXT_NULL =
+ DUMP_PATH_CONTEXT_NULL__;
res_T
dump_path
diff --git a/src/stardis-prog-properties.h.in b/src/stardis-prog-properties.h.in
@@ -36,6 +36,8 @@ struct stardis_vertex {
double P[3]; /* World space position */
double time; /* "Time" of the vertex */
};
+#define STARDIS_VERTEX_NULL__ {{0,0,0}, 0}
+static const struct stardis_vertex STARDIS_VERTEX_NULL = STARDIS_VERTEX_NULL__;
enum stardis_side {
FRONT,
