stardis

Perform coupled heat transfer calculations
git clone git://git.meso-star.fr/stardis.git
Log | Files | Refs | README | LICENSE

commit 9dc4d7f2f779c36fd5eebd36c1bea499d849d411
parent c9519be529e1215d8faa76973067745451e01b07
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Wed, 30 Jun 2021 14:58:37 +0200

Merge branch 'release_0.6'

Diffstat:
MREADME.md | 24+++++++++++++++++++-----
Mcmake/CMakeLists.txt | 28+++++++++++++++-------------
Mcmake/doc/CMakeLists.txt | 2+-
Mdoc/stardis-input.5.txt | 18+++++++++++++++---
Mdoc/stardis-output.5.txt | 20+++++++++++---------
Mdoc/stardis.1.txt.in | 79++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/stardis-app.c | 31++++++++++++++++++++++++++-----
Msrc/stardis-app.h | 92++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----
Msrc/stardis-compute.c | 418++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------
Msrc/stardis-compute.h | 6++++--
Msrc/stardis-fluid.c | 2+-
Msrc/stardis-fluid.h | 2+-
Msrc/stardis-intface.c | 33++++++++++++++++++++++++++++++---
Msrc/stardis-intface.h | 5++++-
Msrc/stardis-main.c | 10+++++++---
Msrc/stardis-output.c | 391+++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------
Msrc/stardis-output.h | 23++++++++++++++++++++++-
Msrc/stardis-parsing.c | 439+++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------
Msrc/stardis-parsing.h | 13++++++++++---
Msrc/stardis-solid.c | 2+-
Msrc/stardis-solid.h | 2+-
21 files changed, 1240 insertions(+), 400 deletions(-)

diff --git a/README.md b/README.md @@ -25,14 +25,21 @@ as well as on the [OpenMP](http://www.openmp.org) 2.0 specification to parallelize its computations. -First ensure that CMake and a C89 compiler that implements the OpenMP 2.0 -specification are installed on your system. Then install the RCMake package as -well as all the aforementioned prerequisites. Finally generate the project from +First ensure that CMake and a C compiler are installed on your system. +Then install the RCMake package as well as all the aforementioned +prerequisites. Finally generate the project from the `cmake/CMakeLists.txt` file by appending to the `CMAKE_PREFIX_PATH` variable the install directories of its dependencies. ## Release notes +### Version 0.6 + +- Add thermal contact resistances. +- Add serialization for random generator's state. +- Bugfixes in arguments parsing. +- Fix Green output file padding. + ### Version 0.5.1 - Fix a memleak @@ -84,4 +91,12 @@ Add radiative transfer computations. To achieve this, media gain 2 new parameter ### Version 0.1 - Allow probe computations on conductive-only thermal systems. -- Allow Dirichlet and h.dT boundary conditions. -\ No newline at end of file +- Allow Dirichlet and h.dT boundary conditions. + +## License + +Copyright (C) 2018-2021 |Meso|Star> (<contact@meso-star.com>). Stardis is free +software released under the GPL v3+ license: GNU GPL version 3 or later. You +are welcome to redistribute it under certain conditions; refer to the COPYING +file for details. + diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) 2018-2020 |Meso|Star> +# Copyright (C) 2018-2021 |Meso|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 @@ -36,7 +36,7 @@ set_property(CACHE STARDIS_DOC PROPERTY STRINGS # Generate files ############################################################################### set(STARDIS_ARGS_DEFAULT_AMBIENT_TEMP "300") -set(STARDIS_ARGS_DEFAULT_COMPUTE_TIME "INF, INF") +set(STARDIS_ARGS_DEFAULT_COMPUTE_TIME "INF") set(STARDIS_ARGS_DEFAULT_RENDERING_FOV "70") # degrees set(STARDIS_ARGS_DEFAULT_RENDERING_IMG_HEIGHT "480") set(STARDIS_ARGS_DEFAULT_RENDERING_IMG_WIDTH "640") @@ -55,8 +55,8 @@ configure_file(${SDIS_SOURCE_DIR}/../doc/stardis.1.txt.in ${CMAKE_CURRENT_BINARY_DIR}/doc/stardis.1.txt @ONLY) set(SDIS_VERSION_MAJOR 0) -set(SDIS_VERSION_MINOR 5) -set(SDIS_VERSION_PATCH 1) +set(SDIS_VERSION_MINOR 6) +set(SDIS_VERSION_PATCH 0) set(SDIS_VERSION ${SDIS_VERSION_MAJOR}.${SDIS_VERSION_MINOR}.${SDIS_VERSION_PATCH}) configure_file(${SDIS_SOURCE_DIR}/stardis-default.h.in @@ -68,13 +68,14 @@ configure_file(${SDIS_SOURCE_DIR}/stardis-version.h.in ############################################################################### # Check dependencies ############################################################################### -find_package(RCMake 0.4 REQUIRED) -find_package(RSys 0.11 REQUIRED) +find_package(RCMake 0.4.1 REQUIRED) +find_package(RSys 0.12 REQUIRED) find_package(StarGeom3D 0.1 REQUIRED) -find_package(Star3D 0.7.3 REQUIRED) -find_package(StarEnc3D 0.4.2 REQUIRED) -find_package(Stardis 0.11 REQUIRED) -find_package(StarSTL 0.3 REQUIRED) +find_package(Star3D 0.8 REQUIRED) +find_package(StarEnc3D 0.5.3 REQUIRED) +find_package(Stardis 0.12 REQUIRED) +find_package(StarSTL 0.3.3 REQUIRED) +find_package(StarSP 0.8.1 REQUIRED) if(MSVC) find_package(MuslGetopt REQUIRED) endif() @@ -86,6 +87,7 @@ include_directories( ${StarEnc3D_INCLUDE_DIR} ${Stardis_INCLUDE_DIR} ${StarSTL_INCLUDE_DIR} + ${StarSP_INCLUDE_DIR} ${CMAKE_CURRENT_BINARY_DIR}) if(MSVC) include_directories(${MuslGetopt_INCLUDE_DIR}) @@ -97,11 +99,11 @@ include(rcmake_runtime) if(CMAKE_COMPILER_IS_GNUCC) rcmake_append_runtime_dirs(_runtime_dirs - RSys Stardis Star3D StarGeom3D StarEnc3D StarSTL) + RSys Stardis Star3D StarGeom3D StarEnc3D StarSTL StarSP) endif() if(MSVC) rcmake_append_runtime_dirs(_runtime_dirs - RSys MuslGetopt Stardis Star3D StarGeom3D StarEnc3D StarSTL) + RSys MuslGetopt Stardis Star3D StarGeom3D StarEnc3D StarSTL StarSP) endif() ############################################################################### @@ -156,7 +158,7 @@ set_target_properties(stardis PROPERTIES VERSION ${SDIS_VERSION}) target_link_libraries(stardis - Stardis Star3D StarGeom3D StarEnc3D StarSTL RSys ${GETOPT_LIB} ${MATH_LIB}) + Stardis Star3D StarGeom3D StarEnc3D StarSTL StarSP RSys ${GETOPT_LIB} ${MATH_LIB}) ############################################################################### # Define output & install directories diff --git a/cmake/doc/CMakeLists.txt b/cmake/doc/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) 2018-2020 |Meso|Star> +# Copyright (C) 2018-2021 |Meso|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 diff --git a/doc/stardis-input.5.txt b/doc/stardis-input.5.txt @@ -99,6 +99,7 @@ _______ | <f-bound-for-solid> <media-connection> ::= <solid-fluid-connect> + | <solid-solid-connect> <comment> ::= "#" Any text introduced by the # character @@ -132,9 +133,13 @@ _______ <solid-fluid-connect> ::= "SOLID_FLUID_CONNECTION" <bound-name> <emissivity> \ <specular-fraction> <hc> <triangles> +<solid-solid-connect> ::= "SOLID_SOLID_CONNECTION" <bound-name> \ + <contact-resistance> <triangles> + ------------------------------------- -<medium-name> ::= STRING # no space allowed +<medium-name> ::= STRING # no space allowed, must not be parsable as a + # number, including INF and others <lambda> ::= REAL # conductivity in W/(m.K); in ]0, INF) @@ -167,6 +172,8 @@ _______ <hc> ::= REAL # in W/(m2.K); in [0, INF) +<contact-resistance> ::= REAL # in m2.K/W in [0, INF) + <flux> ::= REAL # in W/m2; in (-INF , INF) <triangles> ::= <file-name> [ <triangles> ] @@ -192,9 +199,14 @@ NAMES Names, either file names, medium names or boundary names, are a sequence of one or ore ASCII characters, including numbers and special characters like *.* *_* *-* as one may consider using in standard file names, *without any -spacing* either escaped or not. Names are case-sensitive anf two different +spacing* either escaped or not. Names are case-sensitive and two different description lines, either in the same description file or from different -description files, cannot use the same name. +description files, cannot use the same name. Additionaly, medium and boundary +names cannot be parsable as a number, nor be one of the few keywords defined +by the present grammar (AUTO, BACK, BOTH, FLUID, FRONT, F_BOUNDARY_FOR_SOLID, +H_BOUNDARY_FOR_FLUID, H_BOUNDARY_FOR_SOLID, SCALE, SOLID, +SOLID_FLUID_CONNECTION, T_BOUNDARY_FOR_FLUID, T_BOUNDARY_FOR_SOLID, UNKNOWN) or +their lowercase counterparts. EXAMPLES -------- diff --git a/doc/stardis-output.5.txt b/doc/stardis-output.5.txt @@ -139,34 +139,36 @@ _______ <failures-report> ::= <error-count> <success-count> <ext-probe-temp> ::= "Temperature at" <ext-probe-position> <ext-time> \ - <ext-MC-estimate> + <ext-temp-estimate> <ext-failures-report> <ext-mean-temp> ::= "Temperature at boundary" <file-name> <ext-time> \ - <ext-MC-estimate> + <ext-temp-estimate> <ext-failures-report> <ext-mean-flux> ::= "Temperature at boundary" <file-name> <ext-time> \ - <ext-MC-estimate> + <ext-temp-estimate> "Convective flux at boundary " <file-name> \ - <ext-time> <ext-MC-estimate> + <ext-time> <ext-flux-estimate> "Radiative flux at boundary" <file-name> <ext-time> \ - <ext-MC-estimate> + <ext-flux-estimate> "Imposed flux at boundary" <file-name> <ext-time> \ - <ext-MC-estimate> + <ext-flux-estimate> "Total flux at boundary" <file-name> <ext-time> \ - <ext-MC-estimate> + <ext-flux-estimate> <ext-failures-report> <ext-medium-temp> ::= "Temperature in medium" <medium-name> <ext-time> \ - <ext-MC-estimate> + <ext-temp-estimate> <ext-failures-report> <ext-probe-position> ::= "[" <x-value> "," <y-value> "," <z-value> "]" <ext-time> ::= "at t=" <time-value> "=" -<ext-MC-estimate> ::= <expected-value> "+/-" <standard-deviation> +<ext-temp-estimate> ::= <expected-value> "K +/-" <standard-deviation> + +<ext-flux-estimate> ::= <expected-value> "W +/-" <standard-deviation> <ext-failures-report> ::= "#failures:" <error-count> "/" <success-count> diff --git a/doc/stardis.1.txt.in b/doc/stardis.1.txt.in @@ -48,7 +48,7 @@ the evaluation of the propagator (a.k.a the *Green function*). The propagator is of great value for thermicist engineers as it gives some crucial information to analyse heat transfers in the system. It helps engineers answer questions like _"Where from does the heat come at this location?"_. -Propagators seamlessly agregate all the provided geometrical and physical +Propagators seamlessly aggregate all the provided geometrical and physical information on the system in an unbiased and very-fast statistical model. *stardis*(1) also provides two additional functionalities: converting the @@ -60,12 +60,10 @@ to radiative transfer physics (Delatorre [1]) combined with conduction's statistical formulation (Kac [2] and Muller [3]). Thanks to recent advances in 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. While this -capability is part of the StarEngine Star3D library, it is internally powered -by IntelĀ® Rendering Framework: Embree. +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 mecanisms (conductive-convective-radiative) +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. @@ -91,57 +89,64 @@ MANDATORY OPTIONS EXCLUSIVE OPTIONS ----------------- -*-p* _x,y,z[,time-range]_:: +*-p* _x,y,z[,time[,time]]_:: Compute the temperature at the given probe at a given time. By default the - compute time range is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The probe must + compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The probe must be in a medium. The probe coordinates must be in the same system as the geometry. -*-P* _x,y,z[,time-range]_:: +*-P* _x,y,z[,time[,time]][:side_indicator]_:: Compute the temperature at the given probe on an interface at a given time. - By default the compute time range is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The + If the probe is on an interface where a thermal contact resistance is + defined, it is mandatory to provide a side indicator (either FRONT, BACK, or + a medium name), as the temperature differs between the two sides. + By default the compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The probe is supposed to be on an interface and is moved to the closest point of - the closest interface before the computation start. The probe coordinates + the closest interface before the computation starts. The probe coordinates must be in the same system as the geometry. -*-m* _medium_name[,time-range]_:: +*-m* _medium_name[,time[,time]]_:: Compute the mean temperature in a given medium at a given time. The medium name must be part of the system description. By default the compute time - range is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The medium does not need to be - connex. + is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The medium region does not need + to be connex. -*-s* _file[,time-range]_:: +*-s* _file[,time[,time]]_:: Compute the mean temperature on a given 2D region at a given time, the region being defined as the front sides of the triangles in the provided *STL* file. - By default the compute time range is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. + By default the compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. These triangles are not added to the geometry, but must be part of it. The region does not need to be connex. -*-S* _file[,time-range]_:: +*-S* _file[,time[,time]]_:: Compute the by-triangle mean temperature on a given 2D region at a given time, the region being defined as the front sides of the triangles in the - provided *VTK* file. These triangles are not added to the geometry, but must - be part of it. By default the compute time range is + provided *STL* file. 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 region does not need to be connex. -*-F* _file[,time-range]_:: +*-F* _file[,time[,time]]_:: Compute the mean flux on a given 2D region at a given time, the region - being defined as the front sides of the triangles in the provided *VTK* file. + being defined as the front sides of the triangles in the provided *STL* file. These triangles are not added to the geometry, but must be part of it. Flux is accounted positive when going from the front side to the back side, at a - single-triangle level. By default the compute time range is + single-triangle level. By default the compute time is @STARDIS_ARGS_DEFAULT_COMPUTE_TIME@. The region does not need to be connex, but it can currently only include geometry appearing in description lines starting with H_BOUNDARY_FOR_SOLID, H_BOUNDARY_FOR_FLUID, F_BOUNDARY_FOR_SOLID or SOLID_FLUID_CONNECTION (see *stardis-input*(5)). *-R* [__sub-option__:...]:: - Render an infrared image of the system through a pinhole camera and write it - to _standard output_. One can use all-default sub-options by - simply providing the colon character (*:*) alone as an argument. Please note - that the camera position must be outside the geometry or in a fluid. + Render an infrared image of the system through a pinhole camera. One can use + all-default sub-options by simply providing the colon character (*:*) alone + as an argument. Please note that the camera position must be outside the + geometry or in a fluid. Available sub-options are: + **file=**_output_file_;; + File name to use to write the infrared image to. If no file name is + provided, the result is written to _standard output_. + **fmt=**_image_file_format_;; Format of the image file in output. Can be *VTK* or *HT*. Default _image_file_format_ is @STARDIS_ARGS_DEFAULT_RENDERING_OUTPUT_FILE_FMT@. @@ -164,8 +169,8 @@ EXCLUSIVE OPTIONS Number of samples per pixel. By default, use @STARDIS_ARGS_DEFAULT_RENDERING_SPP@ samples per pixel. - **t=**_time-range_;; - Rendering time range. By default _time-range_ is + **t=**_time[,time]_;; + Rendering time. By default the rendering time is @STARDIS_ARGS_DEFAULT_RENDERING_TIME@. **tgt=**_x_**,**_y_**,**_z_;; @@ -202,7 +207,7 @@ options that write to _standard output_ (-g, -h, -R, -v). rank starting at index 00000000, and possibly followed by *_err* for failure paths: prefix00000000.vtk, prefix00000001_err.vtk, ... + -This option can only be used in conjuction with options that compute a +This option can only be used in conjunction with options that compute a result (-F, -m, -P, -p, -R, -S, -s) and cannot be used in conjunction with options -g or -G. @@ -254,6 +259,16 @@ different temperature, flux or volumic power values. and informative messages). All the messages are written to _standard error_. Default verbosity *level* is @STARDIS_ARGS_DEFAULT_VERBOSE_LEVEL@. +*-x* _file_name_:: + Read the provided file and use its content to initialize the random + generator's internal state. Used in conjunction with the *-X* option, this + can be used to ensure statistical independence between subsequent + computations. + +*-X* _file_name_:: + Write the random generator's internal state, as it is at the end of the + computation, to the provided file. + EXAMPLES -------- Preprocess the system as described in *scene 5.txt* when intending to compute @@ -284,6 +299,12 @@ time range. The system is read from the file *model.txt*: $ stardis -M model.txt -p 1,2.5,0,50,5000 +Compute 3 probe temperatures, ensuring statistical independence: + + $ stardis -M model.txt -p 1,1.5,0,50,5000 -Xstate1 + $ stardis -M model.txt -p 1,2.5,0,50,5000 -xstate1 -Xstate2 + $ stardis -M model.txt -p 1,3.5,0,50,5000 -xstate2 + Render the system as described in *scene.txt* with default settings: $ stardis -M scene.txt -R : @@ -299,7 +320,7 @@ post-processed using *htpp*(1) with default settings to obtain a png file. -D error,err_path_ > img.ht $ htpp -o img.pgn -v -m default img.ht -Compute the Green fonction that computes the temperature at the probe point +Compute the Green function that computes the temperature at the probe point *0, 0, 0* at steady state. The system is read from the file *model.txt* and the Green function is written to the *probe.green file* and the heat paths' ends are written to the *probe_ends.csv* file: diff --git a/src/stardis-app.c b/src/stardis-app.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -220,15 +220,18 @@ stardis_init stardis->senc3d_scn = NULL; stardis->mode = args->mode; stardis->counts = COUNTS_NULL; - init_camera(&stardis->camera); + init_camera(stardis->allocator, &stardis->camera); str_init(stardis->allocator, &stardis->solve_name); str_init(stardis->allocator, &stardis->paths_filename); str_init(stardis->allocator, &stardis->bin_green_filename); str_init(stardis->allocator, &stardis->end_paths_filename); + str_init(stardis->allocator, &stardis->rndgen_state_in_filename); + str_init(stardis->allocator, &stardis->rndgen_state_out_filename); str_init(stardis->allocator, &stardis->chunks_prefix); darray_size_t_init(stardis->allocator, &stardis->compute_surface.primitives); darray_sides_init(stardis->allocator, &stardis->compute_surface.sides); darray_uint_init(stardis->allocator, &stardis->compute_surface.err_triangles); + stardis->compute_surface.area = 0; stardis->samples = args->samples; stardis->nthreads = args->nthreads; stardis->scale_factor = -1; /* invalid value */ @@ -260,6 +263,9 @@ 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 & SURFACE_COMPUTE_MODES) { ERR(str_set(&stardis->solve_name, args->solve_filename)); } @@ -309,6 +315,8 @@ stardis_init goto error; } } + logger_print(stardis->logger, LOG_OUTPUT, "Compute surface area is %g m2\n", + stardis->compute_surface.area); } /* Create enclosures */ @@ -356,6 +364,14 @@ stardis_init if(args->end_paths_filename) { ERR(str_set(&stardis->end_paths_filename, args->end_paths_filename)); } + if(args->rndgen_state_in_filename) { + ERR(str_set(&stardis->rndgen_state_in_filename, + args->rndgen_state_in_filename)); + } + if(args->rndgen_state_out_filename) { + ERR(str_set(&stardis->rndgen_state_out_filename, + args->rndgen_state_out_filename)); + } } /* If computation is on a volume, check medium is known */ @@ -372,9 +388,6 @@ stardis_init struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; ASSERT(darray_interface_ptrs_size_get(&stardis->geometry.interf_bytrg) == tcount); - /* Can release enclosures as they are no longer needed if compute */ - SENC3D(scene_ref_put(stardis->senc3d_scn)); - stardis->senc3d_scn = NULL; scn_args.get_indices = sg3d_sdisXd_geometry_get_indices; scn_args.get_interface = sg3d_sdisXd_geometry_get_interface; scn_args.get_position = sg3d_sdisXd_geometry_get_position; @@ -427,6 +440,7 @@ stardis_release SDIS(medium_ref_put(darray_media_ptr_data_get(&stardis->media)[i])); } darray_media_ptr_release(&stardis->media); + release_camera(&stardis->camera); } res_T @@ -572,6 +586,13 @@ validate_properties goto end; } break; + case DESC_SOLID_SOLID_CONNECT: + if(solid_count != 2) { + /*if(soli_count == 1 && fluid_count == 1)*/ + /**properties_conflict_status = SSCONNECT_BETWEEN_SOLID_AND_FLUID;*/ + goto end; + } + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } diff --git a/src/stardis-app.h b/src/stardis-app.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -106,6 +106,7 @@ enum description_type { DESC_BOUND_T_FOR_SOLID, DESC_BOUND_F_FOR_SOLID, DESC_SOLID_FLUID_CONNECT, + DESC_SOLID_SOLID_CONNECT, DESCRIPTION_TYPE_COUNT__, DESC_OUTSIDE }; @@ -415,6 +416,58 @@ cp_release_sf_connect return cp_sf_connect(dst, src); } +struct solid_solid_connect { + struct str name; + double tcr; + unsigned connection_id; +}; + +static FINLINE void +release_ss_connect(struct solid_solid_connect* connect) +{ + str_release(&connect->name); +} + +static FINLINE void +init_ss(struct mem_allocator* allocator, struct solid_solid_connect* dst) +{ + str_init(allocator, &dst->name); + dst->tcr = 0; + dst->connection_id = UINT_MAX; +} + +static res_T +str_print_ss_connect + (struct str* str, + const struct solid_solid_connect* c) +{ + res_T res = RES_OK; + ASSERT(str && c); + STR_APPEND_PRINTF(str, "Solid-Solid connection '%s':", ARG1( str_cget(&c->name) ) ); + STR_APPEND_PRINTF(str, " contact resistance=%g", ARG1( c->tcr ) ); +end: + return res; +error: + goto end; +} + +static FINLINE res_T +cp_ss_connect + (struct solid_solid_connect* dst, const struct solid_solid_connect* src) +{ + dst->connection_id = src->connection_id; + dst->tcr = src->tcr; + return str_copy(&dst->name, &src->name); +} + +static FINLINE res_T +cp_release_ss_connect + (struct solid_solid_connect* dst, struct solid_solid_connect* src) +{ + return cp_ss_connect(dst, src); +} + + struct description { enum description_type type; union { @@ -424,6 +477,7 @@ struct description { struct f_boundary f_boundary; struct h_boundary h_boundary; struct solid_fluid_connect sf_connect; + struct solid_solid_connect ss_connect; } d; }; @@ -460,6 +514,9 @@ release_description(struct description* desc) case DESC_SOLID_FLUID_CONNECT: release_sf_connect(&desc->d.sf_connect); break; + case DESC_SOLID_SOLID_CONNECT: + release_ss_connect(&desc->d.ss_connect); + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } @@ -496,6 +553,9 @@ str_print_description case DESC_SOLID_FLUID_CONNECT: ERR(str_print_sf_connect(str, &desc->d.sf_connect)); break; + case DESC_SOLID_SOLID_CONNECT: + ERR(str_print_ss_connect(str, &desc->d.ss_connect)); + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } @@ -525,6 +585,8 @@ get_description_name return &desc->d.f_boundary.name; case DESC_SOLID_FLUID_CONNECT: return &desc->d.sf_connect.name; + case DESC_SOLID_SOLID_CONNECT: + return &desc->d.ss_connect.name; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } @@ -582,7 +644,14 @@ cp_description } ERR(cp_sf_connect(&dst->d.sf_connect, &src->d.sf_connect)); break; - default: + case DESC_SOLID_SOLID_CONNECT: + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + init_ss(src->d.ss_connect.name.allocator, &dst->d.ss_connect); + } + ERR(cp_ss_connect(&dst->d.ss_connect, &src->d.ss_connect)); + break; + default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } end: @@ -616,6 +685,7 @@ description_get_medium_id *id = desc->d.f_boundary.mat_id; return; case DESC_SOLID_FLUID_CONNECT: /* No medium linked to SF */ + case DESC_SOLID_SOLID_CONNECT: /* No medium linked to SS */ default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } @@ -636,10 +706,12 @@ struct camera { unsigned spp; unsigned img_width, img_height; int auto_look_at; + struct str file_name; }; static INLINE void -init_camera(struct camera* cam) { +init_camera(struct mem_allocator* alloc, struct camera* cam) { + ASSERT(alloc && cam); d3(cam->pos, STARDIS_DEFAULT_RENDERING_POS); d3(cam->tgt, STARDIS_DEFAULT_RENDERING_TGT); d3(cam->up, STARDIS_DEFAULT_RENDERING_UP); @@ -650,6 +722,13 @@ init_camera(struct camera* cam) { cam->img_height = STARDIS_DEFAULT_RENDERING_IMG_HEIGHT; d2(cam->time_range, STARDIS_DEFAULT_RENDERING_TIME); cam->auto_look_at = 1; + str_init(alloc, &cam->file_name); +} + +static INLINE void +release_camera(struct camera* cam) { + ASSERT(cam); + str_release(&cam->file_name); } static INLINE void @@ -690,10 +769,10 @@ log_prt_fn(const char* msg, void* ctx) struct counts { unsigned smed_count, fmed_count, tbound_count, hbound_count, - fbound_count, sfconnect_count; + fbound_count, sfconnect_count, ssconnect_count; }; #define COUNTS_NULL__ {\ - 0, 0, 0, 0, 0, 0\ + 0, 0, 0, 0, 0, 0, 0\ } /* Type to store the primitives of a compute region */ @@ -712,6 +791,7 @@ struct compute_surface { struct darray_size_t primitives; struct darray_sides sides; struct darray_uint err_triangles; + double area; /* in m2, regardless of scale factor */ }; struct stardis { @@ -734,6 +814,8 @@ struct stardis { struct str paths_filename; struct str bin_green_filename; struct str end_paths_filename; + struct str rndgen_state_in_filename; + struct str rndgen_state_out_filename; struct str chunks_prefix; int mode; size_t samples; diff --git a/src/stardis-compute.c b/src/stardis-compute.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -24,16 +24,24 @@ #include <star/s3d.h> #include <star/sg3d_sXd_helper.h> +#include <star/ssp.h> #include <rsys/double3.h> #include <rsys/double2.h> #include <rsys/logger.h> #include <rsys/image.h> +#include <rsys/clock_time.h> + +#ifdef COMPILER_GCC +#include <strings.h> /* strcasecmp */ +#else +#include <string.h> +#define strcasecmp(s1, s2) _stricmp((s1), (s2)) +#endif /******************************************************************************* * Local Functions ******************************************************************************/ - struct filter_ctx { const struct stardis* stardis; unsigned prim; @@ -76,19 +84,20 @@ static int hit_filter (const struct s3d_hit* hit, const float ray_org[3], - const float* invalid_, /* In closest_point queries ray_dir is not informed */ + const float ray_dir[3], + const float ray_range[2], void* ray_data, void* filter_data) { struct filter_ctx* filter_ctx = ray_data; float s, dir[3]; - enum senc3d_side hit_side; + enum sg3d_property_type prop; struct s3d_attrib interf_pos; const struct description* descriptions; unsigned descr[SG3D_PROP_TYPES_COUNT__]; const struct stardis* stardis; - (void)ray_org; (void)invalid_; (void)filter_data; + (void)ray_org; (void)ray_dir; (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)); @@ -143,11 +152,11 @@ hit_filter filter_ctx->probe_on_boundary = 0; } - /* Determine which side was hit + /* Determine which side was hit and the property to look at. * Warning: following Embree 2 convention for geometrical normals, * the Star3D hit normals are left-handed */ - hit_side = (s > 0) ? SG3D_BACK : SG3D_FRONT; - if(descr[hit_side] == SG3D_UNSPECIFIED_PROPERTY) { + prop = (s > 0) ? SG3D_BACK : SG3D_FRONT; + if(descr[prop] == SG3D_UNSPECIFIED_PROPERTY) { filter_ctx->prim = hit->prim.prim_id; filter_ctx->desc = NULL; f3_set(filter_ctx->pos, interf_pos.value); @@ -156,7 +165,7 @@ hit_filter filter_ctx->probe_on_boundary = 0; } else { filter_ctx->prim = hit->prim.prim_id; - filter_ctx->desc = descriptions + descr[hit_side]; + filter_ctx->desc = descriptions + descr[prop]; f3_set(filter_ctx->pos, interf_pos.value); filter_ctx->dist = hit->distance; filter_ctx->outside = 0; @@ -173,7 +182,7 @@ check_probe_conform_to_type (const struct stardis* stardis, const int move2boundary, double pos[3], - size_t* iprim, + unsigned* iprim, double uv[2]) { res_T res = RES_OK; @@ -365,19 +374,57 @@ error: goto end; } +#define READ_RANDOM_STATE(Name) \ + if(!str_is_empty(Name)) { \ + const char* name = str_cget(Name); \ + stream = fopen(name, "r"); \ + if(!stream) { \ + res = RES_IO_ERR; \ + logger_print(stardis->logger, LOG_ERROR, \ + "Could not open generator's state file ('%s').\n", \ + name); \ + goto error; \ + } \ + ERR(ssp_rng_create(stardis->allocator, &ssp_rng_mt19937_64, &args.rng_state)); \ + res = read_random_generator_state(args.rng_state, stream); \ + if(res != RES_OK) { \ + logger_print(stardis->logger, LOG_ERROR, \ + "Could not read random generator's state ('%s').\n", \ + name); \ + goto error; \ + } \ + fclose(stream); stream = NULL; \ + } + +#define WRITE_RANDOM_STATE(Name) \ + if(!str_is_empty(Name)) { \ + const char* name = str_cget(Name); \ + stream = fopen(name, "wb"); \ + if(!stream) { \ + res = RES_IO_ERR; \ + logger_print(stardis->logger, LOG_ERROR, \ + "Could not write random generator's state ('%s').\n", \ + name); \ + goto error; \ + } \ + ERR(write_random_generator_state(estimator, stream)); \ + fclose(stream); stream = NULL; \ + } + static res_T -compute_probe(struct stardis* stardis) +compute_probe(struct stardis* stardis, struct time* start) { res_T res = RES_OK; double uv[2] = { 0,0 }; - size_t iprim = SIZE_MAX; + unsigned iprim = UINT_MAX; struct sdis_green_function* green = NULL; struct sdis_estimator* estimator = NULL; struct dump_path_context dump_ctx; struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; FILE* stream = NULL; + struct time compute_start, compute_end; - ASSERT(stardis && (stardis->mode & MODE_PROBE_COMPUTE)); + ASSERT(stardis && start && (stardis->mode & MODE_PROBE_COMPUTE)); ERR(check_probe_conform_to_type(stardis, 0, stardis->probe, &iprim, uv)); @@ -385,9 +432,15 @@ compute_probe(struct stardis* stardis) d3_set(args.position, stardis->probe); 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)) { + time_current(&compute_start); ERR(sdis_solve_probe_green_function(stardis->sdis_scn, &args, &green)); + time_current(&compute_end); if(stardis->mode & MODE_BIN_GREEN) { + struct time output_end; ASSERT(!str_is_empty(&stardis->bin_green_filename)); stream = fopen(str_cget(&stardis->bin_green_filename), "wb"); if(!stream) { @@ -407,56 +460,223 @@ compute_probe(struct stardis* stardis) ERR(dump_paths_end(green, stardis, stream)); fclose(stream); stream = NULL; } + 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; + time_current(&compute_start); ERR(sdis_solve_probe(stardis->sdis_scn, &args, &estimator)); + time_current(&compute_end); + ERR(print_computation_time(estimator, stardis, + start, &compute_start, &compute_end, NULL)); ERR(print_single_MC_result(estimator, stardis, stdout)); - /* Dump paths recorded according to user settings */ + + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); } + /* Output random state? */ + WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); + end: if(stream) fclose(stream); 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 -compute_probe_on_interface(struct stardis* stardis) +compute_probe_on_interface(struct stardis* stardis, struct time* start) { res_T res = RES_OK; double uv[2] = { 0,0 }; - size_t iprim = SIZE_MAX; + 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 = NULL; - - ASSERT(stardis && (stardis->mode & MODE_PROBE_COMPUTE_ON_INTERFACE)); + 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, &iprim, uv)); - ASSERT(iprim != SIZE_MAX); + 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(d->type == DESC_MAT_SOLID || d->type == DESC_MAT_FLUID); + 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(d->type == DESC_MAT_SOLID || d->type == DESC_MAT_FLUID); + 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(fd->type == DESC_MAT_SOLID || fd->type == DESC_MAT_FLUID); + fmat_id = (fd->type == DESC_MAT_SOLID) + ? fd->d.solid.solid_id : fd->d.fluid.fluid_id; + ASSERT(bd->type == DESC_MAT_SOLID || bd->type == DESC_MAT_FLUID); + bmat_id = (bd->type == DESC_MAT_SOLID) + ? bd->d.solid.solid_id : bd->d.fluid.fluid_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; d3_set(args.uv, uv); - args.side = SDIS_FRONT; + 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)) { + 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; ASSERT(!str_is_empty(&stardis->bin_green_filename)); stream = fopen(str_cget(&stardis->bin_green_filename), "wb"); if(!stream) { @@ -476,24 +696,36 @@ compute_probe_on_interface(struct stardis* stardis) ERR(dump_paths_end(green, stardis, stream)); fclose(stream); stream = NULL; } + 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; + 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)); ERR(print_single_MC_result(estimator, stardis, stdout)); - /* Dump paths recorded according to user settings */ + + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); } + /* Output random state? */ + WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); + end: if(stream) fclose(stream); 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; @@ -578,7 +810,7 @@ error: } static res_T -compute_camera(struct stardis* stardis) +compute_camera(struct stardis* stardis, struct time* start) { res_T res = RES_OK; double proj_ratio; @@ -589,8 +821,10 @@ compute_camera(struct stardis* stardis) struct dump_path_context dump_ctx; size_t definition[2]; size_t ix, iy; + FILE* stream; + struct time compute_start, compute_end, output_end; - ASSERT(stardis + ASSERT(stardis && start && !(stardis->mode & (MODE_BIN_GREEN | MODE_GREEN)) && (stardis->mode & MODE_IR_COMPUTE)); @@ -621,16 +855,29 @@ compute_camera(struct stardis* stardis) args.register_paths = stardis->dump_paths; /* Launch the simulation */ + time_current(&compute_start); ERR(sdis_solve_camera(stardis->sdis_scn, &args, &buf)); + time_current(&compute_end); /* Write the image */ + if(str_is_empty(&stardis->camera.file_name)) + stream = stdout; + else { + stream = fopen(str_cget(&stardis->camera.file_name), "w"); + if(!stream) { + res = RES_IO_ERR; + goto error; + } + } ASSERT(stardis->camera.fmt == STARDIS_RENDERING_OUTPUT_FILE_FMT_VTK || stardis->camera.fmt == STARDIS_RENDERING_OUTPUT_FILE_FMT_HT); if(stardis->camera.fmt == STARDIS_RENDERING_OUTPUT_FILE_FMT_VTK) - ERR(dump_vtk_image(buf, stdout)); - else ERR(dump_ht_image(buf, stdout)); + ERR(dump_vtk_image(buf, stream)); + else ERR(dump_ht_image(buf, stream)); + if(!str_is_empty(&stardis->camera.file_name)) + fclose(stream); - /* Dump paths recorded according to user settings */ + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_buffer_get_definition(buf, definition)); @@ -641,6 +888,9 @@ compute_camera(struct stardis* stardis) ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); } } + time_current(&output_end); + ERR(print_computation_time(NULL, stardis, + start, &compute_start, &compute_end, NULL)); end: if(cam) SDIS(camera_ref_put(cam)); @@ -651,7 +901,7 @@ error: } static res_T -compute_medium(struct stardis* stardis) +compute_medium(struct stardis* stardis, struct time* start) { res_T res = RES_OK; struct sdis_medium* medium = NULL; @@ -660,8 +910,9 @@ compute_medium(struct stardis* stardis) struct sdis_solve_medium_args args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; struct dump_path_context dump_ctx; FILE* stream = NULL; + struct time compute_start, compute_end; - ASSERT(stardis && (stardis->mode & MODE_MEDIUM_COMPUTE)); + ASSERT(stardis && start && (stardis->mode & MODE_MEDIUM_COMPUTE)); /* Find medium */ medium = find_medium_by_name(stardis, &stardis->solve_name, NULL); @@ -677,9 +928,15 @@ compute_medium(struct stardis* stardis) args.medium = medium; 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)) { + time_current(&compute_start); ERR(sdis_solve_medium_green_function(stardis->sdis_scn, &args, &green)); + time_current(&compute_end); if(stardis->mode & MODE_BIN_GREEN) { + struct time output_end; ASSERT(!str_is_empty(&stardis->bin_green_filename)); stream = fopen(str_cget(&stardis->bin_green_filename), "wb"); if(!stream) { @@ -699,24 +956,35 @@ compute_medium(struct stardis* stardis) ERR(dump_paths_end(green, stardis, stream)); fclose(stream); stream = NULL; } + 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; + time_current(&compute_start); ERR(sdis_solve_medium(stardis->sdis_scn, &args, &estimator)); + time_current(&compute_end); + ERR(print_computation_time(estimator, stardis, + start, &compute_start, &compute_end, NULL)); ERR(print_single_MC_result(estimator, stardis, stdout)); - /* Dump paths recorded according to user settings */ + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); } + /* Output random state? */ + WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); + end: if(stream) fclose(stream); 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; @@ -771,7 +1039,7 @@ error: } static res_T -compute_boundary(struct stardis* stardis) +compute_boundary(struct stardis* stardis, struct time* start) { res_T res = RES_OK; struct sdis_green_function* green = NULL; @@ -779,8 +1047,9 @@ compute_boundary(struct stardis* stardis) struct dump_path_context dump_ctx; struct sdis_solve_boundary_args args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; FILE* stream = NULL; + struct time compute_start, compute_end; - ASSERT(stardis && (stardis->mode & MODE_BOUNDARY_COMPUTE)); + ASSERT(stardis && start && (stardis->mode & MODE_BOUNDARY_COMPUTE)); args.nrealisations = stardis->samples; args.primitives @@ -790,9 +1059,15 @@ compute_boundary(struct stardis* stardis) = darray_size_t_size_get(&stardis->compute_surface.primitives); 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)) { + time_current(&compute_start); ERR(sdis_solve_boundary_green_function(stardis->sdis_scn, &args, &green)); + time_current(&compute_end); if(stardis->mode & MODE_BIN_GREEN) { + struct time output_end; ASSERT(!str_is_empty(&stardis->bin_green_filename)); stream = fopen(str_cget(&stardis->bin_green_filename), "wb"); if(!stream) { @@ -812,31 +1087,42 @@ compute_boundary(struct stardis* stardis) ERR(dump_paths_end(green, stardis, stream)); fclose(stream); stream = NULL; } + 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; + time_current(&compute_start); ERR(sdis_solve_boundary(stardis->sdis_scn, &args, &estimator)); + time_current(&compute_end); + ERR(print_computation_time(estimator, stardis, + start, &compute_start, &compute_end, NULL)); ERR(print_single_MC_result(estimator, stardis, stdout)); - /* Dump paths recorded according to user settings */ + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); } + /* Output random state? */ + WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); + end: if(stream) fclose(stream); 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 -compute_flux_boundary(struct stardis* stardis) +compute_flux_boundary(struct stardis* stardis, struct time* start) { res_T res = RES_OK; struct sdis_green_function* green = NULL; @@ -844,8 +1130,10 @@ compute_flux_boundary(struct stardis* stardis) struct sdis_solve_boundary_flux_args args = SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT; struct dump_path_context dump_ctx; + FILE* stream = NULL; + struct time compute_start, compute_end; - ASSERT(stardis && (stardis->mode & MODE_FLUX_BOUNDARY_COMPUTE)); + ASSERT(stardis && start && (stardis->mode & MODE_FLUX_BOUNDARY_COMPUTE)); args.nrealisations = stardis->samples; args.primitives @@ -854,24 +1142,35 @@ compute_flux_boundary(struct stardis* stardis) = darray_size_t_size_get(&stardis->compute_surface.primitives); d2_set(args.time_range, stardis->time_range); + /* Input random state? */ + READ_RANDOM_STATE(&stardis->rndgen_state_in_filename); + + time_current(&compute_start); ERR(sdis_solve_boundary_flux(stardis->sdis_scn, &args, &estimator)); + time_current(&compute_end); + ERR(print_computation_time(estimator, stardis, + start, &compute_start, &compute_end, NULL)); ERR(print_single_MC_result(estimator, stardis, stdout)); - /* Dump paths recorded according to user settings */ + /* Dump recorded paths according to user settings */ dump_ctx.stardis = stardis; dump_ctx.rank = 0; ERR(sdis_estimator_for_each_path(estimator, dump_path, &dump_ctx)); + /* 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 -compute_map(struct stardis* stardis) +compute_map(struct stardis* stardis, struct time* start) { res_T res = RES_OK; struct sdis_green_function* green = NULL; @@ -880,7 +1179,8 @@ compute_map(struct stardis* stardis) int estimators_initialized = 0; size_t p; - ASSERT(stardis + (void)start; + ASSERT(stardis && start && (stardis->mode & MODE_MAP_COMPUTE) && !(stardis->mode & (MODE_BIN_GREEN | MODE_GREEN))); @@ -921,6 +1221,9 @@ error: goto end; } +#undef READ_RANDOM_STATE +#undef WRITE_RANDOM_STATE + /******************************************************************************* * Public Functions ******************************************************************************/ @@ -972,6 +1275,8 @@ read_compute_surface struct add_geom_ctx add_geom_ctx; struct sg3d_geometry_add_callbacks callbacks = SG3D_ADD_CALLBACKS_NULL__; const char* file; + size_t i; + double _2area = 0; ASSERT(stardis); @@ -1006,6 +1311,26 @@ read_compute_surface goto error; } + /* Compute area of the compute surface */ + FOR_EACH(i, 0, add_geom_ctx.stl_desc.triangles_count) { + const unsigned* trg = add_geom_ctx.stl_desc.indices + (i * 3); + const float* v0f = add_geom_ctx.stl_desc.vertices + (trg[0] * 3); + const float* v1f = add_geom_ctx.stl_desc.vertices + (trg[1] * 3); + const float* v2f = add_geom_ctx.stl_desc.vertices + (trg[2] * 3); + double v0[3], v1[3], v2[3], edge0[3], edge1[3], normal[3], norm; + + /* Compute component area and volume */ + d3_sub(edge0, d3_set_f3(v1, v1f), d3_set_f3(v0, v0f)); + d3_sub(edge1, d3_set_f3(v2, v2f), d3_set_f3(v0, v0f)); + d3_cross(normal, edge0, edge1); + norm = d3_normalize(normal, normal); + ASSERT(norm); + _2area += norm; + } + + stardis->compute_surface.area = _2area * 0.5 + * stardis->scale_factor * stardis->scale_factor; + end: if(sstl) SSTL(ref_put(sstl)); return res; @@ -1015,27 +1340,28 @@ error: res_T stardis_compute - (struct stardis* stardis) + (struct stardis* stardis, + struct time* start) { res_T res = RES_OK; - ASSERT(stardis && (stardis->mode & COMPUTE_MODES)); + ASSERT(stardis && start && (stardis->mode & COMPUTE_MODES)); /* Compute */ if(stardis->mode & MODE_PROBE_COMPUTE) - ERR(compute_probe(stardis)); + ERR(compute_probe(stardis, start)); else if(stardis->mode & MODE_PROBE_COMPUTE_ON_INTERFACE) - ERR(compute_probe_on_interface(stardis)); + ERR(compute_probe_on_interface(stardis, start)); else if(stardis->mode & MODE_IR_COMPUTE) - ERR(compute_camera(stardis)); + ERR(compute_camera(stardis, start)); else if(stardis->mode & MODE_MEDIUM_COMPUTE) - ERR(compute_medium(stardis)); + ERR(compute_medium(stardis, start)); else if(stardis->mode & MODE_BOUNDARY_COMPUTE) - ERR(compute_boundary(stardis)); + ERR(compute_boundary(stardis, start)); else if(stardis->mode & MODE_FLUX_BOUNDARY_COMPUTE) - ERR(compute_flux_boundary(stardis)); + ERR(compute_flux_boundary(stardis, start)); else if(stardis->mode & MODE_MAP_COMPUTE) - ERR(compute_map(stardis)); + ERR(compute_map(stardis, start)); else FATAL("Unknown mode.\n"); end: diff --git a/src/stardis-compute.h b/src/stardis-compute.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -20,6 +20,7 @@ #include <rsys/dynamic_array.h> struct stardis; +struct time; struct str; #define DARRAY_NAME estimators @@ -34,7 +35,8 @@ find_medium_by_name extern LOCAL_SYM res_T stardis_compute - (struct stardis* stardis); + (struct stardis* stardis, + struct time* start); extern LOCAL_SYM res_T read_compute_surface diff --git a/src/stardis-fluid.c b/src/stardis-fluid.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 diff --git a/src/stardis-fluid.h b/src/stardis-fluid.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 diff --git a/src/stardis-intface.c b/src/stardis-intface.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -75,6 +75,17 @@ interface_get_alpha return interface_props->alpha; } +static double +interface_get_tcr + (const struct sdis_interface_fragment* frag, + struct sdis_data* data) +{ + const struct intface* interface_props = sdis_data_cget(data); + (void)frag; + return interface_props->tcr; +} + + /******************************************************************************* * Public Functions ******************************************************************************/ @@ -94,7 +105,8 @@ create_intface struct sdis_data* data = NULL; unsigned fd, bd, cd, descr[SG3D_PROP_TYPES_COUNT__]; int fluid_count = 0, solid_count = 0; - int solid_fluid_connection_count = 0, connection_count = 0, boundary_count = 0; + int solid_fluid_connection_count = 0, solid_solid_connection_count = 0; + int connection_count = 0, boundary_count = 0; const struct description* descriptions; struct sdis_medium** media; res_T res = RES_OK; @@ -185,7 +197,8 @@ create_intface struct sdis_medium* def_medium = NULL; unsigned ext_id; if(connect->type != DESC_SOLID_FLUID_CONNECT - && front_defined == back_defined) + && connect->type != DESC_SOLID_SOLID_CONNECT + && front_defined == back_defined) { /* 1 and only 1 defined */ res = RES_BAD_ARG; @@ -335,6 +348,20 @@ create_intface fluid_side_shader->specular_fraction = interface_get_alpha; } break; + case DESC_SOLID_SOLID_CONNECT: + /* Both front and back should be defined */ + if(solid_count != 2 || fluid_count != 0) { + res = RES_BAD_ARG; + goto error; + } + ASSERT(front_defined && back_defined); + connection_count++; + solid_solid_connection_count++; + interface_props->tcr = connect->d.ss_connect.tcr; + if(connect->d.ss_connect.tcr > 0) { + interface_shader.thermal_contact_resistance = interface_get_tcr; + } + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } diff --git a/src/stardis-intface.h b/src/stardis-intface.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -27,9 +27,12 @@ struct dummies; * Interface data ******************************************************************************/ struct intface { + /* fluid - solid */ double hc; double emissivity; double alpha; + /* solid - solid */ + double tcr; /* Imposed compute temperature & flux */ double imposed_temperature; double imposed_flux; diff --git a/src/stardis-main.c b/src/stardis-main.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -21,6 +21,7 @@ #include <rsys/rsys.h> #include <rsys/mem_allocator.h> #include <rsys/logger.h> +#include <rsys/clock_time.h> #include <stdlib.h> #include <stdio.h> @@ -38,8 +39,11 @@ main struct mem_allocator allocator; struct logger logger; int mode = 0; + struct time start; res_T res = RES_OK; + time_current(&start); + ERR(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); allocator_initialized = 1; @@ -52,7 +56,7 @@ main ERR(init_args(&logger, &allocator, &args)); args_initialized = 1; - ERR(parse_args(argc, argv, args)); + ERR(parse_args(argc, argv, args, &allocator)); mode = args->mode; if(mode & MODE_DUMP_HELP) { @@ -99,7 +103,7 @@ main } ASSERT(mode & COMPUTE_MODES); - ERR(stardis_compute(&stardis)); + ERR(stardis_compute(&stardis, &start)); exit: if(args_initialized) release_args(args); diff --git a/src/stardis-output.c b/src/stardis-output.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -25,6 +25,7 @@ #include<star/senc3d.h> #include<star/sg3d.h> +#include<star/ssp.h> #include <rsys/math.h> #include <rsys/mem_allocator.h> @@ -32,6 +33,7 @@ #include <rsys/hash_table.h> #include <rsys/logger.h> #include <rsys/str.h> +#include <rsys/clock_time.h> #include <limits.h> #include <string.h> @@ -428,61 +430,66 @@ dump_sample_end void* ctx) { res_T res = RES_OK; - struct sdis_point pt = SDIS_POINT_NULL; - struct sdis_data* data = NULL; - enum sdis_medium_type type; struct e_ctx* e_ctx = ctx; - const struct description* descs; + enum sdis_green_path_end_type end_type; FILE* stream; double elapsed; - double* pos; - unsigned id; CHK(path && ctx); stream = e_ctx->stream; - ERR(sdis_green_path_get_limit_point(path, &pt)); ERR(sdis_green_path_get_elapsed_time(path, &elapsed)); - - descs = darray_descriptions_cdata_get(e_ctx->desc); - switch (pt.type) { - case SDIS_FRAGMENT: { - struct intface* d__; - data = sdis_interface_get_data(pt.data.itfrag.intface); - pos = pt.data.itfrag.fragment.P; - d__ = sdis_data_get(data); - id = d__->desc_id; - CHK(DESC_IS_T(descs[id].type) || DESC_IS_H(descs[id].type)); - break; - } - case SDIS_VERTEX: - type = sdis_medium_get_type(pt.data.mdmvert.medium); - data = sdis_medium_get_data(pt.data.mdmvert.medium); - pos = pt.data.mdmvert.vertex.P; - if(pt.data.mdmvert.vertex.P[0] == INF) { - /* Radiative output (ambient temperature) */ - size_t sz = darray_descriptions_size_get(e_ctx->desc); - ASSERT(sz <= UINT_MAX); - id = (unsigned)sz; /* Ambient ID */ - } - else if(type == SDIS_FLUID) { - struct fluid* d__ = sdis_data_get(data); - id = d__->desc_id; - } else { - struct solid* d__ = sdis_data_get(data); - ASSERT(type == SDIS_SOLID); - ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */ - id = d__->desc_id; - } - break; - default: FATAL("Unreachable code.\n"); break; - } - - /* End, End ID, X, Y, Z, Elapsed time */ - if(pt.data.mdmvert.vertex.P[0] == INF) { + ERR(sdis_green_path_get_end_type(path, &end_type)); + if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) { + size_t ambient_id = darray_descriptions_size_get(e_ctx->desc); + ASSERT(ambient_id <= UINT_MAX); + /* End, End ID, X, Y, Z, Elapsed time */ fprintf(stream, "AMBIANT, %u, 0, 0, 0, %g\n", - id, elapsed); + (unsigned)ambient_id, elapsed); } else { + struct sdis_point pt = SDIS_POINT_NULL; + struct sdis_data* data = NULL; + enum sdis_medium_type type; + const struct description* descs; + double* pos; + unsigned id; + + ERR(sdis_green_path_get_limit_point(path, &pt)); + + descs = darray_descriptions_cdata_get(e_ctx->desc); + switch(pt.type) { + case SDIS_FRAGMENT: { + struct intface* d__; + data = sdis_interface_get_data(pt.data.itfrag.intface); + pos = pt.data.itfrag.fragment.P; + d__ = sdis_data_get(data); + id = d__->desc_id; + CHK(DESC_IS_T(descs[id].type) || DESC_IS_H(descs[id].type)); + break; + } + case SDIS_VERTEX: + type = sdis_medium_get_type(pt.data.mdmvert.medium); + data = sdis_medium_get_data(pt.data.mdmvert.medium); + pos = pt.data.mdmvert.vertex.P; + if(pt.data.mdmvert.vertex.P[0] == INF) { + /* Radiative output (ambient temperature) */ + size_t sz = darray_descriptions_size_get(e_ctx->desc); + ASSERT(sz <= UINT_MAX); + id = (unsigned)sz; /* Ambient ID */ + } + else if(type == SDIS_FLUID) { + struct fluid* d__ = sdis_data_get(data); + id = d__->desc_id; + } else { + struct solid* d__ = sdis_data_get(data); + ASSERT(type == SDIS_SOLID); + ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */ + id = d__->desc_id; + } + break; + default: FATAL("Unreachable code.\n"); break; + } + /* End, End ID, X, Y, Z, Elapsed time */ fprintf(stream, "%s, %u, %g, %g, %g, %g\n", str_cget(get_description_name(descs + id)), id, SPLIT3(pos), elapsed); } @@ -499,65 +506,75 @@ dump_sample void* ctx) { res_T res = RES_OK; - struct sdis_point pt = SDIS_POINT_NULL; - struct sdis_data* data = NULL; - enum sdis_medium_type type; struct htable_weigth_iterator it, end; struct path_header header; struct w_ctx* w_ctx = ctx; - const struct description* descs; + enum sdis_green_path_end_type end_type; FILE* stream; unsigned* ids = NULL; double* weights = NULL; size_t sz, i; - double t0; CHK(path && ctx); stream = w_ctx->stream; - ERR(sdis_green_path_get_limit_point(path, &pt)); - - /* For each path, dump: - * # end_id #power_terms #flux_terms - * power_id_1 ... power_id_n flux_id_1 ... flux_id_n - * power_factor_1 ... power_factor_n flux_factor_1 ... flux_factor_n - */ - - descs = darray_descriptions_cdata_get(w_ctx->desc); - switch (pt.type) { - case SDIS_FRAGMENT: { - struct intface* d__; - unsigned desc_id; - data = sdis_interface_get_data(pt.data.itfrag.intface); - d__ = sdis_data_get(data); - desc_id = d__->desc_id; - CHK(DESC_IS_T(descs[desc_id].type) || DESC_IS_H(descs[desc_id].type)); - header.id = desc_id; + ERR(sdis_green_path_get_end_type(path, &end_type)); + if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) { + size_t ambient_id = darray_descriptions_size_get(w_ctx->desc); + ASSERT(ambient_id <= UINT_MAX); header.at_initial = 0; - break; - } - case SDIS_VERTEX: - type = sdis_medium_get_type(pt.data.mdmvert.medium); - data = sdis_medium_get_data(pt.data.mdmvert.medium); - t0 = medium_get_t0(pt.data.mdmvert.medium); - header.at_initial = (pt.data.mdmvert.vertex.time <= t0); - if(pt.data.mdmvert.vertex.P[0] == INF) { - /* Radiative output (ambient temperature) */ - sz = darray_descriptions_size_get(w_ctx->desc); - ASSERT(sz <= UINT_MAX); - header.id = (unsigned)sz; /* Ambient ID */ - } - else if(type == SDIS_FLUID) { - struct fluid* d__ = sdis_data_get(data); - header.id = d__->desc_id; - } else { - struct solid* d__ = sdis_data_get(data); - ASSERT(type == SDIS_SOLID); - ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */ - header.id = d__->desc_id; + header.id = (unsigned)ambient_id; + } else { + struct sdis_point pt = SDIS_POINT_NULL; + struct sdis_data* data = NULL; + enum sdis_medium_type type; + const struct description* descs; + double t0; + + ERR(sdis_green_path_get_limit_point(path, &pt)); + + /* For each path, dump: + * # end_id #power_terms #flux_terms + * power_id_1 ... power_id_n flux_id_1 ... flux_id_n + * power_factor_1 ... power_factor_n flux_factor_1 ... flux_factor_n + */ + + descs = darray_descriptions_cdata_get(w_ctx->desc); + switch(pt.type) { + case SDIS_FRAGMENT: { + struct intface* d__; + unsigned desc_id; + data = sdis_interface_get_data(pt.data.itfrag.intface); + d__ = sdis_data_get(data); + desc_id = d__->desc_id; + CHK(DESC_IS_T(descs[desc_id].type) || DESC_IS_H(descs[desc_id].type)); + header.id = desc_id; + header.at_initial = 0; + break; + } + case SDIS_VERTEX: + type = sdis_medium_get_type(pt.data.mdmvert.medium); + data = sdis_medium_get_data(pt.data.mdmvert.medium); + t0 = medium_get_t0(pt.data.mdmvert.medium); + header.at_initial = (pt.data.mdmvert.vertex.time <= t0); + if(pt.data.mdmvert.vertex.P[0] == INF) { + /* Radiative output (ambient temperature) */ + sz = darray_descriptions_size_get(w_ctx->desc); + ASSERT(sz <= UINT_MAX); + header.id = (unsigned)sz; /* Ambient ID */ + } + else if(type == SDIS_FLUID) { + struct fluid* d__ = sdis_data_get(data); + header.id = d__->desc_id; + } else { + struct solid* d__ = sdis_data_get(data); + ASSERT(type == SDIS_SOLID); + ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */ + header.id = d__->desc_id; + } + break; + default: FATAL("Unreachable code.\n"); break; } - break; - default: FATAL("Unreachable code.\n"); break; } /* Merge power and flux terms */ @@ -638,7 +655,13 @@ dump_green_bin char* name_pool = NULL; char* pool_ptr; const char green_string[] = "GREEN_BIN_FILE:"; - const unsigned file_fmt_version = 1; + const unsigned file_fmt_version = 2; + /* The following type must be identical to its stardis-green counterpart! */ + struct bfile_green_counts { + unsigned desc_count, smed_count, fmed_count, tbound_count, hbound_count, + fbound_count, sfconnect_count, ssconnect_count, name_pool_sz; + size_t ok_count, failed_count; + } file_counts; ASSERT(green && stardis && stream); @@ -647,6 +670,11 @@ dump_green_bin sz = darray_descriptions_size_get(&stardis->descriptions); ASSERT(sz <= UINT_MAX); szd = (unsigned)sz; + ASSERT(szd == + (stardis->counts.smed_count + stardis->counts.fmed_count + + stardis->counts.tbound_count + stardis->counts.hbound_count + + stardis->counts.fbound_count + stardis->counts.sfconnect_count + + stardis->counts.ssconnect_count)); descs = darray_descriptions_cdata_get(&stardis->descriptions); /* Save names that do not fit inplace */ @@ -675,11 +703,18 @@ dump_green_bin FW(&file_fmt_version, 1); /* Write counts */ - FW(&szd, 1); - FW(&stardis->counts, 1); - FW(&name_pool_sz, 1); - FW(&ok_count, 1); - FW(&failed_count, 1); + file_counts.desc_count = szd; + file_counts.smed_count = stardis->counts.smed_count; + file_counts.fmed_count = stardis->counts.fmed_count; + file_counts.tbound_count = stardis->counts.tbound_count; + file_counts.hbound_count = stardis->counts.hbound_count; + file_counts.fbound_count = stardis->counts.fbound_count; + file_counts.sfconnect_count = stardis->counts.sfconnect_count; + file_counts.ssconnect_count = stardis->counts.ssconnect_count; + file_counts.name_pool_sz = name_pool_sz; + file_counts.ok_count = ok_count; + file_counts.failed_count = failed_count; + FW(&file_counts, 1); /* Write descriptions*/ FW(descs, szd); @@ -1300,6 +1335,50 @@ error: } res_T +print_computation_time + (struct sdis_estimator* estimator, + struct stardis* stardis, + struct time* start, + struct time* compute_start, + struct time* compute_end, + struct time* output_end) +{ + res_T res = RES_OK; + struct sdis_mc time; + struct time tmp; + char buf[128]; + const int flag = TIME_MSEC | TIME_SEC | TIME_MIN | TIME_HOUR | TIME_DAY; + + ASSERT(stardis && start && compute_start && compute_end); + + time_sub(&tmp, compute_start, start); + time_dump(&tmp, flag, NULL, buf, sizeof(buf)); + logger_print(stardis->logger, LOG_OUTPUT, + "Initialisation time = %s\n", buf); + time_sub(&tmp, compute_end, compute_start); + time_dump(&tmp, flag, NULL, buf, sizeof(buf)); + logger_print(stardis->logger, LOG_OUTPUT, + "Computation time = %s\n", buf); + if(output_end) { + time_sub(&tmp, output_end, compute_end); + time_dump(&tmp, flag, NULL, buf, sizeof(buf)); + logger_print(stardis->logger, LOG_OUTPUT, + "Result output time = %s\n", buf); + } + + if(estimator) { + ERR(sdis_estimator_get_realisation_time(estimator, &time)); + logger_print(stardis->logger, LOG_OUTPUT, + "Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); + } + +exit: + return res; +error: + goto exit; +} + +res_T print_single_MC_result (struct sdis_estimator* estimator, struct stardis* stardis, @@ -1330,12 +1409,13 @@ print_single_MC_result case MODE_PROBE_COMPUTE: 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 +/- %g\n", + 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 +/- %g\n", + 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 */ @@ -1346,12 +1426,14 @@ print_single_MC_result case MODE_PROBE_COMPUTE_ON_INTERFACE: if(stardis->mode & MODE_EXTENDED_RESULTS) { if(stardis->time_range[0] == stardis->time_range[1]) - fprintf(stream, "Boundary temperature at [%g, %g, %g] at t=%g = %g +/- %g\n", + fprintf(stream, + "Boundary 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, "Boundary temperature at [%g, %g, %g] with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Boundary 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 */ @@ -1362,12 +1444,13 @@ print_single_MC_result case MODE_MEDIUM_COMPUTE: if(stardis->mode & MODE_EXTENDED_RESULTS) { if(stardis->time_range[0] == stardis->time_range[1]) - fprintf(stream, "Temperature in medium '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Temperature in medium '%s' at t=%g = %g K +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], result.E, /* Expected value */ result.SE); /* Standard error */ else - fprintf(stream, "Temperature in medium '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Temperature in medium '%s' with t in [%g %g] = %g K +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), result.E, /* Expected value */ result.SE); /* Standard error */ @@ -1378,12 +1461,13 @@ print_single_MC_result case MODE_BOUNDARY_COMPUTE: if(stardis->mode & MODE_EXTENDED_RESULTS) { if(stardis->time_range[0] == stardis->time_range[1]) - fprintf(stream, "Temperature at boundary '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], result.E, /* Expected value */ result.SE); /* Standard error */ else - fprintf(stream, "Temperature at boundary '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), result.E, /* Expected value */ result.SE); /* Standard error */ @@ -1398,69 +1482,82 @@ print_single_MC_result if(stardis->mode & MODE_EXTENDED_RESULTS) { if(stardis->time_range[0] == stardis->time_range[1]) - fprintf(stream, "Temperature at boundary '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], result.E, /* Expected value */ result.SE); /* Standard error */ else - fprintf(stream, "Temperature at boundary '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n", str_cget(&stardis->solve_name), 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 at boundary '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Convective flux at boundary '%s' at t=%g = %g W +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * result.SE); /* Standard error */ else - fprintf(stream, "Convective flux at boundary '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Convective flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * 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 at boundary '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Radiative flux at boundary '%s' at t=%g = %g W +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * result.SE); /* Standard error */ else - fprintf(stream, "Radiative flux at boundary '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Radiative flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * 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 at boundary '%s' at t=%g = %g +/- %g\n", + fprintf(stream, "Imposed flux at boundary '%s' at t=%g = %g W +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * result.SE); /* Standard error */ else - fprintf(stream, "Imposed flux at boundary '%s' with t in [%g %g] = %g +/- %g\n", + fprintf(stream, + "Imposed flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + 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 +/- %g\n", + fprintf(stream, "Total flux Flux at boundary '%s' at t=%g = %g W +/- %g\n", str_cget(&stardis->solve_name), stardis->time_range[0], - result.E, /* Expected value */ - result.SE); /* Standard error */ + 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 +/- %g\n", + fprintf(stream, + "Total flux Flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), - result.E, /* Expected value */ - result.SE); /* Standard error */ + stardis->compute_surface.area * result.E, /* Expected value */ + stardis->compute_surface.area * result.SE); /* Standard error */ } else { - fprintf(stream, "%g %g ", result.E, result.SE); + 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); + fprintf(stream, "%g %g ", + stardis->compute_surface.area * result.E, + stardis->compute_surface.area * result.SE); ERR(sdis_estimator_get_radiative_flux(estimator, &result)); - fprintf(stream, "%g %g ", result.E, result.SE); + fprintf(stream, "%g %g ", + stardis->compute_surface.area * result.E, + stardis->compute_surface.area * result.SE); ERR(sdis_estimator_get_imposed_flux(estimator, &result)); - fprintf(stream, "%g %g ", result.E, result.SE); + fprintf(stream, "%g %g ", + stardis->compute_surface.area * result.E, + stardis->compute_surface.area * result.SE); ERR(sdis_estimator_get_total_flux(estimator, &result)); - fprintf(stream, "%g %g ", result.E, result.SE); + fprintf(stream, "%g %g ", + stardis->compute_surface.area * result.E, + stardis->compute_surface.area * result.SE); fprintf(stream, "%lu %lu\n", nfailures, nsamples); } break; @@ -1727,3 +1824,23 @@ exit: error: goto exit; } + +res_T +write_random_generator_state + (struct sdis_estimator* estimator, + FILE* stream) +{ + res_T res; + struct ssp_rng* state; + res = sdis_estimator_get_rng_state(estimator, &state); + if(res != RES_OK) return res; + return ssp_rng_write(state, stream); +} + +res_T +read_random_generator_state + (struct ssp_rng* state, + FILE* stream) +{ + return ssp_rng_read(state, stream); +} diff --git a/src/stardis-output.h b/src/stardis-output.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -31,6 +31,8 @@ struct stardis; struct geometry; struct vertex; struct darray_estimators; +struct time; +struct ssp_rng; struct dump_path_context { unsigned long rank; @@ -81,6 +83,15 @@ dump_enclosure_related_stuff_at_the_end_of_vtk FILE* stream); extern res_T +print_computation_time + (struct sdis_estimator* estimator, /* Can be NULL */ + struct stardis* stardis, + struct time* start, + struct time* computation_start, + struct time* computation_end, + struct time* output_end); /* Can be NULL */ + +extern res_T print_single_MC_result (struct sdis_estimator* estimator, struct stardis* stardis, @@ -107,4 +118,14 @@ dump_model_as_c_chunks (struct stardis* stardis, FILE* stream); +extern res_T +write_random_generator_state + (struct sdis_estimator* estimator, + FILE* stream); + +extern res_T +read_random_generator_state + (struct ssp_rng* state, + FILE* stream); + #endif diff --git a/src/stardis-parsing.c b/src/stardis-parsing.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -30,26 +30,16 @@ #include <stdio.h> #include <ctype.h> #include <string.h> +#ifdef COMPILER_GCC +#include <strings.h> /* strcasecmp */ +#else +#define strcasecmp(s1, s2) _stricmp((s1), (s2)) +#endif /******************************************************************************* * Local Functions ******************************************************************************/ -#ifdef COMPILER_GCC -static char* -_strupr - (char* s) -{ - char* ptr; - for(ptr = s; *ptr; ++ptr) { - int tmp = toupper(*ptr); - ASSERT(tmp == (char)tmp); - *ptr = (char)tmp; - } - return s; -} -#endif - static char** split_line (char* a_str, @@ -225,17 +215,16 @@ read_sides_and_files } else break; } - _strupr(tk); add_geom_ctx.properties[SG3D_INTFACE] = SG3D_UNSPECIFIED_PROPERTY; - if(0 == strcmp(tk, "FRONT")) { + if(0 == strcasecmp(tk, "FRONT")) { add_geom_ctx.properties[SG3D_FRONT] = description_id; add_geom_ctx.properties[SG3D_BACK] = SG3D_UNSPECIFIED_PROPERTY; } - else if(0 == strcmp(tk, "BACK")) { + else if(0 == strcasecmp(tk, "BACK")) { add_geom_ctx.properties[SG3D_FRONT] = SG3D_UNSPECIFIED_PROPERTY; add_geom_ctx.properties[SG3D_BACK] = description_id; } - else if(0 == strcmp(tk, "BOTH")) { + else if(0 == strcasecmp(tk, "BOTH")) { add_geom_ctx.properties[SG3D_FRONT] = description_id; add_geom_ctx.properties[SG3D_BACK] = description_id; } @@ -346,7 +335,8 @@ init_args /* Set default values */ args->samples = STARDIS_DEFAULT_SAMPLES_COUNT; args->nthreads = SDIS_NTHREADS_DEFAULT; - d2(args->pos_and_time+3, STARDIS_DEFAULT_COMPUTE_TIME); + d2(args->pos_and_time+3, + STARDIS_DEFAULT_COMPUTE_TIME, STARDIS_DEFAULT_COMPUTE_TIME); args->ambient_temp = STARDIS_DEFAULT_AMBIENT_TEMP; args->ref_temp = STARDIS_DEFAULT_REFERENCE_TEMP; args->verbose = STARDIS_DEFAULT_VERBOSE_LEVEL; @@ -512,15 +502,31 @@ short_help fprintf(stream, "\n -V LEVEL\n"); fprintf(stream, " Set the verbosity level.\n"); + + fprintf(stream, "\n -x <FILE>\n"); + fprintf(stream, " Use a random generator's state read from a file.\n"); + + fprintf(stream, "\n -X <FILE>\n"); + fprintf(stream, " Save the final random generator's state in a file.\n"); +} + +#define FREE_AARRAY(ARRAY) \ +if(ARRAY) {\ + int i__ = 0; \ + for(i__=0; *((ARRAY)+i__);i__++){\ + free((ARRAY)[i__]);\ + }\ + free(ARRAY);\ + (ARRAY) = NULL;\ } /* Workaround for a gcc warning when GET_OPTIONAL_TIME_RANGE used with Rank=0 */ -FINLINE int is_less(size_t a, size_t b) { return a < b; } +static FINLINE int is_less(size_t a, size_t b) { return a < b; } /* Get a time range from a coma-separated list of doubles * The first Rank values are mandatory, followed by an optional time range * that can be a single time */ -#define GET_OPTIONAL_TIME_RANGE(Src, Rank, Dst, Logger, OptionString, Option) \ +#define GET_OPTIONAL_TIME_RANGE(Src, Rank, Dst, Logger, OptionString, Option, FullSrc) \ res = cstr_to_list_double((Src), ',', (Dst), &len, (Rank)+2); \ if(res != RES_OK \ || is_less(len, (Rank)) \ @@ -531,7 +537,7 @@ FINLINE int is_less(size_t a, size_t b) { return a < b; } if(res == RES_OK) res = RES_BAD_ARG; \ logger_print((Logger), LOG_ERROR, \ "Invalid argument for option "OptionString": %s\n", \ - (Option), (Src)); \ + (Option), (FullSrc)); \ goto error; \ } else { \ if(len == (Rank)+1) (Dst)[(Rank)+1] = (Dst)[(Rank)];\ @@ -541,29 +547,33 @@ FINLINE int is_less(size_t a, size_t b) { return a < b; } #define GET_STR_AND_OPTIONAL_TIME_RANGE(Str, Time) \ ptr = strchr(optarg, ','); /* First ',' */ \ if(ptr) { /* Time range provided */ \ - GET_OPTIONAL_TIME_RANGE(ptr+1, 0, (Time), args->logger, "-%c", opt); \ + GET_OPTIONAL_TIME_RANGE(ptr+1, 0, (Time), args->logger, "-%c", opt, optarg); \ *ptr = '\0'; \ } \ (Str) = optarg; /* Get a position followed by an optional time range */ -#define GET_POS_AND_OPTIONAL_TIME_RANGE(Dst) \ - GET_OPTIONAL_TIME_RANGE(optarg, 3, (Dst), args->logger, "-%c", opt); +#define GET_POS_AND_OPTIONAL_TIME_RANGE(Src, Dst, FullSrc) \ + GET_OPTIONAL_TIME_RANGE((Src), 3, (Dst), args->logger, "-%c", opt, (FullSrc)); res_T parse_args (const int argc, char** argv, - struct args* args) + struct args* args, + struct mem_allocator* allocator) { int opt = 0, n_used = 0; size_t len = 0; - const char option_list[] = "a:c:dD:eF:gG:hm:M:n:p:P:r:R:s:S:t:vV:"; + const char option_list[] = "a:c:dD:eF:gG:hm:M:n:p:P:r:R:s:S:t:vV:x:X:"; char buf[128]; + struct str keep; + char** line = NULL; res_T res = RES_OK; ASSERT(argv && args); + str_init(allocator, &keep); opterr = 0; /* No default error messages */ while((opt = getopt(argc, argv, option_list)) != -1) { switch (opt) { @@ -629,13 +639,13 @@ parse_args *ptr = '\0'; } if(res == RES_OK) { - if(0 == strcmp(optarg, "all")) { + if(0 == strcasecmp(optarg, "all")) { args->dump_paths = DUMP_ALL; } - else if(0 == strcmp(optarg, "error")) { + else if(0 == strcasecmp(optarg, "error")) { args->dump_paths = DUMP_ERROR; } - else if(0 == strcmp(optarg, "success")) { + else if(0 == strcasecmp(optarg, "success")) { args->dump_paths = DUMP_SUCCESS; } } @@ -752,20 +762,45 @@ parse_args goto error; } args->mode |= MODE_PROBE_COMPUTE; - GET_POS_AND_OPTIONAL_TIME_RANGE(args->pos_and_time); + GET_POS_AND_OPTIONAL_TIME_RANGE(optarg, args->pos_and_time, optarg); break; case 'P': - if(args->mode & EXCLUSIVE_MODES) { - res = RES_BAD_ARG; - logger_print(args->logger, LOG_ERROR, - "Options -%c and -%c are exclusive.\n", - (char)opt, mode_option(args->mode)); - goto error; - } - args->mode |= MODE_PROBE_COMPUTE_ON_INTERFACE; - GET_POS_AND_OPTIONAL_TIME_RANGE(args->pos_and_time); - break; + if(args->mode & EXCLUSIVE_MODES) { + res = RES_BAD_ARG; + logger_print(args->logger, LOG_ERROR, + "Options -%c and -%c are exclusive.\n", + (char)opt, mode_option(args->mode)); + goto error; + } + args->mode |= MODE_PROBE_COMPUTE_ON_INTERFACE; + + ERR(str_set(&keep, optarg)); + line = split_line(optarg, ':'); + if(!line) { + res = RES_MEM_ERR; + str_release(&keep); + goto error; + } + + /* We expect 1 or 2 parts in line */ + if(!line[0] || (line[1] && line[2])) { + logger_print((args->logger), LOG_ERROR, + "Invalid argument for option ""-%c"": %s\n", + opt, str_cget(&keep)); + str_release(&keep); + res = RES_BAD_ARG; + goto error; + } + + /* First part is pos and optional time, optional second part is a + * medium name (OK if NULL) */ + GET_POS_AND_OPTIONAL_TIME_RANGE(line[0], args->pos_and_time, + str_cget(&keep)); + if(line[1]) + args->medium_name = optarg + strlen(line[0]) + 1; + + break; case 'r': res = cstr_to_double(optarg, &args->ref_temp); @@ -856,6 +891,30 @@ parse_args goto error; } break; + + case 'x': + if(!(args->mode & RANDOM_RW_MODES)) { + res = RES_BAD_ARG; + print_multiple_modes(buf, sizeof(buf), RANDOM_RW_MODES, 0); + logger_print(args->logger, LOG_ERROR, + "Option -%c can only be used in conjunction with one of the following options: %s.\n", + (char)opt, buf); + goto error; + } + args->rndgen_state_in_filename = optarg; + break; + + case 'X': + if(!(args->mode & RANDOM_RW_MODES)) { + res = RES_BAD_ARG; + print_multiple_modes(buf, sizeof(buf), RANDOM_RW_MODES, 0); + logger_print(args->logger, LOG_ERROR, + "Option -%c can only be used in conjunction with one of the following options: %s.\n", + (char)opt, buf); + goto error; + } + args->rndgen_state_out_filename = optarg; + break; } } @@ -944,13 +1003,36 @@ parse_args } } + if(args->rndgen_state_in_filename && !(args->mode & COMPUTE_MODES)) { + res = RES_BAD_ARG; + print_multiple_modes(buf, sizeof(buf), COMPUTE_MODES, 0); + logger_print(args->logger, LOG_ERROR, + "Option -x can only be used in conjunction with an option" + " that launch a MC computation (%s).\n", + buf); + goto error; + } + + if(args->rndgen_state_out_filename && !(args->mode & COMPUTE_MODES)) { + res = RES_BAD_ARG; + print_multiple_modes(buf, sizeof(buf), COMPUTE_MODES, 0); + logger_print(args->logger, LOG_ERROR, + "Option -X can only be used in conjunction with an option" + " that launch a MC computation (%s).\n", + buf); + goto error; + } + end: + FREE_AARRAY(line); + str_release(&keep); return res; error: logger_print(args->logger, LOG_ERROR, "Use option -h to print help.\n"); goto end; } + res_T parse_camera (struct logger* logger, @@ -961,87 +1043,83 @@ parse_camera char** opt = NULL; struct camera* cam; struct str keep; + int i = 0; res_T res = RES_OK; ASSERT(cam_param && stardis); cam = &stardis->camera; line = split_line(cam_param, ':'); - str_init(stardis->allocator, &keep); + if(!line) { + res = RES_MEM_ERR; + goto error; + } - if(line) { - int i = 0; - for(i = 0; *(line + i); i++) { - size_t len = 0; - ERR(str_set(&keep, line[i])); - opt = split_line(line[i], '='); - if(!opt[0] || ! opt[1] || opt[2]) { - if(res == RES_OK) res = RES_BAD_ARG; - logger_print((logger), LOG_ERROR, - "Invalid option syntax: %s\n", str_cget(&keep)); - goto error; - } - ERR(str_set(&keep, opt[0])); - _strupr(opt[0]); - if(strcmp(opt[0], "T") == 0) { - GET_OPTIONAL_TIME_RANGE(opt[1], 0, cam->time_range, logger, "%s", - str_cget(&keep)); - } - else if(strcmp(opt[0], "FMT") == 0) { - _strupr(opt[1]); - if(strcmp(opt[1], "VTK") == 0) - cam->fmt = STARDIS_RENDERING_OUTPUT_FILE_FMT_VTK; - else if(strcmp(opt[1], "HT") == 0) - cam->fmt = STARDIS_RENDERING_OUTPUT_FILE_FMT_HT; - else { - logger_print(logger, LOG_ERROR, - "Unexpected value for rendering option %s: %s.\n", - opt[0], opt[1]); - res = RES_BAD_ARG; - goto error; - } - } - else if(strcmp(opt[0], "FOV") == 0) { - ERR(cstr_to_double(opt[1], &cam->fov)); - } - else if(strcmp(opt[0], "UP") == 0) { - ERR(cstr_to_list_double(opt[1], ',', cam->up, &len, 3)); - } - else if(strcmp(opt[0], "TGT") == 0) { - ERR(cstr_to_list_double(opt[1], ',', cam->tgt, &len, 3)); - cam->auto_look_at = 0; - } - else if(strcmp(opt[0], "POS") == 0) { - ERR(cstr_to_list_double(opt[1], ',', cam->pos, &len, 3)); - cam->auto_look_at = 0; - } - else if(strcmp(opt[0], "IMG") == 0) { - unsigned img_sz[2]; - ERR(cstr_to_list_uint(opt[1], 'x', img_sz, &len, 2)); - cam->img_width = img_sz[0]; - cam->img_height = img_sz[1]; - } - else if(strcmp(opt[0], "SPP") == 0) { - ERR(cstr_to_uint(opt[1], &cam->spp)); - } else { + str_init(stardis->allocator, &keep); + for(i = 0; *(line + i); i++) { + size_t len = 0; + ERR(str_set(&keep, line[i])); + opt = split_line(line[i], '='); + if(!opt[0] || !opt[1] || opt[2]) { /* We expect 2 parts */ + if(res == RES_OK) res = RES_BAD_ARG; + logger_print((logger), LOG_ERROR, + "Invalid option syntax: %s\n", str_cget(&keep)); + goto error; + } + if(strcasecmp(opt[0], "T") == 0) { + GET_OPTIONAL_TIME_RANGE(opt[1], 0, cam->time_range, logger, "%s", opt[0], + str_cget(&keep)); + } + else if(strcasecmp(opt[0], "FILE") == 0) { + ERR(str_set(&cam->file_name, opt[1])); + } + else if(strcasecmp(opt[0], "FMT") == 0) { + if(strcasecmp(opt[1], "VTK") == 0) + cam->fmt = STARDIS_RENDERING_OUTPUT_FILE_FMT_VTK; + else if(strcasecmp(opt[1], "HT") == 0) + cam->fmt = STARDIS_RENDERING_OUTPUT_FILE_FMT_HT; + else { logger_print(logger, LOG_ERROR, - "Unexpected option for rendering mode: %s.\n", - opt[0]); + "Unexpected value for rendering option %s: %s.\n", + opt[0], opt[1]); res = RES_BAD_ARG; goto error; } } + else if(strcasecmp(opt[0], "FOV") == 0) { + ERR(cstr_to_double(opt[1], &cam->fov)); + } + else if(strcasecmp(opt[0], "UP") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->up, &len, 3)); + } + else if(strcasecmp(opt[0], "TGT") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->tgt, &len, 3)); + cam->auto_look_at = 0; + } + else if(strcasecmp(opt[0], "POS") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->pos, &len, 3)); + cam->auto_look_at = 0; + } + else if(strcasecmp(opt[0], "IMG") == 0) { + unsigned img_sz[2]; + ERR(cstr_to_list_uint(opt[1], 'x', img_sz, &len, 2)); + cam->img_width = img_sz[0]; + cam->img_height = img_sz[1]; + } + else if(strcasecmp(opt[0], "SPP") == 0) { + ERR(cstr_to_uint(opt[1], &cam->spp)); + } else { + logger_print(logger, LOG_ERROR, + "Unexpected option for rendering mode: %s.\n", + opt[0]); + res = RES_BAD_ARG; + goto error; + } + FREE_AARRAY(opt); } end: -#define FREE_AARRAY(ARRAY) if(ARRAY) {\ - int i = 0; \ - for(i=0; *(ARRAY+i);i++){\ - free(ARRAY[i]);\ - }\ - free(ARRAY);\ -} - FREE_AARRAY(line) - FREE_AARRAY(opt) + FREE_AARRAY(line); + FREE_AARRAY(opt); #undef FREE_AARRAY str_release(&keep); @@ -1056,6 +1134,46 @@ error: #undef GET_POS_AND_OPTIONAL_TIME_RANGE #undef GET_OPTIONAL_TIME_RANGE +static res_T +description_set_name + (struct stardis* stardis, + struct str* name, + const char* tk) +{ + res_T res = RES_OK; + double foo; + const char* keywords[] = { + "AUTO", "BACK", "BOTH", "FLUID", "FRONT", "F_BOUNDARY_FOR_SOLID", + "H_BOUNDARY_FOR_FLUID", "H_BOUNDARY_FOR_SOLID", "SCALE", "SOLID", + "SOLID_FLUID_CONNECTION", "T_BOUNDARY_FOR_FLUID", "T_BOUNDARY_FOR_SOLID", + "UNKNOWN" }; + int i; + ASSERT(name && tk); + + /* Use name before uppercasing it */ + ERR(str_set(name, tk)); + + if(RES_OK == cstr_to_double(tk, &foo)) { + /* A number is not a sensible choice for a name! */ + res = RES_BAD_ARG; + goto error; + } + FOR_EACH(i, 0, sizeof(keywords) / sizeof(*keywords)) { + if(0 == strcasecmp(tk, keywords[i])) { + /* A keyword is not a sensible choice for a name! */ + res = RES_BAD_ARG; + goto error; + } + } + /* Name is OK */ + +end: + return res; +error: + logger_print(stardis->logger, LOG_ERROR, "Invalid name: %s\n", tk); + goto end; +} + static struct description* find_description_by_name (struct stardis* stardis, @@ -1102,7 +1220,7 @@ process_h desc->type = type; CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "h boundary name"); - ERR(str_set(&desc->d.h_boundary.name, tk)); + ERR(description_set_name(stardis, &desc->d.h_boundary.name, tk)); if(find_description_by_name(stardis, &desc->d.h_boundary.name, NULL) != desc) { @@ -1209,7 +1327,7 @@ process_t ? get_dummy_solid_id(stardis, dummies) : get_dummy_fluid_id(stardis, dummies); CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "temperature boundary name"); - ERR(str_set(&desc->d.t_boundary.name, tk)); + ERR(description_set_name(stardis, &desc->d.t_boundary.name, tk)); if(find_description_by_name(stardis, &desc->d.t_boundary.name, NULL) != desc) { @@ -1295,7 +1413,7 @@ process_flx desc->d.f_boundary.mat_id = get_dummy_fluid_id(stardis, dummies); CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "flux boundary name"); - ERR(str_set(&desc->d.f_boundary.name, tk)); + ERR(description_set_name(stardis, &desc->d.f_boundary.name, tk)); if(find_description_by_name(stardis, &desc->d.f_boundary.name, NULL) != desc) { @@ -1351,7 +1469,7 @@ process_sfc desc->d.sf_connect.connection_id = allocate_stardis_medium_id(stardis); CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "solid fluid connection name"); - ERR(str_set(&desc->d.sf_connect.name, tk)); + ERR(description_set_name(stardis, &desc->d.sf_connect.name, tk)); if(find_description_by_name(stardis, &desc->d.sf_connect.name, NULL) != desc) { @@ -1401,6 +1519,68 @@ error: goto end; } +/* SOLID_SOLID_CONNECTION Name contact-resitance STL_filenames */ +static res_T +process_ssc + (struct stardis* stardis, + char** tok_ctx) +{ + char* tk = NULL; + struct description* desc; + size_t sz; + res_T res = RES_OK; + + ASSERT(stardis && tok_ctx); + + stardis->counts.ssconnect_count++; + + sz = darray_descriptions_size_get(&stardis->descriptions); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_ss(stardis->allocator, &desc->d.ss_connect); + + /* Use a medium ID even if there is no medium here + * As other cases use media IDs as unique IDs for read_sides_and_files calls + * we continue the trend to ensure connection ID is OK */ + desc->type = DESC_SOLID_SOLID_CONNECT; + desc->d.ss_connect.connection_id = allocate_stardis_medium_id(stardis); + + CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "solid solid connection name"); + ERR(description_set_name(stardis, &desc->d.ss_connect.name, tk)); + if(find_description_by_name(stardis, &desc->d.ss_connect.name, NULL) + != desc) + { + logger_print(stardis->logger, LOG_ERROR, + "Name already used: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + + CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "contact resistance"); + res = cstr_to_double(tk, &desc->d.ss_connect.tcr); + if(res != RES_OK + || desc->d.ss_connect.tcr < 0) + { + logger_print(stardis->logger, LOG_ERROR, + "Invalid contact resistance: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + else if(desc->d.ss_connect.tcr == 0) { + logger_print(stardis->logger, LOG_WARNING, + "Solid-solid connection %s: defining a contact resistance to 0 has " + "no effect\n", str_cget(&desc->d.ss_connect.name)); + } + + ASSERT(sz <= UINT_MAX); + ERR(read_sides_and_files(stardis, 1, (unsigned)sz, tok_ctx)); + +end: + return res; +error: + goto end; +} + static res_T read_imposed_temperature (struct stardis* stardis, @@ -1423,8 +1603,7 @@ read_imposed_temperature } } else { /* Could be 'unknown' */ - _strupr(tk); - if(0 == strcmp(tk, "UNKNOWN")) { + if(0 == strcasecmp(tk, "UNKNOWN")) { *imposed_temperature = UNKNOWN_MEDIUM_TEMPERATURE; } else { res = RES_BAD_ARG; @@ -1462,8 +1641,7 @@ read_delta } } else { /* Could be 'auto' */ - _strupr(tk); - if(0 == strcmp(tk, "AUTO")) { + if(0 == strcasecmp(tk, "AUTO")) { /* Set to DELTA_AUTO until actual value is substituted */ *delta = DELTA_AUTO; } else { @@ -1508,7 +1686,7 @@ process_solid desc->d.solid.desc_id = (unsigned)sz; CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "solid name"); - ERR(str_set(&desc->d.solid.name, tk)); + ERR(description_set_name(stardis, &desc->d.solid.name, tk)); if(find_description_by_name(stardis, &desc->d.solid.name, NULL) != desc) { @@ -1614,7 +1792,7 @@ process_fluid desc->d.fluid.desc_id = (unsigned)sz; CHK_TOK(strtok_r(NULL, " \t", tok_ctx), "fluid name"); - ERR(str_set(&desc->d.fluid.name, tk)); + ERR(description_set_name(stardis, &desc->d.fluid.name, tk)); if(find_description_by_name(stardis, &desc->d.fluid.name, NULL) != desc) { @@ -1738,25 +1916,26 @@ process_model_line str_init(stardis->allocator, &keep); ERR(str_set(&keep, line)); CHK_TOK(strtok_r(line, " \t", &tok_ctx), "model line type"); - _strupr(tk); - if(0 == strcmp(tk, "H_BOUNDARY_FOR_SOLID")) + if(0 == strcasecmp(tk, "H_BOUNDARY_FOR_SOLID")) ERR(process_h(stardis, dummies, DESC_BOUND_H_FOR_SOLID, &tok_ctx)); - else if(0 == strcmp(tk, "H_BOUNDARY_FOR_FLUID")) + else if(0 == strcasecmp(tk, "H_BOUNDARY_FOR_FLUID")) ERR(process_h(stardis, dummies, DESC_BOUND_H_FOR_FLUID, &tok_ctx)); - else if(0 == strcmp(tk, "T_BOUNDARY_FOR_SOLID")) + else if(0 == strcasecmp(tk, "T_BOUNDARY_FOR_SOLID")) ERR(process_t(stardis, dummies, DESC_BOUND_T_FOR_SOLID, &tok_ctx)); - else if(0 == strcmp(tk, "T_BOUNDARY_FOR_FLUID")) + else if(0 == strcasecmp(tk, "T_BOUNDARY_FOR_FLUID")) ERR(process_t(stardis, dummies, DESC_BOUND_T_FOR_FLUID, &tok_ctx)); - else if(0 == strcmp(tk, "F_BOUNDARY_FOR_SOLID")) + else if(0 == strcasecmp(tk, "F_BOUNDARY_FOR_SOLID")) ERR(process_flx(stardis, dummies, &tok_ctx)); - else if(0 == strcmp(tk, "SOLID_FLUID_CONNECTION")) + else if(0 == strcasecmp(tk, "SOLID_FLUID_CONNECTION")) ERR(process_sfc(stardis, &tok_ctx)); - else if(0 == strcmp(tk, "SOLID")) + else if(0 == strcasecmp(tk, "SOLID_SOLID_CONNECTION")) + ERR(process_ssc(stardis, &tok_ctx)); + else if(0 == strcasecmp(tk, "SOLID")) ERR(process_solid(stardis, &tok_ctx)); - else if(0 == strcmp(tk, "FLUID")) + else if(0 == strcasecmp(tk, "FLUID")) ERR(process_fluid(stardis, &tok_ctx)); - else if(0 == strcmp(tk, "SCALE")) + else if(0 == strcasecmp(tk, "SCALE")) ERR(process_scale(stardis, &tok_ctx)); else { logger_print(stardis->logger, LOG_ERROR, diff --git a/src/stardis-parsing.h b/src/stardis-parsing.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 @@ -81,7 +81,11 @@ enum stardis_mode { USE_STDOUT_MODES = MODE_DUMP_C_CHUNKS | MODE_DUMP_VTK | MODE_DUMP_HELP | MODE_DUMP_VERSION - | MODE_IR_COMPUTE | MODE_GREEN + | MODE_IR_COMPUTE | MODE_GREEN, + + RANDOM_RW_MODES + = MODE_PROBE_COMPUTE | MODE_PROBE_COMPUTE_ON_INTERFACE | MODE_MEDIUM_COMPUTE + | MODE_BOUNDARY_COMPUTE | MODE_FLUX_BOUNDARY_COMPUTE }; STATIC_ASSERT(GREEN_COMPATIBLE_MODES == (COMPUTE_MODES & GREEN_COMPATIBLE_MODES), @@ -103,6 +107,8 @@ struct args { char* bin_green_filename; char* end_paths_filename; char* paths_filename; + char* rndgen_state_in_filename; + char* rndgen_state_out_filename; char* chunks_prefix; size_t samples; unsigned nthreads; @@ -175,7 +181,8 @@ extern res_T parse_args (const int argc, char** argv, - struct args* args); + struct args* args, + struct mem_allocator* allocator); extern res_T parse_camera diff --git a/src/stardis-solid.c b/src/stardis-solid.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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 diff --git a/src/stardis-solid.h b/src/stardis-solid.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2018-2020 |Meso|Star> (contact@meso-star.com) +/* Copyright (C) 2018-2021 |Meso|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