commit 3c8d91d5b8021a655bcff307a1aae645d5248a78
parent 21fae3ba4f4bde64a1554c1cecb4d909e596ba4c
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Mon, 24 Feb 2025 17:19:32 +0100
Add the -f option to compute flux density probes
Diffstat:
9 files changed, 320 insertions(+), 81 deletions(-)
diff --git a/doc/stardis-output.5 b/doc/stardis-output.5
@@ -112,6 +112,7 @@ or an extended format that mixes numbers and their descriptions
.It Ta \& \& | Ta Ao Va medium-temp Ac # Option Fl m
.It Ta \& \& | Ta Ao Va mean-temp Ac # Option Fl s
.It Ta \& \& | Ta Ao Va mean-flux Ac # Option Fl F
+.It Ta \& \& | Ta Ao Va flux-density Ac # Option Fl f
.It \ Ta Ta
.\" Probe temperature
.It Ao Va probe-temp Ac Ta ::= Ta Ao Va probe-temp-raw Ac | Ao Va probe-temp-ext Ac
@@ -133,8 +134,11 @@ or an extended format that mixes numbers and their descriptions
.It \ Ta Ta
.\" Mean flux
.It Ao Va mean-flux Ac Ta ::= Ta Ao Va mean-flux-raw Ac | Ao Va mean-flux-ext Ac
+.It Ao Va flux-density Ac Ta ::= Ta Ao Va flux-density-raw Ac | Ao Va flux-density-ext Ac
.It Ao Va mean-flux-raw Ac Ta ::= Ta Ao Va estimate Ac Ao Va estimate Ac Ao Va estimate Ac \e
.It Ta Ta Ao Va estimate Ac Ao Va estimate Ac Ao Va failures Ac
+.It Ao Va flux-density-raw Ac Ta ::= Ta Ao Va estimate Ac Ao Va estimate Ac Ao Va estimate Ac \e
+.It Ta Ta Ao Va estimate Ac Ao Va estimate Ac Ao Va failures Ac
.It Ao Va mean-flux-ext Ac Ta ::= Ta Li Temperature at boundary Ao Va stl-path Ac \e
.It Ta Ta Ao Va time Ac Ao Va estimate-temp-ext Ac
.It Ta Ta Li Convective flux at boundary Ao Va stl-path Ac \e
@@ -146,11 +150,23 @@ or an extended format that mixes numbers and their descriptions
.It Ta Ta Li Total flux at boundary Ao Va stl-path Ac \e
.It Ta Ta Ao Va time Ac Ao Va estimate-flux-ext Ac
.It Ta Ta Ao Va failures-ext Ac
+.It Ao Va flux-density-ext Ac Ta ::= Ta Li Temperature at Ao Va position Ac Ao Va time Ac \e
+.It Ta Ta Ao Va estimate-temp-ext Ac
+.It Ta Ta Li Convective flux density at Ao Va position Ac Ao Va time Ac \e
+.It Ta Ta Ao Va estimate-f-dsity-ext Ac
+.It Ta Ta Li Radiative flux density at Ao Va position Ac Ao Va time Ac \e
+.It Ta Ta Ao Va estimate-f-dsity-ext Ac
+.It Ta Ta Li Imposed flux density at Ao Va position Ac Ao Va time Ac \e
+.It Ta Ta Ao Va estimate-f-dsity-ext Ac
+.It Ta Ta Li Total flux density at Ao Va position Ac Ao Va time Ac \e
+.It Ta Ta Ao Va estimate-f-dsity-ext Ac
+.It Ta Ta Ao Va failures-ext Ac
.It \ Ta Ta
.\" Miscellaneous
.It Ao Va estimate Ac Ta ::= Ta Ao Va expected-value Ac Ao Va standard-error Ac
.It Ao Va estimate-temp-ext Ac Ta ::= Ta Ao Va expected-value Ac Li K +/- Ao Va standard-error Ac
.It Ao Va estimate-flux-ext Ac Ta ::= Ta Ao Va expected-value Ac Li W +/- Ao Va standard-error Ac
+.It Ao Va estimate-f-dsity-ext Ac Ta ::= Ta Ao Va expected-value Ac Li W/m² +/- Ao Va standard-error Ac
.It Ao Va expected-value Ac Ta ::= Ta Vt real
.It Ao Va standard-error Ac Ta ::= Ta Vt real
.It \ Ta Ta
diff --git a/doc/stardis.1.in b/doc/stardis.1.in
@@ -27,6 +27,7 @@
.Op Fl D Ar path_type , Ns Ar files_name_prefix
.Op Fl d Ar file_base_name
.Op Fl F Pa surface Ns Op , Ns Ar time Ns Op , Ns Ar time
+.Op Fl f Ar x , Ns Ar y , Ns Ar z Ns Op , Ns Ar time Ns Op , Ns Ar time
.Op Fl G Pa green_bin Ns Op , Ns Pa green_ascii
.Op Fl I Ar initial_time
.Op Fl L Pa interface_probes
@@ -209,6 +210,14 @@ positive from the solid to the fluid.
These triangles are not added to the geometry, but must be part of it.
By default the compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@.
The surface does not need to be connex.
+.It Fl f Ar x , Ns Ar y , Ns Ar z Ns Op , Ns Ar time Ns Op , Ns Ar time
+Compute the flux density at a given probe on an interface at a given time.
+The probe is supposed to be on an interface and is moved to the closest point of
+the closest interface before the computation starts.
+The flux density can only be computed on a solid-fluid interface and is
+accounted positive from the solid to the fluid.
+The probe coordinates must be in the same system as the geometry.
+By default the compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@.
.It Fl G Pa green_bin Ns Op , Ns Pa green_ascii
Compute the Green function at the specified time and write it to a binary file.
If a
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -328,6 +328,7 @@ stardis_init
ERR(str_set(&stardis->solve_name, args->medium_name));
}
else if((args->mode & MODE_PROBE_COMPUTE_ON_INTERFACE)
+ || (args->mode & MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE)
|| (args->mode & MODE_PROBE_LIST_COMPUTE_ON_INTERFACE)) {
ERR(darray_probe_boundary_copy
(&stardis->probe_boundary_list, &args->probe_boundary_list));
diff --git a/src/stardis-args.c b/src/stardis-args.c
@@ -113,6 +113,7 @@ mode_option(const int m)
if(m & MODE_DUMP_PATHS) { found++; res = 'D'; }
if(m & MODE_EXTENDED_RESULTS) { found++; res = 'e'; }
if(m & MODE_FLUX_BOUNDARY_COMPUTE) { found++; res = 'F'; }
+ if(m & MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE) { found++; res = 'f'; }
if(m & MODE_GREEN) { found++; res = 'g'; }
if(m & MODE_BIN_GREEN) { found++; res = 'G'; }
if(m & MODE_DUMP_HELP) { found++; res = 'h'; }
@@ -459,10 +460,10 @@ usage(FILE* stream)
#define PRINT(Msg) fprintf(stream, Msg)
PRINT("stardis [-eghiv] [-a diff_algo] [-D path_type,files_name_prefix]\n");
PRINT(" [-d file_base_name] [-F surface[,time[,time]]]\n");
- PRINT(" [-G green_bin[,green_ascii]] [-I initial_time]\n");
- PRINT(" [-L interface_probes] [-m medium_name[,time[,time]]]\n");
- PRINT(" [-n samples_count] [-o picard_order]\n");
- PRINT(" [-P x,y,z[,time[,time]][:side_indicator]]\n");
+ PRINT(" [-f x,y,z[,time[,time]]] [-G green_bin[,green_ascii]]\n");
+ PRINT(" [-I initial_time] [-L interface_probes]\n");
+ PRINT(" [-m medium_name[,time[,time]]] [-n samples_count]\n");
+ PRINT(" [-o picard_order] [-P x,y,z[,time[,time]][:side_indicator]]\n");
PRINT(" [-p x,y,z[,time[,time]]] [-R rendering_opt[:rendering_opt ...]]\n");
PRINT(" [-S surface[,time[,time]]] [-s surface[,time[,time]]]\n");
PRINT(" [-t threads_count] [-V verbosity_level] [-X output_rng]\n");
@@ -525,7 +526,7 @@ parse_args
{
int opt = 0, n_used = 0, o_used = 0;
size_t len = 0;
- const char option_list[] = "a:c:d:D:eF:gG:hiI:L:m:M:n:o:p:P:R:s:S:t:vV:x:X:";
+ const char option_list[] = "a:c:d:D:ef:F:gG:hiI:L:m:M:n:o:p:P:R:s:S:t:vV:x:X:";
char buf[128];
struct str keep;
char** line = NULL;
@@ -609,6 +610,25 @@ parse_args
args->mode |= MODE_EXTENDED_RESULTS;
break;
+ case 'f':
+ {
+ struct stardis_probe_boundary* probe = NULL;
+
+ if(args->mode & EXCLUSIVE_MODES) {
+ logger_print(args->logger, LOG_ERROR,
+ "Options -%c and -%c are exclusive.\n",
+ (char)opt, mode_option(args->mode));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ args->mode |= MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE;
+ res = allocate_probe_boundary(args, &probe);
+ if(res != RES_OK) goto error;
+ res = parse_position_and_time(optarg, probe->position, probe->time);
+ if(res != RES_OK) goto error;
+ break;
+ }
+
/*case 'F': see 's' */
case 'g':
diff --git a/src/stardis-args.h.in b/src/stardis-args.h.in
@@ -50,18 +50,19 @@ enum stardis_mode {
MODE_DUMP_MODEL = BIT(2), /* -d */
MODE_EXTENDED_RESULTS = BIT(3), /* -e */
MODE_FLUX_BOUNDARY_COMPUTE = BIT(4), /* -F */
- MODE_BIN_GREEN = BIT(5), /* -G */
- MODE_GREEN = BIT(6), /* -g */
- MODE_DUMP_HELP = BIT(7), /* -h */
- MODE_PROBE_LIST_COMPUTE_ON_INTERFACE = BIT(8), /* -L */
- MODE_MEDIUM_COMPUTE = BIT(9), /* -m */
- MODE_PROBE_COMPUTE_ON_INTERFACE = BIT(10), /* -P */
- MODE_PROBE_COMPUTE = BIT(11), /* -p */
- MODE_IR_COMPUTE = BIT(12), /* -R */
- MODE_MAP_COMPUTE = BIT(13), /* -S */
- MODE_BOUNDARY_COMPUTE = BIT(14), /* -s */
- MODE_VERBOSITY = BIT(15), /* -V */
- MODE_DUMP_VERSION = BIT(16), /* -v */
+ MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE = BIT(5), /* -f */
+ MODE_BIN_GREEN = BIT(6), /* -G */
+ MODE_GREEN = BIT(7), /* -g */
+ MODE_DUMP_HELP = BIT(8), /* -h */
+ MODE_PROBE_LIST_COMPUTE_ON_INTERFACE = BIT(9), /* -L */
+ MODE_MEDIUM_COMPUTE = BIT(10), /* -m */
+ MODE_PROBE_COMPUTE_ON_INTERFACE = BIT(11), /* -P */
+ MODE_PROBE_COMPUTE = BIT(12), /* -p */
+ MODE_IR_COMPUTE = BIT(13), /* -R */
+ MODE_MAP_COMPUTE = BIT(14), /* -S */
+ MODE_BOUNDARY_COMPUTE = BIT(15), /* -s */
+ MODE_VERBOSITY = BIT(16), /* -V */
+ MODE_DUMP_VERSION = BIT(17), /* -v */
GREEN_COMPATIBLE_MODES =
MODE_PROBE_COMPUTE
@@ -77,7 +78,8 @@ enum stardis_mode {
EXT_COMPATIBLE_MODES =
GREEN_COMPATIBLE_MODES
| MODE_MEDIUM_COMPUTE
- | MODE_FLUX_BOUNDARY_COMPUTE,
+ | MODE_FLUX_BOUNDARY_COMPUTE
+ | MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE,
REGION_COMPUTE_MODES =
SURFACE_COMPUTE_MODES
@@ -87,7 +89,8 @@ enum stardis_mode {
GREEN_COMPATIBLE_MODES
| MODE_IR_COMPUTE
| SURFACE_COMPUTE_MODES
- | MODE_PROBE_LIST_COMPUTE_ON_INTERFACE,
+ | MODE_PROBE_LIST_COMPUTE_ON_INTERFACE
+ | MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE,
EXCLUSIVE_MODES = COMPUTE_MODES,
diff --git a/src/stardis-compute-probe-boundary.c b/src/stardis-compute-probe-boundary.c
@@ -235,7 +235,8 @@ check_move_to_solid_boundary
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 */
+ const double distance, /* Move distance */
+ const int advice)
{
struct stardis_vertex vtx = STARDIS_VERTEX_NULL;
const char* prefix = "";
@@ -261,7 +262,7 @@ check_move_to_solid_boundary
vtx.P[2] = pos[2];
vtx.time = time;
delta = desc->d.solid_prog->delta(&vtx, desc->d.solid_prog->prog_data);
- prefix = "programmed";
+ prefix = "programmed ";
break;
default: FATAL("Unreachable code.\n");
@@ -271,7 +272,7 @@ check_move_to_solid_boundary
logger_print(stardis->logger, LOG_OUTPUT, "Probe was in %ssolid '%s'.\n",
prefix, solid_name);
- /* The position is closed from the triangle */
+ /* The position is close from the triangle */
if(distance < 0.5*delta) {
logger_print(stardis->logger, LOG_OUTPUT,
"Probe was %g delta from closest boundary.\n",
@@ -280,9 +281,9 @@ check_move_to_solid_boundary
/* 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);
+ "Probe was %g delta from closest boundary.%s\n",
+ distance/delta,
+ (advice ? " Consider using -p instead of -P.\n" : ""));
/* The position is too far from the triangle */
} else {
@@ -316,7 +317,7 @@ check_move_to_fluid_boundary
switch(desc->type) {
case DESC_MAT_FLUID: prefix = ""; break;
- case DESC_MAT_FLUID_PROG: prefix = "programmed"; break;
+ case DESC_MAT_FLUID_PROG: prefix = "programmed "; break;
default: FATAL("Unreachable code.\n");
}
@@ -334,6 +335,7 @@ move_to_boundary
(struct stardis* stardis,
const double pos[3],
const double time, /* [s] */
+ const int is_probe_temp_computation,
size_t* out_iprim,
double uv[2])
{
@@ -342,7 +344,6 @@ move_to_boundary
double proj_pos[3] = {0,0,0};
size_t iprim = 0;
-
/* Miscellaneous */
size_t nvertices_close = 0;
res_T res = RES_OK;
@@ -376,7 +377,8 @@ move_to_boundary
case DESC_MAT_SOLID:
case DESC_MAT_SOLID_PROG:
ERR(check_move_to_solid_boundary
- (stardis, pos, time, desc, iprim, uv, filter_ctx.distance));
+ (stardis, pos, time, desc, iprim, uv, filter_ctx.distance,
+ is_probe_temp_computation));
break;
case DESC_MAT_FLUID:
case DESC_MAT_FLUID_PROG:
@@ -826,6 +828,91 @@ error:
}
static res_T
+setup_solve_probe_boundary_flux_args
+ (struct stardis* stardis,
+ const struct stardis_probe_boundary* probe,
+ struct sdis_solve_probe_boundary_flux_args* args)
+{
+ double uv[2] = {0, 0};
+ size_t iprim = SIZE_MAX;
+ unsigned desc_ids[SG3D_PROP_TYPES_COUNT__];
+ res_T res = RES_OK;
+ ASSERT(stardis && probe && args);
+
+ /* Calculate the probe position on the boundary */
+ ERR(move_to_boundary(stardis, probe->position, probe->time[0], 0, &iprim, uv));
+
+ ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d,
+ (unsigned)iprim, desc_ids));
+
+ d3_set(stardis->probe, probe->position);
+ /* Setup of solve input parameters */
+
+ args->nrealisations = stardis->samples;
+ args->iprim = iprim;
+ d2_set(args->uv, uv);
+ args->picard_order = stardis->picard_order;
+ d2_set(args->time_range, stardis->time_range);
+ args->diff_algo = stardis->diff_algo;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+compute_single_flux_probe_on_interface
+ (struct stardis* stardis,
+ struct time* start,
+ const struct stardis_probe_boundary* probe)
+{
+ res_T res = RES_OK;
+ struct sdis_green_function* green = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_solve_probe_boundary_flux_args args
+ = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT;
+ struct dump_path_context dump_ctx;
+ FILE* stream_r = NULL;
+ struct time compute_start, compute_end;
+
+ ASSERT(stardis && start && probe);
+
+ ERR(setup_solve_probe_boundary_flux_args(stardis, probe, &args));
+
+ /* Input random state? */
+ READ_RANDOM_STATE(&stardis->rndgen_state_in_filename);
+
+ if(stardis->mpi_initialized && stardis->mpi_rank != 0) {
+ ERR(sdis_solve_probe_boundary_flux(stardis->sdis_scn, &args, &estimator));
+ } else {
+ struct sdis_mc time = SDIS_MC_NULL;
+ time_current(&compute_start);
+ ERR(sdis_solve_probe_boundary_flux(stardis->sdis_scn, &args, &estimator));
+ time_current(&compute_end);
+ ERR(sdis_estimator_get_realisation_time(estimator, &time));
+ ERR(print_computation_time
+ (&time, stardis, start, &compute_start, &compute_end, NULL));
+
+ ERR(print_single_MC_result(estimator, stardis, stdout));
+
+ /* TODO: warning if dump paths requested */
+
+ }
+
+ /* Output random state? */
+ WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename);
+
+end:
+ 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
setup_solve_probe_boundary_args
(struct stardis* stardis,
const struct stardis_probe_boundary* probe,
@@ -836,10 +923,10 @@ setup_solve_probe_boundary_args
size_t iprim = SIZE_MAX;
unsigned desc_ids[SG3D_PROP_TYPES_COUNT__];
res_T res = RES_OK;
- ASSERT(stardis && probe && args);
+ ASSERT(stardis && probe && args && (stardis->mode && MODE_MAP_COMPUTE));
/* Calculate the probe position on the boundary */
- ERR(move_to_boundary(stardis, probe->position, probe->time[0], &iprim, uv));
+ ERR(move_to_boundary(stardis, probe->position, probe->time[0], 1, &iprim, uv));
ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d,
(unsigned)iprim, desc_ids));
@@ -962,7 +1049,8 @@ compute_probe_on_interface(struct stardis* stardis, struct time* start)
res_T res = RES_OK;
ASSERT(stardis && start);
ASSERT((stardis->mode & MODE_PROBE_COMPUTE_ON_INTERFACE)
- || (stardis->mode & MODE_PROBE_LIST_COMPUTE_ON_INTERFACE));
+ || (stardis->mode & MODE_PROBE_LIST_COMPUTE_ON_INTERFACE)
+ || (stardis->mode & MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE));
/* Multiple probes */
if(stardis->mode & MODE_PROBE_LIST_COMPUTE_ON_INTERFACE) {
@@ -979,6 +1067,13 @@ compute_probe_on_interface(struct stardis* stardis, struct time* start)
probe = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list);
res = compute_single_probe_on_interface(stardis, start, probe);
if(res != RES_OK) goto error;
+ } else if(stardis->mode & MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE) {
+ const struct stardis_probe_boundary* probe = NULL;
+ ASSERT(darray_probe_boundary_size_get(&stardis->probe_boundary_list) == 1);
+
+ probe = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list);
+ res = compute_single_flux_probe_on_interface(stardis, start, probe);
+ if(res != RES_OK) goto error;
}
exit:
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -333,50 +333,6 @@ error:
goto end;
}
-#define READ_RANDOM_STATE(Name) \
- if(!stardis->mpi_initialized || stardis->mpi_rank == 0) { \
- if(str_is_empty(Name)) { \
- /* Force using threefry independently of the default RNG type */ \
- args.rng_type = SSP_RNG_THREEFRY; \
- } else { \
- const char* name = str_cget(Name); \
- stream_r = fopen(name, "r"); \
- if(!stream_r) { \
- res = RES_IO_ERR; \
- logger_print(stardis->logger, LOG_ERROR, \
- "Cannot open generator's state file ('%s').\n", \
- name); \
- goto error; \
- } \
- ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &args.rng_state)); \
- res = read_random_generator_state(args.rng_state, stream_r); \
- if(res != RES_OK) { \
- logger_print(stardis->logger, LOG_ERROR, \
- "Cannot read random generator's state ('%s').\n", \
- name); \
- goto error; \
- } \
- fclose(stream_r); stream_r = NULL; \
- } \
- }
-
-#define WRITE_RANDOM_STATE(Name) \
- if(!stardis->mpi_initialized || stardis->mpi_rank == 0) { \
- if(!str_is_empty(Name)) { \
- const char* name = str_cget(Name); \
- stream_r = fopen(name, "wb"); \
- if(!stream_r) { \
- res = RES_IO_ERR; \
- logger_print(stardis->logger, LOG_ERROR, \
- "Cannot write random generator's state ('%s').\n", \
- name); \
- goto error; \
- } \
- ERR(write_random_generator_state(estimator, stream_r)); \
- fclose(stream_r); stream_r = NULL; \
- } \
- }
-
static res_T
compute_probe(struct stardis* stardis, struct time* start)
{
@@ -1065,9 +1021,6 @@ error:
goto end;
}
-#undef READ_RANDOM_STATE
-#undef WRITE_RANDOM_STATE
-
/*******************************************************************************
* Public Functions
******************************************************************************/
@@ -1188,6 +1141,8 @@ stardis_compute
ERR(compute_probe_on_interface(stardis, start));
else if(stardis->mode & MODE_PROBE_LIST_COMPUTE_ON_INTERFACE)
ERR(compute_probe_on_interface(stardis, start));
+ else if(stardis->mode & MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE)
+ ERR(compute_probe_on_interface(stardis, start));
else if(stardis->mode & MODE_IR_COMPUTE)
ERR(compute_camera(stardis, start));
else if(stardis->mode & MODE_MEDIUM_COMPUTE)
diff --git a/src/stardis-compute.h b/src/stardis-compute.h
@@ -27,6 +27,50 @@ struct str;
#define DARRAY_DATA struct sdis_estimator*
#include <rsys/dynamic_array.h>
+#define READ_RANDOM_STATE(Name) \
+ if(!stardis->mpi_initialized || stardis->mpi_rank == 0) { \
+ if(str_is_empty(Name)) { \
+ /* Force using threefry independently of the default RNG type */ \
+ args.rng_type = SSP_RNG_THREEFRY; \
+ } else { \
+ const char* name = str_cget(Name); \
+ stream_r = fopen(name, "r"); \
+ if(!stream_r) { \
+ res = RES_IO_ERR; \
+ logger_print(stardis->logger, LOG_ERROR, \
+ "Cannot open generator's state file ('%s').\n", \
+ name); \
+ goto error; \
+ } \
+ ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &args.rng_state)); \
+ res = read_random_generator_state(args.rng_state, stream_r); \
+ if(res != RES_OK) { \
+ logger_print(stardis->logger, LOG_ERROR, \
+ "Cannot read random generator's state ('%s').\n", \
+ name); \
+ goto error; \
+ } \
+ fclose(stream_r); stream_r = NULL; \
+ } \
+ }
+
+#define WRITE_RANDOM_STATE(Name) \
+ if(!stardis->mpi_initialized || stardis->mpi_rank == 0) { \
+ if(!str_is_empty(Name)) { \
+ const char* name = str_cget(Name); \
+ stream_r = fopen(name, "wb"); \
+ if(!stream_r) { \
+ res = RES_IO_ERR; \
+ logger_print(stardis->logger, LOG_ERROR, \
+ "Cannot write random generator's state ('%s').\n", \
+ name); \
+ goto error; \
+ } \
+ ERR(write_random_generator_state(estimator, stream_r)); \
+ fclose(stream_r); stream_r = NULL; \
+ } \
+ }
+
extern LOCAL_SYM struct sdis_medium*
find_medium_by_name
(struct stardis* stardis,
@@ -47,4 +91,9 @@ compute_probe_on_interface
(struct stardis* stardis,
struct time* start);
+extern LOCAL_SYM res_T
+compute_flux_density_boundary
+ (struct stardis* stardis,
+ struct time* start);
+
#endif
diff --git a/src/stardis-output.c b/src/stardis-output.c
@@ -1690,7 +1690,98 @@ print_single_MC_result
else fprintf(stream, "%g %g %lu %lu\n",
result.E, result.SE, nfailures, nsamples);
break;
- case MODE_FLUX_BOUNDARY_COMPUTE: {
+
+ case MODE_PROBE_COMPUTE_FLUX_ON_INTERFACE:
+ {
+ enum sdis_estimator_type type;
+ ERR(sdis_estimator_get_type(estimator, &type));
+ ASSERT(type == SDIS_ESTIMATOR_FLUX);
+
+ if(stardis->mode & MODE_EXTENDED_RESULTS) {
+ if(stardis->time_range[0] == stardis->time_range[1])
+ fprintf(stream, "Temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n",
+ SPLIT3(stardis->probe), stardis->time_range[0],
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ else
+ fprintf(stream,
+ "Temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n",
+ SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ ERR(sdis_estimator_get_convective_flux(estimator, &result));
+ if(stardis->time_range[0] == stardis->time_range[1])
+ fprintf(stream, "Convective flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), stardis->time_range[0],
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ else
+ fprintf(stream,
+ "Convective flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ ERR(sdis_estimator_get_radiative_flux(estimator, &result));
+ if(stardis->time_range[0] == stardis->time_range[1])
+ fprintf(stream, "Radiative flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), stardis->time_range[0],
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ else
+ fprintf(stream,
+ "Radiative flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ ERR(sdis_estimator_get_imposed_flux(estimator, &result));
+ if(stardis->time_range[0] == stardis->time_range[1])
+ fprintf(stream, "Imposed flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), stardis->time_range[0],
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ else
+ fprintf(stream,
+ "Imposed flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ ERR(sdis_estimator_get_total_flux(estimator, &result));
+ if(stardis->time_range[0] == stardis->time_range[1])
+ fprintf(stream, "Total flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), stardis->time_range[0],
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ else
+ fprintf(stream,
+ "Total flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n",
+ SPLIT3(stardis->probe), SPLIT2(stardis->time_range),
+ result.E, /* Expected value */
+ result.SE); /* Standard error */
+ } else {
+ fprintf(stream, "%g %g ", result.E, result.SE); /* Temperature */
+ ERR(sdis_estimator_get_convective_flux(estimator, &result));
+ fprintf(stream, "%g %g ",
+ result.E,
+ result.SE);
+ ERR(sdis_estimator_get_radiative_flux(estimator, &result));
+ fprintf(stream, "%g %g ",
+ result.E,
+ result.SE);
+ ERR(sdis_estimator_get_imposed_flux(estimator, &result));
+ fprintf(stream, "%g %g ",
+ result.E,
+ result.SE);
+ ERR(sdis_estimator_get_total_flux(estimator, &result));
+ fprintf(stream, "%g %g ",
+ result.E,
+ result.SE);
+ fprintf(stream, "%lu %lu\n", nfailures, nsamples);
+ }
+ break;
+ }
+
+ case MODE_FLUX_BOUNDARY_COMPUTE:
+ {
enum sdis_estimator_type type;
ERR(sdis_estimator_get_type(estimator, &type));
ASSERT(type == SDIS_ESTIMATOR_FLUX);
@@ -1745,13 +1836,13 @@ print_single_MC_result
stardis->compute_surface.area * result.SE); /* Standard error */
ERR(sdis_estimator_get_total_flux(estimator, &result));
if(stardis->time_range[0] == stardis->time_range[1])
- fprintf(stream, "Total flux Flux at boundary '%s' at t=%g = %g W +/- %g\n",
+ fprintf(stream, "Total flux at boundary '%s' at t=%g = %g W +/- %g\n",
str_cget(&stardis->solve_name), stardis->time_range[0],
stardis->compute_surface.area * result.E, /* Expected value */
stardis->compute_surface.area * result.SE); /* Standard error */
else
fprintf(stream,
- "Total flux Flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
+ "Total flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n",
str_cget(&stardis->solve_name), SPLIT2(stardis->time_range),
stardis->compute_surface.area * result.E, /* Expected value */
stardis->compute_surface.area * result.SE); /* Standard error */
