commit 9d1152d78345c0db8905fa2c1bf6ecee66f106f0
parent 1ad54ad0656369b371f26136e4f40423149597d5
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 3 Dec 2021 11:43:08 +0100
Merge branch 'feature_non_linear' into develop
Diffstat:
9 files changed, 110 insertions(+), 28 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -38,6 +38,7 @@ set_property(CACHE STARDIS_DOC PROPERTY STRINGS
set(STARDIS_ARGS_DEFAULT_TRAD "300")
set(STARDIS_ARGS_DEFAULT_TRAD_REFERENCE "300")
set(STARDIS_ARGS_DEFAULT_COMPUTE_TIME "INF")
+set(STARDIS_ARGS_DEFAULT_PICARD_ORDER "1")
set(STARDIS_ARGS_DEFAULT_RENDERING_FOV "70") # degrees
set(STARDIS_ARGS_DEFAULT_RENDERING_IMG_HEIGHT "480")
set(STARDIS_ARGS_DEFAULT_RENDERING_IMG_WIDTH "640")
diff --git a/doc/stardis.1.txt.in b/doc/stardis.1.txt.in
@@ -29,12 +29,14 @@ SYNOPSIS
*stardis* *-M* <__file__> [_option_]
DESCRIPTION
------------
-*stardis* solves coupled thermal systems under the linear assumption. Here
-coupled refers to conductive, convective and radiative transfers, and linear
-means that each phenomena is represented using a model that is linear
-with temperature. *stardis* can deal with complex geometries as well as
-high-frequency external solicitations over a very long period of time,
+*stardis* solves coupled thermal systems: conductive, convective and
+radiative transfers are solved together. The physical model used for
+conduction is the local unstationary heat conduction equation.
+Convection fluxes are assumed to be linear with temperature, and radiation
+is assumed to be integrated over the whole thermal spectral range,
+therefore radiative heat fluxes are proportionnal to a difference of
+temperatures to the power 4. *stardis* can deal with complex geometries as
+well as high-frequency external solicitations over a very long period of time,
relative to the characteristic time of the system. The provided system
description should comply with the *stardis-input*(5) format.
@@ -62,13 +64,23 @@ computer graphics technology which has already been a game changer in the
cinema industry (FX and animated movies), this theoretical framework can now
be practically used on the most geometrically complex systems.
-Everytime the linear assumption is relevant, this theoretical framework allows
-to encompass all the heat transfer mechanisms (conductive-convective-radiative)
-in an unified statistical model. Such systems can be solved by a Monte-Carlo
-approach just by sampling heat paths. This can be seen as an extension of
-Monte-Carlo algorithms that solve radiative transfer by sampling optical paths.
-A main property of this approach is that the resulting algorithms does not rely
-on a volume mesh of the system.
+Monte-Carlo algorithms associated with convective and conductive processes
+consist in sampling heat paths: this can be seen as an extension of
+Monte-Carlo algorithms that solve monochromatic radiative transfer.
+The radiative transfer algorithm, based on the Picard method, is also based
+on sampling radiative paths. However, since stardis solves the spectrally
+integrated radiative transfer, the process can be recursive: secondary heat
+paths (convective, conductive and radiative) may be necessary along the
+sampling of an initial radiative path.
+
+The solution may not be sufficiently converged with a Picard order equal
+to 1 in the presence of high temperature gradients.
+Increasing the Picard order may be necessary in this case, until the
+required convergence is reached.
+
+A main property of this approach is that the resulting algorithms do
+not rely on a volumic mesh of the system: only the representation
+of interfaces is necessary.
[1] Delatorre et al., Monte Carlo advances and concentrated solar applications,
Solar Energy, 2014
@@ -238,6 +250,14 @@ different temperature, flux or volumic power values.
Number of Monte-Carlo samples. By default *samples-count* is set to
@STARDIS_ARGS_DEFAULT_SAMPLES_COUNT@.
+*-o* _Picard_order_::
+ Determine the iteration level used with the Picard method to deal with
+ non-linear radiative transfer accross the model.
+ By default *Picard_order* is set to @STARDIS_ARGS_DEFAULT_PICARD_ORDER@.
+ Note that a Picard order greater than 1 is incompatible both with Green
+ computations and models including volumic power sources or non zero flux
+ at a boundary.
+
*-t* _threads-count_::
Hint on the number of threads to use. By default use as many threads as CPU
cores.
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -83,7 +83,9 @@ read_model
stardis->t_range[0] = MMIN(stardis->t_range[0], stardis->trad_ref);
stardis->t_range[1] = MMAX(stardis->t_range[1], stardis->trad_ref);
logger_print(stardis->logger, LOG_OUTPUT,
- "System Tref range is [%g %g]\n", SPLIT2(stardis->t_range));
+ "System T range is [%g %g]\n", SPLIT2(stardis->t_range));
+ logger_print(stardis->logger, LOG_OUTPUT,
+ "Picard order is %u\n", stardis->picard_order);
ASSERT(!f && !txtrdr);
exit:
@@ -239,6 +241,7 @@ stardis_init
stardis->compute_surface.area = 0;
stardis->samples = args->samples;
stardis->nthreads = args->nthreads;
+ stardis->picard_order = args->picard_order;
stardis->scale_factor = -1; /* invalid value */
stardis->trad = STARDIS_DEFAULT_TRAD;
stardis->trad_ref = STARDIS_DEFAULT_TRAD_REFERENCE;
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -793,16 +793,17 @@ struct stardis {
struct str rndgen_state_in_filename;
struct str rndgen_state_out_filename;
struct str chunks_prefix;
- int mode;
- int trad_def;
+ struct mem_allocator* allocator;
+ struct logger* logger;
+ struct sdis_device* dev;
size_t samples;
- unsigned nthreads;
double scale_factor;
double trad, trad_ref;
double t_range[2];
- struct mem_allocator* allocator;
- struct logger* logger;
- struct sdis_device* dev;
+ int mode;
+ int trad_def;
+ unsigned nthreads;
+ unsigned picard_order;
unsigned next_medium_id;
unsigned undefined_medium_behind_boundary_id;
int dump_paths;
diff --git a/src/stardis-args.c b/src/stardis-args.c
@@ -181,9 +181,10 @@ init_args
args->logger = logger;
args->allocator = allocator;
darray_str_init(allocator, &args->model_files);
- /* Set default values */
+ /* Set non-zero default values */
args->samples = STARDIS_DEFAULT_SAMPLES_COUNT;
args->nthreads = SDIS_NTHREADS_DEFAULT;
+ args->picard_order = STARDIS_DEFAULT_PICARD_ORDER;
d2(args->pos_and_time+3,
STARDIS_DEFAULT_COMPUTE_TIME, STARDIS_DEFAULT_COMPUTE_TIME);
args->verbose = STARDIS_DEFAULT_VERBOSE_LEVEL;
@@ -228,13 +229,13 @@ short_help
print_version(stream);
fprintf(stream, "\nMandatory options\n");
- fprintf(stream, "-------------------\n");
+ fprintf(stream, "-----------------\n");
fprintf(stream, "\n -M <FILE>\n");
fprintf(stream, " Read a text file that contains (partial) description of the model.\n");
fprintf(stream, "\nExclusive options\n");
- fprintf(stream, "-------------------\n");
+ fprintf(stream, "-----------------\n");
fprintf(stream, "\n -F STL_FILE[,TIME-RANGE]\n");
fprintf(stream, " Compute the mean flux on a given 2D region at a given time.\n");
@@ -258,7 +259,7 @@ short_help
fprintf(stream, " Compute the by-triangle mean temperature on a given 2D region.\n");
fprintf(stream, "\nOther options\n");
- fprintf(stream, "-------------------\n");
+ fprintf(stream, "-------------\n");
fprintf(stream, "\n -c NAMES_PREFIX\n");
fprintf(stream, " Dump the geometry and property ids to stdout as C chunks.\n");
@@ -284,6 +285,9 @@ short_help
fprintf(stream, "\n -n SAMPLE_COUNT\n");
fprintf(stream, " Set the number of Monte-Carlo samples.\n");
+ fprintf(stream, "\n -o\n");
+ fprintf(stream, " Set the order for the Picard linearization of radiative transfer.\n");
+
fprintf(stream, "\n -t NUM_OF_THREADS\n");
fprintf(stream, " Hint on the number of threads.\n");
@@ -359,9 +363,9 @@ parse_args
struct args* args,
struct mem_allocator* allocator)
{
- int opt = 0, n_used = 0;
+ int opt = 0, n_used = 0, o_used = 0;
size_t len = 0;
- const char option_list[] = "c:dD:eF:gG:hm:M:n:p:P:R:s:S:t:vV:x:X:";
+ const char option_list[] = "c:dD:eF:gG:hm:M:n:o:p:P:R:s:S:t:vV:x:X:";
char buf[128];
struct str keep;
char** line = NULL;
@@ -536,6 +540,23 @@ parse_args
break;
}
+ case 'o': {
+ int order;
+ res = cstr_to_int(optarg, &order);
+ if(res != RES_OK
+ || order <= 0)
+ {
+ if(res == RES_OK) res = RES_BAD_ARG;
+ logger_print(args->logger, LOG_ERROR,
+ "Invalid argument for option -%c: %s\n",
+ opt, optarg);
+ goto error;
+ }
+ args->picard_order = (unsigned)order;
+ o_used = 1;
+ break;
+ }
+
case 'p':
if(args->mode & EXCLUSIVE_MODES) {
res = RES_BAD_ARG;
@@ -728,6 +749,17 @@ parse_args
goto error;
}
+ if(args->mode & (MODE_BIN_GREEN | MODE_GREEN)
+ && o_used && args->picard_order > 1)
+ {
+ logger_print(args->logger, LOG_ERROR,
+ "Option -%c cannot be used if Picard order is not 1 (here order is %u)\n",
+ mode_option(args->mode & (MODE_BIN_GREEN | MODE_GREEN)),
+ args->picard_order);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
if(args->mode & MODE_IR_COMPUTE && n_used) {
logger_print(args->logger, LOG_ERROR,
"The -n option has no effect in rendering mode;"
diff --git a/src/stardis-args.h b/src/stardis-args.h
@@ -96,11 +96,12 @@ struct args {
char* rndgen_state_in_filename;
char* rndgen_state_out_filename;
char* chunks_prefix;
+ char* camera;
size_t samples;
- unsigned nthreads;
double pos_and_time[5];
+ unsigned nthreads;
+ unsigned picard_order;
enum stardis_mode mode;
- char* camera;
enum dump_path_type dump_paths;
int verbose;
};
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -469,6 +469,7 @@ compute_probe(struct stardis* stardis, struct time* start)
}
} else {
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
time_current(&compute_start);
ERR(sdis_solve_probe(stardis->sdis_scn, &args, &estimator));
time_current(&compute_end);
@@ -703,6 +704,7 @@ compute_probe_on_interface(struct stardis* stardis, struct time* start)
}
} else {
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
time_current(&compute_start);
ERR(sdis_solve_probe_boundary(stardis->sdis_scn, &args, &estimator));
time_current(&compute_end);
@@ -851,6 +853,7 @@ compute_camera(struct stardis* stardis, struct time* start)
args.image_resolution[1] = height;
args.spp = (size_t)stardis->camera.spp;
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
/* Launch the simulation */
time_current(&compute_start);
@@ -963,6 +966,7 @@ compute_medium(struct stardis* stardis, struct time* start)
}
} else {
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
time_current(&compute_start);
ERR(sdis_solve_medium(stardis->sdis_scn, &args, &estimator));
time_current(&compute_end);
@@ -1094,6 +1098,7 @@ compute_boundary(struct stardis* stardis, struct time* start)
}
} else {
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
time_current(&compute_start);
ERR(sdis_solve_boundary(stardis->sdis_scn, &args, &estimator));
time_current(&compute_end);
@@ -1138,6 +1143,7 @@ compute_flux_boundary(struct stardis* stardis, struct time* start)
= darray_size_t_cdata_get(&stardis->compute_surface.primitives);
args.nprimitives
= darray_size_t_size_get(&stardis->compute_surface.primitives);
+ args.picard_order = stardis->picard_order;
d2_set(args.time_range, stardis->time_range);
/* Input random state? */
@@ -1199,6 +1205,7 @@ compute_map(struct stardis* stardis, struct time* start)
= darray_size_t_size_get(&stardis->compute_surface.primitives);
d2_set(args.time_range, stardis->time_range);
args.register_paths = stardis->dump_paths;
+ args.picard_order = stardis->picard_order;
ERR(sdis_solve_boundary(stardis->sdis_scn, &args,
darray_estimators_data_get(&estimators) + p));
diff --git a/src/stardis-default.h.in b/src/stardis-default.h.in
@@ -20,6 +20,7 @@
#define STARDIS_DEFAULT_TRAD @STARDIS_ARGS_DEFAULT_TRAD@
#define STARDIS_DEFAULT_TRAD_REFERENCE @STARDIS_ARGS_DEFAULT_TRAD_REFERENCE@
#define STARDIS_DEFAULT_COMPUTE_TIME @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@
+#define STARDIS_DEFAULT_PICARD_ORDER @STARDIS_ARGS_DEFAULT_PICARD_ORDER@
#define STARDIS_DEFAULT_RENDERING_FOV @STARDIS_ARGS_DEFAULT_RENDERING_FOV@
#define STARDIS_DEFAULT_RENDERING_IMG_HEIGHT @STARDIS_ARGS_DEFAULT_RENDERING_IMG_HEIGHT@
#define STARDIS_DEFAULT_RENDERING_IMG_WIDTH @STARDIS_ARGS_DEFAULT_RENDERING_IMG_WIDTH@
diff --git a/src/stardis-parsing.c b/src/stardis-parsing.c
@@ -523,6 +523,14 @@ process_flx
logger_print(stardis->logger, LOG_ERROR, "Invalid flux: %s\n", tk);
goto end;
}
+ if(f_boundary->imposed_flux != 0 && stardis->picard_order > 1) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot have a flux defined at a boundary (here %f) if Picard order "
+ "is not 1 (here order is %u)\n",
+ f_boundary->imposed_flux, stardis->picard_order);
+ res = RES_BAD_ARG;
+ goto end;
+ }
ASSERT(sz <= UINT_MAX);
ERR(read_sides_and_files(stardis, 1, (unsigned)sz, tok_ctx));
@@ -859,6 +867,14 @@ process_solid
logger_print(stardis->logger, LOG_ERROR, "Invalid volumic power: %s\n", tk);
goto end;
}
+ if(solid->vpower != 0 && stardis->picard_order > 1) {
+ logger_print(stardis->logger, LOG_ERROR,
+ "Cannot have volumic power (here %f) if Picard order is not 1 "
+ "(here order is %u)\n",
+ solid->vpower, stardis->picard_order);
+ res = RES_BAD_ARG;
+ goto end;
+ }
/* Actual solid creation is defered until geometry is read to allow
* enclosure shape VS delta analysis (and auto delta computation) */
