commit 4350f88693db871b487a1482edd84d1fe7b9e17b
parent 580c13150738b78675bcf046dbd9b98a35d8ac63
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Mon, 20 May 2019 15:42:33 +0200
Add Green computation and output
Diffstat:
| M | src/args.h | | | 92 | +++++++++++++++++++++++++++++++++++++++++++++++-------------------------------- |
| M | src/main.c | | | 2 | +- |
| M | src/stardis-app.c | | | 76 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------- |
| M | src/stardis-app.h | | | 58 | ++++++++++++++++++++++++++++++++++++++++++++++++++-------- |
| M | src/stardis-compute.c | | | 798 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------ |
5 files changed, 779 insertions(+), 247 deletions(-)
diff --git a/src/args.h b/src/args.h
@@ -11,17 +11,27 @@
#include <sdis_version.h>
enum stardis_mode {
- PROBE_COMPUTE,
- MEDIUM_COMPUTE,
- BOUNDARY_COMPUTE,
- IR_COMPUTE,
- DUMP_VTK,
- UNDEF_MODE
+ UNDEF_MODE = 0,
+ PROBE_COMPUTE = BIT(0),
+ MEDIUM_COMPUTE = BIT(1),
+ BOUNDARY_COMPUTE = BIT(2),
+ GREEN_MODE = BIT(3),
+ IR_COMPUTE = BIT(4),
+ DUMP_VTK = BIT(5),
+ COMPUTE_MODES = PROBE_COMPUTE | MEDIUM_COMPUTE | BOUNDARY_COMPUTE | IR_COMPUTE | DUMP_VTK
};
-static const char MODE_OPTION[UNDEF_MODE] = {
- 'p', 'm', 'b', 'R', 'd'
-};
+static char mode_option(const enum stardis_mode m) {
+ int found = 0;
+ char res = '?';
+ if (m & PROBE_COMPUTE) { found++; res = 'p'; }
+ if (m & MEDIUM_COMPUTE) { found++; res = 'm'; }
+ if (m & BOUNDARY_COMPUTE) { found++; res = 'b'; }
+ if (m & IR_COMPUTE) { found++; res = 'R'; }
+ if (m & DUMP_VTK) { found++; res = 'd'; }
+ ASSERT(found == 1);
+ return res;
+}
enum dump_path_type {
DUMP_NONE = 0,
@@ -67,7 +77,7 @@ print_help(char* prog)
printf("usage: \n%s -M MEDIUM.txt -B BOUNDARY.txt\n", prog);
printf("[ -p X,Y,Z,TIME | -m MEDIUM_NAME,TIME | -d | -b SURFACE.vtk,TIME\n");
printf(" | -R t=TIME:spp=SPP:fov=FOV:up=XUP,YUP,ZUP:pos=X,Y,Z:tgt=XT,YT,ZT:img=WxH ]\n");
- printf("[-s SCALE_FACTOR] [-D {all | error | success}] -d\n");
+ printf("[ -g] [-s SCALE_FACTOR] [-D {all | error | success}] -d\n");
printf("[-n NUM_OF_REALIZATIONS] [-t NUM_OF_THREADS]\n");
printf("[-r AMBIENT_RAD_TEMP:REFERENCE_TEMP]\n");
printf("\n -h: print this help.\n");
@@ -82,21 +92,21 @@ print_help(char* prog)
printf("\nExclusive arguments\n");
printf("-------------------\n");
printf("\n -p X:Y:Z:TIME:\n");
- printf(" compute the temperature at the given probe.\n");
+ printf(" Compute the temperature at the given probe.\n");
printf(" The probe must be in a medium. If some boundary conditions are\n");
- printf(" time-dependant, TIME can not be INF.\n");
+ printf(" time-dependant, TIME cannot be INF.\n");
printf("\n -m MEDIUM_NAME,TIME:\n");
- printf(" compute the mean temperature in a given medium at a given time.\n");
- printf(" If some boundary conditions are time-dependant, TIME can not be INF.\n");
+ printf(" Compute the mean temperature in a given medium at a given time.\n");
+ printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n");
printf("\n -d:\n");
- printf(" dump the geometry in VTK format with medium front/back id and\n");
+ printf(" Dump the geometry in VTK format with medium front/back id and\n");
printf(" boundary id.\n");
printf("\n -b SURFACE.vtk,TIME:\n");
- printf(" compute the mean temperature on a given 2D shape at a given time.\n");
+ printf(" Compute the mean temperature on a given 2D shape at a given time.\n");
printf(" Mean temperature is computed on the front sides of the primitives\n");
printf(" listed in SURFACE.vtk. These primitives are not added to the geometry,\n");
printf(" but must be already known. The shape doesn't need to be connex.\n");
- printf(" If some boundary conditions are time-dependant, TIME can not be INF.\n");
+ printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n");
printf("\n -R t=TIME:spp=SPP:fov=FOV:up=XUP,YUP,ZUP:pos=X,Y,Z:tgt=XT,YT,ZT:img=WxH:\n");
printf(" The infra-red rendering mode.\n");
printf(" t is the rendering time (default: INF).\n");
@@ -108,20 +118,24 @@ print_help(char* prog)
printf(" img is the resolution (default: 640x480).\n");
printf("\nOptionnal arguments\n");
printf("-------------------\n");
+ printf("\n -g:\n");
+ printf(" Compute the green function.\n");
+ printf(" Must be used in conjonction with one the -p, -m or -b options.\n");
+ printf(" Provided TIME is silently ignored.\n");
printf("\n -s SCALE_FACTOR:\n");
printf(" default value: 1.0.\n");
printf("\n -D { all | error | success }:\n");
printf(" dump paths in VTK format.\n");
printf("\n -n NUM_OF_REALIZATIONS:\n");
- printf(" number of Monte-Carlo realizations.\n");
+ printf(" Number of Monte-Carlo realizations.\n");
printf(" default value: 10000.\n");
printf("\n -t NUM_OF_THREADS:\n");
- printf(" default value: all threads available.\n");
+ printf(" Number of threads; default value is all threads available.\n");
printf("\n -r AMBIENT_RAD_TEMP,REFERENCE_TEMP:\n");
printf(" Parameters for the radiative transfer.\n");
printf(" AMBIENT_RAD_TEMP is the ambient radiative temperature .\n");
- printf(" REFERENCE_TEMP is the temperature used for the linearization.\n");
- printf(" of theradiative transfer.\n");
+ printf(" REFERENCE_TEMP is the temperature used for the linearization\n");
+ printf(" of the radiative transfer.\n");
}
@@ -138,13 +152,17 @@ parse_args(const int argc, char** argv, struct args* args)
goto error;
}
- while((opt = getopt(argc, argv, "hn:t:b:B:M:m:p:dD:s:r:R:")) != -1) {
+ while((opt = getopt(argc, argv, "hgn:t:b:B:M:m:p:dD:s:r:R:")) != -1) {
switch(opt) {
case 'h':
print_help(argv[0]);
args->just_help = 1;
goto exit;
+ case 'g':
+ args->mode |= GREEN_MODE;
+ break;
+
case 'n': {
args->N = strtoull(optarg, NULL, 10);
if (args->N <= 0) {
@@ -174,13 +192,13 @@ parse_args(const int argc, char** argv, struct args* args)
break;
case 'p':
- if (args->mode != UNDEF_MODE) {
+ if (args->mode & COMPUTE_MODES) {
res = RES_BAD_ARG;
fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, MODE_OPTION[args->mode]);
+ (char)opt, mode_option(args->mode));
goto error;
}
- args->mode = PROBE_COMPUTE;
+ args->mode |= PROBE_COMPUTE;
cstr_to_list_double(optarg, ',', args->u.probe, &len, 4);
if(len != 4
|| res != RES_OK) {
@@ -193,13 +211,13 @@ parse_args(const int argc, char** argv, struct args* args)
case 'b': {
int err = 0;
char *ptr;
- if (args->mode != UNDEF_MODE) {
+ if (args->mode & COMPUTE_MODES) {
res = RES_BAD_ARG;
fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, MODE_OPTION[args->mode]);
+ (char)opt, mode_option(args->mode));
goto error;
}
- args->mode = BOUNDARY_COMPUTE;
+ args->mode |= BOUNDARY_COMPUTE;
ptr = strrchr(optarg, ',');
/* Single ',' allowed */
if (ptr != strchr(optarg, ','))
@@ -219,13 +237,13 @@ parse_args(const int argc, char** argv, struct args* args)
}
case 'm': {
char* ptr;
- if (args->mode != UNDEF_MODE) {
+ if (args->mode & COMPUTE_MODES) {
res = RES_BAD_ARG;
fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, MODE_OPTION[args->mode]);
+ (char)opt, mode_option(args->mode));
goto error;
}
- args->mode = MEDIUM_COMPUTE;
+ args->mode |= MEDIUM_COMPUTE;
ptr = strpbrk(optarg, ",");
if (! ptr || ptr == optarg)
res = RES_BAD_ARG;
@@ -240,13 +258,13 @@ parse_args(const int argc, char** argv, struct args* args)
}
case 'd':
- if (args->mode != UNDEF_MODE) {
+ if (args->mode & COMPUTE_MODES) {
res = RES_BAD_ARG;
fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, MODE_OPTION[args->mode]);
+ (char)opt, mode_option(args->mode));
goto error;
}
- args->mode = DUMP_VTK;
+ args->mode |= DUMP_VTK;
break;
case 'D':
@@ -286,13 +304,13 @@ parse_args(const int argc, char** argv, struct args* args)
break;
case 'R':
- if (args->mode != UNDEF_MODE) {
+ if (args->mode & COMPUTE_MODES) {
res = RES_BAD_ARG;
fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, MODE_OPTION[args->mode]);
+ (char)opt, mode_option(args->mode));
goto error;
}
- args->mode = IR_COMPUTE;
+ args->mode |= IR_COMPUTE;
args->camera = optarg;
break;
}
diff --git a/src/main.c b/src/main.c
@@ -20,7 +20,7 @@ int main(int argc, char** argv){
res = stardis_init(&args, &stardis);
if (res != RES_OK) goto error;
- if (args.mode == DUMP_VTK) {
+ if (args.mode & DUMP_VTK) {
res = dump_vtk(stdout, &stardis.geometry);
if (res != RES_OK) goto error;
goto exit;
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -135,12 +135,12 @@ parse_medium_line(char* line, char** stl_filename, struct description* desc)
desc->type = DESC_MAT_SOLID;
CHK_TOK(strtok(NULL, " "), "medium name");
- if (strlen(tk) >= sizeof(desc->name)) {
+ if (strlen(tk) >= sizeof(desc->d.solid.name)) {
fprintf(stderr, "Medium name is too long: %s\n", tk);
res = RES_BAD_ARG;
goto exit;
}
- strncpy(desc->name, tk, sizeof(desc->name));
+ strncpy(desc->d.solid.name, tk, sizeof(desc->d.solid.name));
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
@@ -177,8 +177,9 @@ parse_medium_line(char* line, char** stl_filename, struct description* desc)
/* Depending of the number of spaces before '"' strtok should
* be called once or twice */
CHK_TOK(strtok(NULL, "\""), "Tinit");
- if(*(tk + strspn(tk, " \t")) == '\0')
+ if (*(tk + strspn(tk, " \t")) == '\0') {
CHK_TOK(strtok(NULL, "\""), "Tinit");
+ }
desc->d.solid.Tinit = malloc(strlen(tk) + 1);
strcpy(desc->d.solid.Tinit, tk);
/* Closing " fot Tinit; can return NULL if line ends just after "Tinit" */
@@ -203,12 +204,12 @@ parse_medium_line(char* line, char** stl_filename, struct description* desc)
desc->type = DESC_MAT_FLUID;
CHK_TOK(strtok(NULL, " "), "medium name");
- if (strlen(tk) >= sizeof(desc->name)) {
+ if (strlen(tk) >= sizeof(desc->d.fluid.name)) {
fprintf(stderr, "Medium name is too long: %s\n", tk);
res = RES_BAD_ARG;
goto exit;
}
- strncpy(desc->name, tk, sizeof(desc->name));
+ strncpy(desc->d.fluid.name, tk, sizeof(desc->d.fluid.name));
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
@@ -231,8 +232,9 @@ parse_medium_line(char* line, char** stl_filename, struct description* desc)
/* Depending of the number of spaces before '"' strtok should
* be called once or twice */
CHK_TOK(strtok(NULL, "\""), "Tinit");
- if (*(tk + strspn(tk, " \t")) == '\0')
+ if (*(tk + strspn(tk, " \t")) == '\0') {
CHK_TOK(strtok(NULL, "\""), "Tinit");
+ }
desc->d.fluid.Tinit = malloc(strlen(tk) + 1);
strcpy(desc->d.fluid.Tinit, tk);
}
@@ -247,12 +249,13 @@ exit:
}
/* Read boundary line; should be one of:
-# H_BOUNDARY_FOR_SOLID STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
-# | H_BOUNDARY_FOR_FLUID STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
-# | T_BOUNDARY_FOR_SOLID STL_filename "T_value(x,y,z,t)"
-# | T_BOUNDARY_FOR_FLUID STL_filename "T_value(x,y,z,t)" hc hc_max
-# | F_BOUNDARY_FOR_SOLID STL_filename "F_value(x,y,z,t)"
+# H_BOUNDARY_FOR_SOLID Boundary_name STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
+# | H_BOUNDARY_FOR_FLUID Boundary_name STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
+# | T_BOUNDARY_FOR_SOLID Boundary_name STL_filename "T_value(x,y,z,t)"
+# | T_BOUNDARY_FOR_FLUID Boundary_name STL_filename "T_value(x,y,z,t)" hc hc_max
+# | F_BOUNDARY_FOR_SOLID Boundary_name STL_filename "F_value(x,y,z,t)"
# ; no F_BOUNDARY_FOR_FLUID
+# | SOLID_FLUID_CONNECTION Connection_name STL_filename emissivity specular_fraction hc hc_max
*/
static res_T
parse_boundary_line(char* line, char** stl_filename, struct description* desc)
@@ -286,9 +289,18 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
}
}
}
+
if (h_solid || h_fluid) {
desc->type = h_solid ? DESC_BOUND_H_FOR_SOLID : DESC_BOUND_H_FOR_FLUID;
+ CHK_TOK(strtok(NULL, " "), "boundary name");
+ if (strlen(tk) >= sizeof(desc->d.h_boundary.name)) {
+ fprintf(stderr, "Boundary name is too long: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto exit;
+ }
+ strncpy(desc->d.h_boundary.name, tk, sizeof(desc->d.h_boundary.name));
+
/* The description is parsed the same way for both types */
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
@@ -332,14 +344,23 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
/* Depending of the number of spaces before '"' strtok should
* be called once or twice */
CHK_TOK(strtok(NULL, "\""), "temperature");
- if (*(tk + strspn(tk, " \t")) == '\0')
+ if (*(tk + strspn(tk, " \t")) == '\0') {
CHK_TOK(strtok(NULL, "\""), "temperature");
+ }
desc->d.h_boundary.T = malloc(strlen(tk) + 1);
strcpy(desc->d.h_boundary.T, tk);
}
else if (t_solid || t_fluid) {
desc->type = t_solid ? DESC_BOUND_T_FOR_SOLID : DESC_BOUND_T_FOR_FLUID;
+ CHK_TOK(strtok(NULL, " "), "boundary name");
+ if (strlen(tk) >= sizeof(desc->d.t_boundary.name)) {
+ fprintf(stderr, "Boundary name is too long: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto exit;
+ }
+ strncpy(desc->d.t_boundary.name, tk, sizeof(desc->d.t_boundary.name));
+
/* This part of the description is parsed the same way for both types */
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
@@ -348,8 +369,9 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
/* Depending of the number of spaces before '"' strtok should
* be called once or twice */
CHK_TOK(strtok(NULL, "\""), "temperature");
- if (*(tk + strspn(tk, " \t")) == '\0')
+ if (*(tk + strspn(tk, " \t")) == '\0') {
CHK_TOK(strtok(NULL, "\""), "temperature");
+ }
desc->d.t_boundary.T = malloc(strlen(tk) + 1);
strcpy(desc->d.t_boundary.T, tk);
@@ -378,6 +400,14 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
else if (f_solid) {
desc->type = DESC_BOUND_F_FOR_SOLID;
+ CHK_TOK(strtok(NULL, " "), "boundary name");
+ if (strlen(tk) >= sizeof(desc->d.f_boundary.name)) {
+ fprintf(stderr, "Boundary name is too long: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto exit;
+ }
+ strncpy(desc->d.f_boundary.name, tk, sizeof(desc->d.f_boundary.name));
+
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
strcpy(*stl_filename, tk);
@@ -385,8 +415,9 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
/* Depending of the number of spaces before '"' strtok should
* be called once or twice */
CHK_TOK(strtok(NULL, "\""), "flux");
- if (*(tk + strspn(tk, " \t")) == '\0')
+ if (*(tk + strspn(tk, " \t")) == '\0') {
CHK_TOK(strtok(NULL, "\""), "flux");
+ }
desc->d.f_boundary.flux = malloc(strlen(tk) + 1);
strcpy(desc->d.f_boundary.flux, tk);
@@ -414,6 +445,14 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
else if (sf_connect) {
desc->type = DESC_SOLID_FLUID_CONNECT;
+ CHK_TOK(strtok(NULL, " "), "boundary name");
+ if (strlen(tk) >= sizeof(desc->d.sf_connect.name)) {
+ fprintf(stderr, "Boundary name is too long: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto exit;
+ }
+ strncpy(desc->d.sf_connect.name, tk, sizeof(desc->d.sf_connect.name));
+
CHK_TOK(strtok(NULL, " "), "file name");
*stl_filename = malloc(strlen(tk) + 1);
strcpy(*stl_filename, tk);
@@ -453,8 +492,7 @@ parse_boundary_line(char* line, char** stl_filename, struct description* desc)
res = RES_BAD_ARG;
goto exit;
}
- }
- else {
+ } else {
fprintf(stderr, "Invalid boundary type: %s", tk);
res = RES_BAD_ARG;
goto exit;
@@ -897,14 +935,14 @@ stardis_init
if (args->dump_paths & DUMP_SUCCESS)
stardis->dump_paths |= SDIS_HEAT_PATH_SUCCEED;
- if (args->mode == IR_COMPUTE){
+ if (args->mode & IR_COMPUTE){
res = parse_camera(args->camera, &stardis->camera);
if (res != RES_OK) goto error;
}
- else if (args->mode == MEDIUM_COMPUTE) {
+ else if (args->mode & MEDIUM_COMPUTE) {
strcpy(stardis->solve_name, args->medium_name);
}
- else if (args->mode == BOUNDARY_COMPUTE) {
+ else if (args->mode & BOUNDARY_COMPUTE) {
strcpy(stardis->solve_name, args->solve_boundary_filename);
}
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -122,7 +122,7 @@ release_geometry(struct geometry* geom)
}
struct mat_fluid {
- char name[16];
+ char name[32];
unsigned fluid_id;
double rho;
double cp;
@@ -217,7 +217,7 @@ hash_fluid(const struct mat_fluid* key)
}
struct mat_solid {
- char name[16];
+ char name[32];
unsigned solid_id;
double lambda;
double rho;
@@ -355,6 +355,7 @@ hash_solid(const struct mat_solid* key)
}
struct h_boundary {
+ char name[32];
unsigned mat_id;
double emissivity;
double specular_fraction;
@@ -373,7 +374,8 @@ print_h_boundary
ASSERT(stream && b
&& (type == DESC_BOUND_H_FOR_SOLID || type == DESC_BOUND_H_FOR_FLUID));
fprintf(stream,
- "H boundary for %s: emissivity=%g specular_fraction=%g hc=%g hc_max=%g T='%s'\n",
+ "H boundary %s for %s: emissivity=%g specular_fraction=%g hc=%g hc_max=%g T='%s'\n",
+ b->name,
(type == DESC_BOUND_H_FOR_SOLID ? "solid" : "fluid"),
(b->has_emissivity ? b->emissivity : 0),
(b->has_emissivity ? b->specular_fraction : 0),
@@ -388,6 +390,7 @@ eq_h_boundary
const struct h_boundary* b)
{
if (a->mat_id != b->mat_id
+ || strcmp(a->name, b->name)
|| a->specular_fraction != b->specular_fraction
|| strcmp(a->T, b->T)
|| a->has_emissivity != b->has_emissivity
@@ -409,6 +412,10 @@ hash_h_boundary(const struct h_boundary* key)
uint32_t hash = OFFSET32_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV32_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV32_PRIME;
@@ -451,6 +458,10 @@ hash_h_boundary(const struct h_boundary* key)
uint64_t hash = OFFSET64_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV64_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV64_PRIME;
@@ -491,6 +502,7 @@ hash_h_boundary(const struct h_boundary* key)
}
struct t_boundary {
+ char name[32];
unsigned mat_id;
char* T;
struct te_expr* te_temperature;
@@ -508,7 +520,8 @@ print_t_boundary
ASSERT(stream && b
&& (type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID));
fprintf(stream,
- "T boundary for %s: hc=%g hc_max=%g T='%s'\n",
+ "T boundary %s for %s: hc=%g hc_max=%g T='%s'\n",
+ b->name,
(type == DESC_BOUND_T_FOR_SOLID ? "solid" : "fluid"),
(b->has_hc ? b->hc : 0),
(b->has_hc ? b->hc_max : 0),
@@ -523,6 +536,8 @@ eq_t_boundary
{
ASSERT(type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID);
/* These fields are meaningful by both types */
+ if (strcmp(a->name, b->name))
+ return 0;
if (a->mat_id != b->mat_id || strcmp(a->T, b->T))
return 0;
if (type != DESC_BOUND_T_FOR_FLUID)
@@ -550,6 +565,10 @@ hash_t_boundary
uint32_t hash = OFFSET32_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV32_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV32_PRIME;
@@ -580,6 +599,10 @@ hash_t_boundary
uint64_t hash = OFFSET64_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV64_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV64_PRIME;
@@ -608,6 +631,7 @@ hash_t_boundary
}
struct f_boundary {
+ char name[32];
unsigned mat_id;
char* flux;
struct te_expr* te_flux;
@@ -624,7 +648,8 @@ print_f_boundary
{
ASSERT(stream && b && type == DESC_BOUND_F_FOR_SOLID); (void)type;
fprintf(stream,
- "F boundary for SOLID: hc=%g hc_max=%g flux='%s'\n",
+ "F boundary %s for SOLID: hc=%g hc_max=%g flux='%s'\n",
+ b->name,
(b->has_hc ? b->hc : 0),
(b->has_hc ? b->hc_max : 0),
(b->flux ? b->flux : "0"));
@@ -634,6 +659,7 @@ static char
eq_f_boundary(const struct f_boundary* a, const struct f_boundary* b)
{
if (a->mat_id != b->mat_id
+ || strcmp(a->name, b->name)
|| strcmp(a->flux, b->flux)
|| a->has_hc != b->has_hc)
return 0;
@@ -652,6 +678,10 @@ hash_f_boundary(const struct f_boundary* key)
uint32_t hash = OFFSET32_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV32_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV32_PRIME;
@@ -681,6 +711,10 @@ hash_f_boundary(const struct f_boundary* key)
uint64_t hash = OFFSET64_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV64_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->mat_id)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
hash = hash * FNV64_PRIME;
@@ -708,6 +742,7 @@ hash_f_boundary(const struct f_boundary* key)
}
struct solid_fluid_connect {
+ char name[32];
double emissivity;
double specular_fraction;
double hc;
@@ -718,7 +753,8 @@ struct solid_fluid_connect {
static char
eq_sf_connect(const struct solid_fluid_connect* a, const struct solid_fluid_connect* b)
{
- return (char) (a->emissivity == b->emissivity
+ return (char)(!strcmp(a->name, b->name)
+ && a->emissivity == b->emissivity
&& a->specular_fraction == b->specular_fraction
&& a->hc == b->hc
&& a->hc_max == b->hc_max
@@ -737,6 +773,10 @@ hash_sf_connect(const struct solid_fluid_connect* key)
uint32_t hash = OFFSET32_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV32_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->emissivity)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]);
hash = hash * FNV32_PRIME;
@@ -770,6 +810,10 @@ hash_sf_connect(const struct solid_fluid_connect* key)
uint64_t hash = OFFSET64_BASIS;
size_t i;
ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->name)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]);
+ hash = hash * FNV64_PRIME;
+ }
FOR_EACH(i, 0, sizeof(key->emissivity)) {
hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]);
hash = hash * FNV64_PRIME;
@@ -801,7 +845,6 @@ hash_sf_connect(const struct solid_fluid_connect* key)
}
struct description {
- char name[16];
enum description_type type;
union d {
struct mat_fluid fluid;
@@ -816,7 +859,6 @@ struct description {
FINLINE void
init_description(struct description* desc) {
ASSERT(desc);
- desc->name[0] = '\0';
desc->type = DESCRIPTION_TYPE_COUNT__;
}
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -15,6 +15,13 @@
#define HTABLE_KEY unsigned
#include <rsys/hash_table.h>
+/* A type to limit variable use in tinyexpr expressions */
+enum var_prohibited_t {
+ ALL_VARS_ALLOWED = 0,
+ T_PROHIBITED = BIT(1),
+ XYZ_PROHIBITED = BIT(2)
+};
+
static void
geometry_get_position
(const size_t ivert,
@@ -53,7 +60,7 @@ static res_T
compile_expr_to_fn
(struct te_expr** f,
const char* math_expr,
- const int t_allowed,
+ const enum var_prohibited_t prohibited,
int* is_zero)
{
te_variable vars[4] = {
@@ -62,12 +69,15 @@ compile_expr_to_fn
TE_DEF_OFFSET("z", offsetof(struct sdis_rwalk_vertex, P[2])),
TE_DEF_OFFSET("t", offsetof(struct sdis_rwalk_vertex, time))
};
-
+ int fst = 0, lst = 4;
ASSERT(math_expr);
- *f = te_compile(math_expr, vars, t_allowed ? 4 : 3, NULL);
+ if (prohibited & XYZ_PROHIBITED) fst = 3;
+ if(prohibited & T_PROHIBITED) lst = 3;
+ *f = te_compile(math_expr, vars+fst, lst-fst, NULL);
if (!*f) {
- if (!t_allowed && te_compile(math_expr, vars, 4, NULL))
- fprintf(stderr, "Expression use variable t and should not.\n");
+ if ((prohibited != ALL_VARS_ALLOWED)
+ && te_compile(math_expr, vars, 4, NULL))
+ fprintf(stderr, "Expression %s use a probibited variable.\n", math_expr);
return RES_BAD_ARG;
}
if (is_zero) *is_zero = !((*f)->type == TE_CONSTANT && (*f)->v.value == 0);
@@ -80,8 +90,13 @@ compile_expr_to_fn
struct fluid {
double cp; /* Calorific capacity */
double rho; /* Volumic mass */
+ /* Compute mode */
+ int is_green, is_outside;
+ double t0;
/* TinyExpr stuff to compute temperature */
struct te_expr *temperature;
+ /* ID */
+ unsigned id;
};
static double
@@ -115,7 +130,7 @@ fluid_get_tinit
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
const struct fluid* fluid_props = sdis_data_cget(data);
- if(vtx->time > 0) {
+ if (fluid_props->is_green || vtx->time > fluid_props->t0) {
return -1;
}
return te_eval(fluid_props->temperature, vtx);
@@ -136,9 +151,14 @@ struct solid {
double lambda; /* Conductivity */
double rho; /* Volumic mass */
double delta; /* Numerical parameter */
+ /* Compute mode */
+ int is_green, is_outside;
+ double t0;
/* TinyExpr stuff to compute temperature & power */
struct te_expr *t_init;
struct te_expr *power;
+ /* ID */
+ unsigned id;
};
static double
@@ -213,7 +233,7 @@ solid_get_tinit
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
const struct solid* solid_props = sdis_data_cget(data);
- if (vtx->time > 0) {
+ if (solid_props->is_green || vtx->time > solid_props->t0) {
return -1;
}
return te_eval(solid_props->t_init, vtx);
@@ -245,6 +265,8 @@ struct intface {
/* TinyExpr stuff to compute temperature & flux */
struct te_expr *temperature;
struct te_expr *flux;
+ /* IDs */
+ unsigned front_boundary_id, back_boundary_id;
};
static double
@@ -426,7 +448,9 @@ create_fluid
(struct sdis_device* dev,
const double rho,
const double cp,
- const int t_allowed,
+ const enum var_prohibited_t prohibited,
+ const int is_green,
+ const int is_outside,
const char* t_expr,
struct sdis_medium*** media_ptr,
sdis_medium_getter_T get_temp,
@@ -445,10 +469,17 @@ create_fluid
res = sdis_data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid),
release_fluid_data, &data);
if (res != RES_OK) goto error;
+
+ sz = sa_size(*media_ptr);
+ ASSERT(sz < INT_MAX);
fluid_props = sdis_data_get(data); /* Fetch the allocated memory space */
fluid_props->cp = cp;
fluid_props->rho = rho;
- res = compile_expr_to_fn(&fluid_props->temperature, t_expr, t_allowed, NULL);
+ fluid_props->is_green = is_green;
+ fluid_props->is_outside = is_outside;
+ fluid_props->t0 = 0;
+ fluid_props->id = (unsigned)sz;
+ res = compile_expr_to_fn(&fluid_props->temperature, t_expr, prohibited, NULL);
if (res != RES_OK) {
fprintf(stderr, "Invalid initial temperature expression: %s\n",
t_expr);
@@ -456,9 +487,7 @@ create_fluid
}
res = sdis_fluid_create(dev, &fluid_shader, data, sa_add(*media_ptr, 1));
if (res != RES_OK) goto error;
- sz = sa_size(*media_ptr) - 1;
- ASSERT(sz < INT_MAX);
- *out_id = (unsigned)sz;
+ *out_id = fluid_props->id;
end:
if (data) SDIS(data_ref_put(data));
@@ -474,7 +503,9 @@ create_solid
const double rho,
const double cp,
const double delta,
- const int t_allowed,
+ const enum var_prohibited_t prohibited,
+ const int is_green,
+ const int is_outside,
const char* t_expr, /* Can be NULL if not used by getter */
const char* power_expr, /* Can be NULL */
struct sdis_medium*** media_ptr,
@@ -502,6 +533,8 @@ create_solid
release_solid_data, &data);
if (res != RES_OK) goto error;
+ sz = sa_size(*media_ptr);
+ ASSERT(sz < INT_MAX);
solid_props = sdis_data_get(data); /* Fetch the allocated memory space */
solid_props->lambda = lambda;
solid_props->rho = rho;
@@ -509,8 +542,12 @@ create_solid
solid_props->delta = delta;
solid_props->t_init = NULL;
solid_props->power = NULL;
+ solid_props->is_green = is_green;
+ solid_props->is_outside = is_outside;
+ solid_props->t0 = 0;
+ solid_props->id = (unsigned)sz;
if (t_expr) {
- res = compile_expr_to_fn(&solid_props->t_init, t_expr, t_allowed, NULL);
+ res = compile_expr_to_fn(&solid_props->t_init, t_expr, prohibited, NULL);
if (res != RES_OK) {
fprintf(stderr, "Invalid initial temperature expression: %s\n",
t_expr);
@@ -533,9 +570,7 @@ create_solid
}
res = sdis_solid_create(dev, &solid_shader, data, sa_add(*media_ptr, 1));
if (res != RES_OK) goto error;
- sz = sa_size(*media_ptr) - 1;
- ASSERT(sz < INT_MAX);
- *out_id = (unsigned)sz;
+ *out_id = solid_props->id;
end:
if (data) SDIS(data_ref_put(data));
@@ -751,6 +786,188 @@ dump_path
return RES_OK;
}
+static res_T
+print_power_term
+ (struct sdis_medium* mdm,
+ const double power_term,
+ void* ctx)
+{
+ struct sdis_data* data = NULL;
+ enum sdis_medium_type type;
+ unsigned id;
+ struct description* desc = ctx;
+
+ CHK(mdm && ctx);
+
+ data = sdis_medium_get_data(mdm);
+ type = sdis_medium_get_type(mdm);
+
+ switch (type) {
+ case SDIS_FLUID: {
+ FATAL("Unexpected power term in fluid");
+#if 0
+ struct fluid* d__ = sdis_data_get(data);
+ id = d__->id;
+ desc = (struct description*)ctx + id;
+ ASSERT(id == desc->d.fluid.fluid_id);
+ printf("\tF\t%u\t%g", id, power_term);
+ break;
+#endif
+ }
+ case SDIS_SOLID: {
+ struct solid* d__ = sdis_data_get(data);
+ id = d__->id;
+ desc = (struct description*)ctx + id;
+ ASSERT(id == desc->d.solid.solid_id);
+ printf("\tS\t%u\t%g", id, power_term);
+ break;
+ }
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return RES_OK;
+}
+
+static res_T
+print_flux_term
+ (struct sdis_interface* interf,
+ const enum sdis_side side,
+ const double flux_term,
+ void* ctx)
+{
+ struct sdis_data* data = NULL;
+ struct intface* d__;
+ unsigned id;
+ struct description* desc;
+
+ CHK(interf && ctx);
+
+ data = sdis_interface_get_data(interf);
+ d__ = sdis_data_get(data);
+ id = (side == SDIS_FRONT) ? d__->front_boundary_id : d__->back_boundary_id;
+ desc = (struct description*)ctx + id;
+
+ switch (desc->type) {
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ case DESC_BOUND_H_FOR_SOLID:
+ case DESC_BOUND_H_FOR_FLUID:
+ FATAL("Cannot have a flux term here.\n"); break;
+ case DESC_BOUND_F_FOR_SOLID:
+ ASSERT(id == desc->d.f_boundary.mat_id);
+ printf("\t%u\t%g", id, flux_term);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return RES_OK;
+}
+
+static res_T
+print_sample
+ (struct sdis_green_path* path,
+ void* ctx)
+{
+ struct sdis_point pt = SDIS_POINT_NULL;
+ struct sdis_data* data = NULL;
+ enum sdis_medium_type type;
+ struct description* desc;
+ unsigned id;
+ size_t pcount, fcount;
+ res_T res = RES_OK;
+ CHK(path && ctx);
+
+ res = sdis_green_path_get_limit_point(path, &pt);
+ if (res != RES_OK) goto error;
+
+
+ /* For each path, prints:
+ * # end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n
+ * with:
+ * - end = end_type end_id; end_type = T | H | X | R | F | S
+ * - power_term_i = power_type_i power_id_i factor_i
+ * - flux_term_i = flux_id_i factor_i
+ */
+
+ switch (pt.type) {
+ case SDIS_FRAGMENT: {
+ struct intface* d__;
+ data = sdis_interface_get_data(pt.data.itfrag.intface);
+ d__ = sdis_data_get(data);
+ id = (pt.data.itfrag.fragment.side == SDIS_FRONT)
+ ? d__->front_boundary_id : d__->back_boundary_id;
+ ASSERT(id != UINT_MAX);
+ desc = (struct description*)ctx + id;
+ switch (desc->type) {
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ ASSERT(id == desc->d.t_boundary.mat_id);
+ printf("T\t%u", id);
+ break;
+ case DESC_BOUND_H_FOR_SOLID:
+ case DESC_BOUND_H_FOR_FLUID:
+ ASSERT(id == desc->d.h_boundary.mat_id);
+ printf("H\t%u", id);
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ ASSERT(id == desc->d.f_boundary.mat_id);
+ printf("X\t%u", id);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ break;
+ }
+ case SDIS_VERTEX:
+ type = sdis_medium_get_type(pt.data.mdmvert.medium);
+ data = sdis_medium_get_data(pt.data.mdmvert.medium);
+ if (pt.data.mdmvert.vertex.P[0] == INF) {
+ /* Radiative output */
+ printf("R %u", (unsigned)sa_size(ctx));
+ }
+ else if (type == SDIS_FLUID) {
+ struct fluid* d__ = sdis_data_get(data);
+ id = d__->id;
+ desc = (struct description*)ctx + id;
+ ASSERT(id == desc->d.fluid.fluid_id);
+ if (d__->is_outside)
+ /* If outside the model and in a fluid with known temperature,
+ * its a fluid attached to a H boundary */
+ printf("H\t%u", id);
+ /* In a standard fluid with known temperature */
+ else printf("F\t%u", id);
+
+ } else {
+ struct solid* d__ = sdis_data_get(data);
+ ASSERT(type == SDIS_SOLID);
+ id = d__->id;
+ desc = (struct description*)ctx + id;
+ ASSERT(id == desc->d.solid.solid_id);
+ desc = (struct description*)ctx + id;
+ ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */
+ printf("S\t%u", id);
+ }
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ res = sdis_green_function_get_power_terms_count(path, &pcount);
+ if (res != RES_OK) goto error;
+
+ res = sdis_green_function_get_flux_terms_count(path, &fcount);
+ if (res != RES_OK) goto error;
+
+ printf("\t%zu\t%zu", pcount, fcount);
+
+ res = sdis_green_path_for_each_power_term(path, print_power_term, ctx);
+ if (res != RES_OK) goto error;
+
+ res = sdis_green_path_for_each_flux_term(path, print_flux_term, ctx);
+ if (res != RES_OK) goto error;
+
+ printf("\n");
+
+error:
+ return res;
+}
+
res_T
stardis_compute(struct stardis* stardis, enum stardis_mode mode)
{
@@ -772,8 +989,9 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
double pos[3] = { 0,0,0 };
double time[2] = { 0, 0 };
size_t nfailures;
+ size_t smed_count = 0, fmed_count = 0,
+ tbound_count = 0, hbound_count = 0, fbound_count = 0, sfconnect_count = 0;
unsigned i = 0;
- const int t_allowed = stardis->probe[3] < INF; /* Can condition use the variable t? */
res = sdis_device_create(NULL, &stardis->allocator, stardis->nthreads, 1, &dev);
if (res != RES_OK) goto error;
@@ -781,22 +999,31 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
/* Create media and property holders found in descriptions */
for (i = 0; i < sa_size(stardis->descriptions); ++i) {
struct description* desc = stardis->descriptions + i;
+ /* Base prohibition; additional restrictions apply for media */
+ enum var_prohibited_t prohibited =
+ ((stardis->probe[3] < INF) | (mode & GREEN_MODE))
+ ? T_PROHIBITED : ALL_VARS_ALLOWED;
switch (desc->type) {
case DESC_BOUND_H_FOR_SOLID:
/* Create an external fluid */
- res = create_fluid(dev, 1, 1, t_allowed, desc->d.h_boundary.T, &media,
- fluid_get_temperature, &desc->d.h_boundary.mat_id);
+ hbound_count++;
+ res = create_fluid(dev, 1, 1, prohibited, (mode & GREEN_MODE), 1,
+ desc->d.h_boundary.T, &media, fluid_get_temperature,
+ &desc->d.h_boundary.mat_id);
if (res != RES_OK) goto error;
break;
case DESC_BOUND_H_FOR_FLUID:
/* Create an external solid */
- res = create_solid(dev, INF, 1, 1, 1, t_allowed, desc->d.h_boundary.T,
- NULL, &media, solid_get_temperature, NULL, &desc->d.h_boundary.mat_id);
+ hbound_count++;
+ res = create_solid(dev, INF, 1, 1, 1, prohibited, (mode & GREEN_MODE), 1,
+ desc->d.h_boundary.T, NULL, &media, solid_get_temperature, NULL,
+ &desc->d.h_boundary.mat_id);
if (res != RES_OK) goto error;
break;
case DESC_BOUND_T_FOR_SOLID:
case DESC_BOUND_T_FOR_FLUID:
+ tbound_count++;
ASSERT(desc->d.t_boundary.T != NULL);
res = compile_expr_to_fn(&desc->d.t_boundary.te_temperature,
desc->d.t_boundary.T, 1, NULL);
@@ -806,11 +1033,13 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
goto error;
}
/* Create an external solid */
- res = create_solid(dev, 1, 1, 1, 1, t_allowed, NULL, NULL, &media,
- solid_dont_get_temperature, NULL, &desc->d.t_boundary.mat_id);
+ res = create_solid(dev, 1, 1, 1, 1, prohibited, (mode & GREEN_MODE), 1,
+ NULL, NULL, &media, solid_dont_get_temperature, NULL,
+ &desc->d.t_boundary.mat_id);
if (res != RES_OK) goto error;
break;
case DESC_BOUND_F_FOR_SOLID:
+ fbound_count++;
ASSERT(desc->d.f_boundary.flux != NULL);
res = compile_expr_to_fn(&desc->d.f_boundary.te_flux,
desc->d.f_boundary.flux, 1, NULL);
@@ -820,42 +1049,48 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
goto error;
}
/* Create an external solid */
- res = create_solid(dev, 1, 1, 1, 1, t_allowed, NULL, NULL, &media, NULL,
- NULL, &desc->d.f_boundary.mat_id);
+ res = create_solid(dev, 1, 1, 1, 1, prohibited, (mode & GREEN_MODE), 1,
+ NULL, NULL, &media, solid_dont_get_temperature, NULL,
+ &desc->d.f_boundary.mat_id);
if (res != RES_OK) goto error;
break;
case DESC_SOLID_FLUID_CONNECT:
+ sfconnect_count++;
ASSERT(desc->d.sf_connect.specular_fraction >= 0
&& desc->d.sf_connect.specular_fraction <= 1);
ASSERT(desc->d.sf_connect.emissivity >= 0);
ASSERT(desc->d.sf_connect.hc >= 0);
break;
case DESC_MAT_SOLID:
+ smed_count++;
+ if (mode & GREEN_MODE) prohibited |= T_PROHIBITED | XYZ_PROHIBITED;
res = create_solid(dev,
stardis->descriptions[i].d.solid.lambda,
stardis->descriptions[i].d.solid.rho,
stardis->descriptions[i].d.solid.cp,
stardis->descriptions[i].d.solid.delta,
- t_allowed,
+ prohibited,
+ (mode & GREEN_MODE), 0,
stardis->descriptions[i].d.solid.Tinit,
stardis->descriptions[i].d.solid.power,
&media,
solid_get_tinit,
&stardis->descriptions[i].d.solid.has_power,
&desc->d.solid.solid_id);
- strcpy(desc->d.solid.name, desc->name);
if (res != RES_OK) goto error;
break;
case DESC_MAT_FLUID:
+ fmed_count++;
+ if (mode & GREEN_MODE) prohibited |= T_PROHIBITED | XYZ_PROHIBITED;
res = create_fluid(dev,
stardis->descriptions[i].d.fluid.rho,
stardis->descriptions[i].d.fluid.cp,
- t_allowed,
+ prohibited,
+ (mode & GREEN_MODE), 0,
stardis->descriptions[i].d.fluid.Tinit,
&media,
fluid_get_tinit,
&desc->d.fluid.fluid_id);
- strcpy(desc->d.fluid.name, desc->name);
if (res != RES_OK) goto error;
break;
default: FATAL("Invalid type.\n");
@@ -902,6 +1137,8 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
interface_props = sdis_data_get(data);
interface_props->temperature = NULL;
interface_props->flux = NULL;
+ interface_props->front_boundary_id = UINT_MAX;
+ interface_props->back_boundary_id = UINT_MAX;
if (front_defined) {
switch (stardis->descriptions[fd].type) {
@@ -957,6 +1194,9 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
goto error;
}
def_medium = front_defined ? front_med : back_med;
+ if (front_defined)
+ interface_props->front_boundary_id = cd;
+ else interface_props->back_boundary_id = cd;
}
switch (connect->type) {
case DESC_BOUND_H_FOR_FLUID:
@@ -1128,179 +1368,373 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
res = create_edge_file_if(scn, &stardis->allocator);
if (res != RES_OK) goto error;
- if (mode == IR_COMPUTE) {
- size_t width = (size_t)stardis->camera.img[0];
- size_t height = (size_t)stardis->camera.img[1];
-
- /* Setup the camera */
- SDIS(camera_create(dev, &cam));
- SDIS(camera_set_proj_ratio(cam, (double)width / (double)height));
- SDIS(camera_set_fov(cam, MDEG2RAD(stardis->camera.fov)));
- SDIS(camera_look_at(cam,
- stardis->camera.pos,
- stardis->camera.tgt,
- stardis->camera.up));
-
- /* Create the accum buffer */
- SDIS(accum_buffer_create(dev, width, height, &buf));
-
- /* Launch the simulation */
- res = sdis_solve_camera(scn,
- cam,
- stardis->camera.u.time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- width,
- height,
- (size_t)stardis->camera.spp,
- sdis_accum_buffer_write,
- buf);
- if (res != RES_OK) goto error;
+ if (mode & GREEN_MODE) {
+ struct sdis_green_function* green = NULL;
+ size_t ok_count, failed_count;
+ ASSERT((mode & PROBE_COMPUTE)
+ || (mode & BOUNDARY_COMPUTE) || (mode & MEDIUM_COMPUTE));
+
+ if (mode & PROBE_COMPUTE) {
+ double uv[2] = { 0,0 };
+ size_t iprim = SIZE_MAX;
+ /* Launch the probe simulation */
+ pos[0] = stardis->probe[0];
+ pos[1] = stardis->probe[1];
+ pos[2] = stardis->probe[2];
+
+ res = select_probe_type(scn,
+ stardis->geometry.triangle,
+ stardis->descriptions,
+ pos,
+ &iprim,
+ uv);
+ if (res != RES_OK) goto error;
- /* Write the image */
- dump_image(buf);
- }
- else if (mode == PROBE_COMPUTE) {
- double uv[2] = { 0,0 };
- size_t iprim = SIZE_MAX;
- /* Launch the probe simulation */
- pos[0] = stardis->probe[0];
- pos[1] = stardis->probe[1];
- pos[2] = stardis->probe[2];
- time[0] = time[1] = stardis->probe[3];
-
- res = select_probe_type(scn,
- stardis->geometry.triangle,
- stardis->descriptions,
- pos,
- &iprim,
- uv);
- if (res != RES_OK) goto error;
+ if (iprim == SIZE_MAX) {
+ res = sdis_solve_probe_green_function(scn,
+ stardis->N,
+ pos,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ } else {
+ res = sdis_solve_probe_boundary_green_function(scn,
+ stardis->N,
+ iprim,
+ uv,
+ SDIS_FRONT,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ }
+ }
+ else if (mode & BOUNDARY_COMPUTE) {
+ ASSERT(sa_size(stardis->boundary.primitives)
+ == sa_size(stardis->boundary.sides));
+ if (sa_size(stardis->boundary.sides) == 0) {
+ fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
+ stardis->solve_name);
+ }
- if (iprim == SIZE_MAX) {
- res = sdis_solve_probe(scn,
+ res = sdis_solve_boundary_green_function(scn,
stardis->N,
- pos,
- time,
+ stardis->boundary.primitives,
+ stardis->boundary.sides,
+ sa_size(stardis->boundary.sides),
stardis->scale_factor,
stardis->radiative_temp[0],
stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- } else {
- res = sdis_solve_probe_boundary(scn,
+ &green);
+ }
+ else if (mode & MEDIUM_COMPUTE) {
+ size_t sz, ii;
+ struct sdis_medium* medium = NULL;
+ /* Find medium */
+ sz = sa_size(stardis->descriptions);
+ FOR_EACH(ii, 0, sz) {
+ struct description* desc = stardis->descriptions + ii;
+ if (desc->type == DESC_MAT_SOLID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
+ }
+ else if (desc->type == DESC_MAT_FLUID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
+ }
+ }
+ if (medium == NULL) {
+ /* Not found */
+ fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
+ stardis->solve_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ res = sdis_solve_medium_green_function(scn,
stardis->N,
- iprim,
- uv,
- time,
- SDIS_FRONT,
+ medium,
stardis->scale_factor,
stardis->radiative_temp[0],
stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
+ &green);
}
- }
- else if (mode == MEDIUM_COMPUTE) {
- size_t sz, ii;
- struct sdis_medium* medium = NULL;
- /* Find medium */
- sz = sa_size(stardis->descriptions);
- FOR_EACH(ii, 0, sz) {
- struct description* desc = stardis->descriptions + ii;
- if (desc->type == DESC_MAT_SOLID) {
- if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
+ if (res != RES_OK) goto error;
+
+ /* Output the green function */
+ CHK(sdis_green_function_get_invalid_paths_count(green, &failed_count) == RES_OK);
+ CHK(sdis_green_function_get_paths_count(green, &ok_count) == RES_OK);
+
+ /* Output counts */
+ printf("---BEGIN GREEN---\n");
+ printf("# #solids #fluids #t_boundaries #h_boundaries #f_boundaries #ok #failures\n");
+ printf("%zu %zu %zu %zu %zu %zu %zu\n",
+ smed_count, fmed_count,
+ tbound_count, hbound_count, fbound_count,
+ ok_count, failed_count);
+
+ /* List Media */
+ if (smed_count) {
+ printf("# Solids\n");
+ printf("# ID Name lambda rho cp power\n");
+ FOR_EACH(i, 0, (unsigned)sa_size(stardis->descriptions)) {
+ struct description* desc = stardis->descriptions + i;
+ struct mat_solid* sl;
+ if (desc->type != DESC_MAT_SOLID) continue;
+ sl = &desc->d.solid;
+ printf("%u\t%s\t%g\t%g\t%g\t%s\n",
+ i, sl->name, sl->lambda, sl->rho, sl->cp, (sl->has_power ? sl->power : "0"));
}
- else if (desc->type == DESC_MAT_FLUID) {
- if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
+ }
+ if (fmed_count) {
+ printf("# Fluids\n");
+ printf("# ID Name rho cp\n");
+ FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
+ struct description* desc = stardis->descriptions + i;
+ struct mat_fluid* fl;
+ fl = &desc->d.fluid;
+ printf("%u\t%s\t%g\t%g\n", i, fl->name, fl->rho, fl->cp);
}
}
- if (medium == NULL) {
- /* Not found */
- fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
- stardis->solve_name);
- res = RES_BAD_ARG;
- goto error;
+
+ /* List Boundaries */
+ if (tbound_count) {
+ printf("# T Boundaries\n");
+ printf("# ID Name temperature\n");
+ FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
+ struct description* desc = stardis->descriptions + i;
+ struct t_boundary* bd;
+ if (desc->type != DESC_BOUND_T_FOR_SOLID
+ && desc->type != DESC_BOUND_T_FOR_FLUID) continue;
+ bd = &desc->d.t_boundary;
+ printf("%u\t%s\t%s\n", i, bd->name, bd->T);
+ }
}
- time[0] = time[1] = stardis->probe[3];
- res = sdis_solve_medium(scn,
- stardis->N,
- medium,
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- } else if (mode == BOUNDARY_COMPUTE) {
- ASSERT(sa_size(stardis->boundary.primitives)
- == sa_size(stardis->boundary.sides));
- if (sa_size(stardis->boundary.sides) == 0) {
- fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
- stardis->solve_name);
+ if (hbound_count) {
+ printf("# H Boundaries\n");
+ printf("# ID Name emissivity specular_fraction hc hc_max T_env\n");
+ FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
+ struct description* desc = stardis->descriptions + i;
+ struct h_boundary* bd;
+ if (desc->type != DESC_BOUND_H_FOR_SOLID
+ && desc->type != DESC_BOUND_H_FOR_FLUID) continue;
+ bd = &desc->d.h_boundary;
+ printf("%u\t%s\t%g\t%g\t%g\t%g\t%s\n",
+ i, bd->name, (bd->has_emissivity ? bd->emissivity : 0),
+ (bd->has_emissivity ? bd->specular_fraction : 0),
+ (bd->has_hc ? bd->hc : 0), (bd->has_hc ? bd->hc_max : 0), bd->T);
+ }
}
+ if (fbound_count) {
+ printf("# F Boundaries\n");
+ printf("# ID Name hc hc_max\n");
+ FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
+ struct description* desc = stardis->descriptions + i;
+ struct f_boundary* bd;
+ if (desc->type != DESC_BOUND_F_FOR_SOLID) continue;
+ bd = &desc->d.f_boundary;
+ printf("%u\t%s\t%g\t%g\n",
+ i, bd->name, (bd->has_hc ? bd->hc : 0), (bd->has_hc ? bd->hc_max : 0));
+ }
+ }
+
+ /* Radiative Temperatures */
+ printf("# Radiative Temperatures\n");
+ printf("# ID Rad_Temp Lin_Temp\n");
+ printf("%u\t%g\t%g\n", (unsigned)sa_size(stardis->descriptions),
+ stardis->radiative_temp[0], stardis->radiative_temp[1]);
- time[0] = time[1] = stardis->probe[3];
- res = sdis_solve_boundary(scn,
- stardis->N,
- stardis->boundary.primitives,
- stardis->boundary.sides,
- sa_size(stardis->boundary.sides),
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
+ printf("# Samples\n");
+ printf("# end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n\n");
+ printf("# end = end_type end_id; end_type = T | H | X | R | F | S\n");
+ printf("# power_term_i = power_type_i power_id_i factor_i\n");
+ printf("# flux_term_i = flux_id_i factor_i\n");
+ sdis_green_function_for_each_path(green, print_sample, stardis->descriptions);
+ printf("---END GREEN---\n");
+
+ CHK(sdis_green_function_ref_put(green) == RES_OK);
} else {
- CHK(0);
- }
- if (res != RES_OK) goto error;
- if (mode == PROBE_COMPUTE
- || mode == MEDIUM_COMPUTE
- || mode == BOUNDARY_COMPUTE) {
- /* Fetch the estimation data */
- SDIS(estimator_get_temperature(estimator, &temperature));
- SDIS(estimator_get_failure_count(estimator, &nfailures));
-
- /* Print the results */
- switch (mode) {
- case PROBE_COMPUTE:
- printf("Temperature at t=[%g, %g, %g, %g] = %g +/- %g\n",
- pos[0], pos[1], pos[2], time[0],
- temperature.E, /* Expected value */
- temperature.SE); /* Standard error */
- break;
- case MEDIUM_COMPUTE:
- printf("Temperature in medium %s at t=%g = %g +/- %g\n",
- stardis->solve_name, time[0],
- temperature.E, /* Expected value */
- temperature.SE); /* Standard error */
- break;
- case BOUNDARY_COMPUTE:
- printf("Temperature at boundary %s at t=%g = %g +/- %g\n",
- stardis->solve_name, time[0],
- temperature.E, /* Expected value */
- temperature.SE); /* Standard error */
- break;
- default: FATAL("Invalid mode.");
+ if (mode & IR_COMPUTE) {
+ size_t width = (size_t)stardis->camera.img[0];
+ size_t height = (size_t)stardis->camera.img[1];
+
+ /* Setup the camera */
+ SDIS(camera_create(dev, &cam));
+ SDIS(camera_set_proj_ratio(cam, (double)width / (double)height));
+ SDIS(camera_set_fov(cam, MDEG2RAD(stardis->camera.fov)));
+ SDIS(camera_look_at(cam,
+ stardis->camera.pos,
+ stardis->camera.tgt,
+ stardis->camera.up));
+
+ /* Create the accum buffer */
+ SDIS(accum_buffer_create(dev, width, height, &buf));
+
+ /* Launch the simulation */
+ res = sdis_solve_camera(scn,
+ cam,
+ stardis->camera.u.time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ width,
+ height,
+ (size_t)stardis->camera.spp,
+ sdis_accum_buffer_write,
+ buf);
+ if (res != RES_OK) goto error;
+
+ /* Write the image */
+ dump_image(buf);
+ }
+ else if (mode & PROBE_COMPUTE) {
+ double uv[2] = { 0,0 };
+ size_t iprim = SIZE_MAX;
+ /* Launch the probe simulation */
+ pos[0] = stardis->probe[0];
+ pos[1] = stardis->probe[1];
+ pos[2] = stardis->probe[2];
+ time[0] = time[1] = stardis->probe[3];
+
+ res = select_probe_type(scn,
+ stardis->geometry.triangle,
+ stardis->descriptions,
+ pos,
+ &iprim,
+ uv);
+ if (res != RES_OK) goto error;
+
+ if (iprim == SIZE_MAX) {
+ res = sdis_solve_probe(scn,
+ stardis->N,
+ pos,
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ } else {
+ res = sdis_solve_probe_boundary(scn,
+ stardis->N,
+ iprim,
+ uv,
+ time,
+ SDIS_FRONT,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ }
+ }
+ else if (mode & MEDIUM_COMPUTE) {
+ size_t sz, ii;
+ struct sdis_medium* medium = NULL;
+ /* Find medium */
+ sz = sa_size(stardis->descriptions);
+ FOR_EACH(ii, 0, sz) {
+ struct description* desc = stardis->descriptions + ii;
+ if (desc->type == DESC_MAT_SOLID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
+ }
+ else if (desc->type == DESC_MAT_FLUID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
+ }
+ }
+ if (medium == NULL) {
+ /* Not found */
+ fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
+ stardis->solve_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ time[0] = time[1] = stardis->probe[3];
+ res = sdis_solve_medium(scn,
+ stardis->N,
+ medium,
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
}
- printf("#failures: %lu/%lu\n",
- (unsigned long)nfailures,
- (unsigned long)stardis->N);
+ else if (mode & BOUNDARY_COMPUTE) {
+ ASSERT(sa_size(stardis->boundary.primitives)
+ == sa_size(stardis->boundary.sides));
+ if (sa_size(stardis->boundary.sides) == 0) {
+ fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
+ stardis->solve_name);
+ }
- /* Dump paths according to user settings */
- sdis_estimator_for_each_path(estimator, dump_path, stdout);
+ time[0] = time[1] = stardis->probe[3];
+ res = sdis_solve_boundary(scn,
+ stardis->N,
+ stardis->boundary.primitives,
+ stardis->boundary.sides,
+ sa_size(stardis->boundary.sides),
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ } else {
+ CHK(0);
+ }
+ if (res != RES_OK) goto error;
+ if (mode & PROBE_COMPUTE
+ || mode & MEDIUM_COMPUTE
+ || mode & BOUNDARY_COMPUTE) {
+ /* Fetch the estimation data */
+ SDIS(estimator_get_temperature(estimator, &temperature));
+ SDIS(estimator_get_failure_count(estimator, &nfailures));
+
+ /* Print the results */
+ switch (mode) {
+ case PROBE_COMPUTE:
+ printf("Temperature at t=[%g, %g, %g, %g] = %g +/- %g\n",
+ pos[0], pos[1], pos[2], time[0],
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ case MEDIUM_COMPUTE:
+ printf("Temperature in medium %s at t=%g = %g +/- %g\n",
+ stardis->solve_name, time[0],
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ case BOUNDARY_COMPUTE:
+ printf("Temperature at boundary %s at t=%g = %g +/- %g\n",
+ stardis->solve_name, time[0],
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ default: FATAL("Invalid mode.");
+ }
+ printf("#failures: %lu/%lu\n",
+ (unsigned long)nfailures,
+ (unsigned long)stardis->N);
+
+ /* Dump paths according to user settings */
+ sdis_estimator_for_each_path(estimator, dump_path, stdout);
+ }
}
end:
