commit 9179878168be78ab80fa39bd16532345557ab760
parent 931d1ca82ae4bbc34a39072d2f707c9c322ac4b2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 10 Oct 2025 16:17:20 +0200
stardis: test Tsrf taking into account the sun
The script calculates the surface temperature over several days and now
performs the calculation in two ways: the previous method treats the sun
as a flux imposed on the ground, and the new method defines the sun as
an external source whose position, power, and diffuse radiation depend
on the meteorological data.
The plugin's data creation is enhanced by arguments to allow or
prohibit the use of meteorological data, depending on the desired
connection or boundary condition to be defined. For example, describing
the sun as an external source prohibits the use of a Neumann+Robin
boundary condition for the ground. The imposed flux must be disabled and
the interface becomes a connection with the atmosphere whose temperature
is known. New plugin functions are thus provided to describe the
atmosphere.
Finally, the graph of the test results now compares:
- the ground temperature provided in the input meteorological data,
- the ground temperature computed by Monte Carlo using an imposed flux
to represent the contribution of the sun,
- and the ground temperature, also calculated by Monte Carlo but
modeling the sun as an external source.
Diffstat:
5 files changed, 168 insertions(+), 38 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -11,5 +11,5 @@ test*
test.txt
smeteo
*.stl
-stardis_model.txt
+stardis_model*.txt
Tsrf.*
diff --git a/src/plot.gp b/src/plot.gp
@@ -33,14 +33,16 @@ set xtics 86400 # 1 day in seconds
# Define the plot style
set style line 1 lc "black" lw 2 # GCM
-set style line 2 lc "black" pt 6 ps 1 # MC
+set style line 2 lc "black" pt 6 ps 1 # MC Flux
+set style line 3 lc "black" pt 8 ps 1 # MC Sun
set errorbars 3 # Width of the error bars
# Move the labels above the graph
-set key at graph 1, 1.15
+set key at graph 1, 1.17
# Increase the right margin to prevent the last date from being cut off
set rmargin 6
plot Tsrf u 1:3 title "GCM" w line ls 1, \
- Tsrf u 1:4:($5*3) title "MC" w errorbars ls 2
+ Tsrf u 1:4:5 title "MC Flux" w errorbars ls 2, \
+ Tsrf u 1:6:7 title "MC Sun" w errorbars ls 3
diff --git a/src/stardis_smeteo.c b/src/stardis_smeteo.c
@@ -22,6 +22,7 @@
#include <star/scem.h>
#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
#include <rsys/math.h>
#include <rsys/mem_allocator.h>
@@ -30,29 +31,70 @@
#include <limits.h>
#include <time.h> /* localtime_r */
+#include <unistd.h> /* getopt */
#define EARTH_TO_SUN_DST 149.5978707e9 /* [m] */
+struct args {
+ double temperature; /* Take precedence on meteorological data [K] */
+ int convection; /* Enable convection */
+ int imposed_flux; /* Enable imposed flux */
+};
+#define ARGS_DEFAULT__ {STARDIS_TEMPERATURE_NONE, 0, 0}
+static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
+
struct boundary_condition {
+ struct args args;
struct stardis_smeteo_lib* lib;
struct stardis_smeteo_lib_desc lib_desc;
};
-
-struct args {
- char dummy; /* No argument yet */
+static const struct boundary_condition BOUNDARY_CONDITION_NULL = {
+ ARGS_DEFAULT__, NULL, STARDIS_SMETEO_LIB_DESC_NULL__
};
-#define ARGS_DEFAULT__ {0}
-static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
/*******************************************************************************
* Helper functions
******************************************************************************/
+static void
+usage(FILE* stream, const char* name)
+{
+ ASSERT(stream && name);
+ fprintf(stream, "usage: %s [-fh] [-t temperature]\n", name);
+
+}
+
static res_T
args_init(struct args* args, int argc, char* argv[])
{
- (void)argc, (void)argv;
+ int opt = 0;
+ res_T res = RES_OK;
+ ASSERT(args && argc && argv);
+
*args = ARGS_DEFAULT;
- return RES_OK;
+
+ optind = 1;
+ while((opt=getopt(argc, argv, "hft:")) != -1) {
+ switch(opt) {
+ case 'h': args->convection = 1; break;
+ case 'f': args->imposed_flux = 1; break;
+ case 't':
+ res = cstr_to_double(optarg, &args->temperature);
+ break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+exit:
+ return res;
+error:
+ usage(stderr, argv[0]);
+ goto exit;
}
static int
@@ -206,15 +248,15 @@ stardis_create_data
size_t argc,
char* argv[])
{
- struct args args = ARGS_DEFAULT;
struct boundary_condition* bcond = NULL;
res_T res = RES_OK;
if(!ctx || !lib_data || !argc || !argv) goto error;
- if((res = args_init(&args, (int)argc, argv)) != RES_OK) goto error;
if(!(bcond = mem_calloc(1, sizeof(*bcond)))) goto error;
+ *bcond = BOUNDARY_CONDITION_NULL;
+ if((res = args_init(&bcond->args, (int)argc, argv)) != RES_OK) goto error;
bcond->lib = lib_data;
stardis_smeteo_lib_ref_get(bcond->lib);
stardis_smeteo_lib_get_desc(bcond->lib, &bcond->lib_desc);
@@ -242,11 +284,17 @@ stardis_convection_coefficient
void* data)
{
struct boundary_condition* bcond = data;
- size_t i = 0; /* Index of the meteo entry including fragment time */
ASSERT(frag && data);
- i = get_meteo_entry_id(bcond, frag->time);
- return bcond->lib_desc.smeteo_desc.entries[i].H;
+ if(!bcond->args.convection) {
+ return 0;
+
+ } else {
+ size_t i = 0; /* Index of the meteo entry including fragment time */
+
+ i = get_meteo_entry_id(bcond, frag->time);
+ return bcond->lib_desc.smeteo_desc.entries[i].H;
+ }
}
double /* [W/K/m^2] */
@@ -263,21 +311,28 @@ stardis_boundary_flux
void* data)
{
struct boundary_condition* bcond = data;
- double net_flux = 0; /* [W/m^2] */
- double SWdn = 0; /* shortwave downward flux [W/m^2] */
- size_t i = 0;
ASSERT(frag && data);
- i = get_meteo_entry_id(bcond, frag->time);
+ if(!bcond->args.imposed_flux) {
+ return STARDIS_FLUX_NONE;
+
+ } else {
+ double net_flux = 0; /* [W/m^2] */
+ double SWdn = 0; /* shortwave downward flux [W/m^2] */
+ size_t i = 0;
- SWdn = bcond->lib_desc.smeteo_desc.entries[i].SWdn_direct
- + bcond->lib_desc.smeteo_desc.entries[i].SWdn_diffuse;
+ i = get_meteo_entry_id(bcond, frag->time);
- /* Net flux on the ground side, i.e. the downward flux - upward flux */
- net_flux = SWdn
- - bcond->lib_desc.smeteo_desc.entries[i].SWup
- - bcond->lib_desc.smeteo_desc.entries[i].LE;
- return net_flux;
+ SWdn = bcond->lib_desc.smeteo_desc.entries[i].SWdn_direct
+ + bcond->lib_desc.smeteo_desc.entries[i].SWdn_diffuse;
+
+ /* Net flux on the ground side, i.e. the downward flux - upward flux */
+ net_flux = SWdn
+ - bcond->lib_desc.smeteo_desc.entries[i].SWup
+ - bcond->lib_desc.smeteo_desc.entries[i].LE;
+
+ return net_flux;
+ }
}
double /* [K] */
@@ -330,6 +385,20 @@ stardis_reference_temperature
return stardis_boundary_temperature(frag, data);
}
+double
+stardis_calorific_capacity(const struct stardis_vertex* vtx, void* data)
+{
+ (void)vtx, (void)data;
+ return 1; /* Dummy */
+}
+
+double
+stardis_volumic_mass(const struct stardis_vertex* vtx, void* data)
+{
+ (void)vtx, (void)data;
+ return 1; /* Dummy */
+}
+
/* Ground temperature */
double /* [K] */
stardis_boundary_temperature
@@ -337,11 +406,16 @@ stardis_boundary_temperature
void* data)
{
struct boundary_condition* bcond = data;
- size_t i = 0; /* Index of the meteo entry including fragment time */
ASSERT(frag && data);
- i = get_meteo_entry_id(bcond, frag->time);
- return bcond->lib_desc.smeteo_desc.entries[i].Tsrf;
+ if(STARDIS_TEMPERATURE_IS_KNOWN(bcond->args.temperature)) {
+ return bcond->args.temperature;
+ } else {
+ size_t i = 0; /* Index of the meteo entry including fragment time */
+
+ i = get_meteo_entry_id(bcond, frag->time);
+ return bcond->lib_desc.smeteo_desc.entries[i].Tsrf;
+ }
}
double /* [K] */
diff --git a/src/stardis_smeteo.h b/src/stardis_smeteo.h
@@ -119,7 +119,6 @@ stardis_reference_temperature
(const struct stardis_interface_fragment* frag,
void* data);
-
/* Radiative temperature obtained from meteorological data. This is the to
* calculate radiative exchanges between the ground and the sky */
STARDIS_API double
@@ -143,6 +142,29 @@ stardis_t_range
double range[2]); /* >0 [K] */
/*******************************************************************************
+ * Functions for the sky fluid.
+ * They are used to determine the sky temperature from meteorological data and
+ * enable the evaluation of radiative transfer using the Monte Carlo method.
+ ******************************************************************************/
+/* Dummy function since fluid as fixed temperature */
+STARDIS_API double
+stardis_calorific_capacity
+ (const struct stardis_vertex* vtx,
+ void* data);
+
+/* Dummy function since fluid as fixed temperature */
+STARDIS_API double
+stardis_volumic_mass
+ (const struct stardis_vertex* vtx,
+ void* data);
+
+/* Atmospheric temperature retrieve from meteorological data */
+STARDIS_API double /* >0 [K] */
+stardis_medium_temperature
+ (const struct stardis_vertex* vtx,
+ void* data);
+
+/*******************************************************************************
* Functions for a Dirichlet boundary condition.
* This condition is used as a very basic coupling between Stardis and the
* meteorological data, of which only its ground temperature is used.
diff --git a/src/test_stardis_smeteo_ground_temperature.sh b/src/test_stardis_smeteo_ground_temperature.sh
@@ -101,7 +101,7 @@ ground() # list of facets id
printf 'endsolid\n'
}
-stardis_input()
+stardis_input_imposed_flux()
{
# The plugin
printf 'PROGRAM Meteo libstardis_smeteo.so %s\n' "${smeteo}"
@@ -114,10 +114,39 @@ stardis_input()
# Limit condition
echo 'H_BOUNDARY_FOR_SOLID adiabatic 0 0 0 0 0 ground_xXyY.stl'
echo 'T_BOUNDARY_FOR_SOLID underground 284.0 ground_z.stl'
- echo 'HF_BOUNDARY_FOR_SOLID_PROG atmosphere Meteo ground_Z.stl'
+ echo 'HF_BOUNDARY_FOR_SOLID_PROG atm Meteo ground_Z.stl PROG_PARAMS -hf'
echo 'TRAD_PROG Meteo'
}
+stardis_input_external_source()
+{
+ # The plugin
+ printf 'PROGRAM Meteo libstardis_smeteo.so %s\n' "${smeteo}"
+ echo ''
+
+ # Media
+ echo 'SOLID ground 1 1500 1500 0.02 281.85 UNKNOWN 0 FRONT ground.stl'
+ echo 'FLUID_PROG sky Meteo BACK ground.stl'
+ echo ''
+
+ # Connection
+ echo 'SOLID_FLUID_CONNECTION_PROG ground_sky_connect Meteo ground_Z.stl PROG_PARAMS -h'
+ echo 'SOLID_FLUID_CONNECTION adiabatic 300 0 0 0 ground_xXyY.stl'
+ echo ''
+
+ # Limit condition
+ echo 'SOLID_FLUID_CONNECTION_PROG underground Meteo ground_z.stl PROG_PARAMS -t 284.0'
+ echo ''
+ echo 'SPHERICAL_SOURCE_PROG 6.987e8 Meteo'
+ echo ''
+ echo 'TRAD_PROG Meteo'
+}
+
+run_stardis() # model, time
+{
+ stardis -a wos -M "$1" -s ground_Z.stl,"$2"
+}
+
########################################################################
# The test
########################################################################
@@ -128,10 +157,11 @@ ground 3 4 > ground_Z.stl
ground 5 6 7 8 9 10 11 12 > ground_xXyY.stl
# Generate the Stardis input file
-stardis_input > stardis_model.txt
+stardis_input_imposed_flux > stardis_model_imposed_flux.txt
+stardis_input_external_source > stardis_model_external_source.txt
# Use Stardis to calculate the ground temperature
-date="28-AUG-1850 01:30:00" # 1st date on which temperature is calculated
+date="30-AUG-1850 01:30:00" # 1st date on which temperature is calculated
ndates="32" # Overall number of dates
# Check that at least the 1st date exists in the meteorlogical file
@@ -174,11 +204,13 @@ printf 'Calculating 0 %%\n'
# Run stardis to retrieve the expected value and standard error of the
# ground temperature estimate.
- mc="$(stardis -a wos -M stardis_model.txt -s ground_Z.stl,"${time}" \
- | cut -d' ' -f1-2)"
+ mc0="$(run_stardis stardis_model_imposed_flux.txt "${time}" \
+ | cut -d' ' -f1-2)"
+ mc1="$(run_stardis stardis_model_external_source.txt "${time}" \
+ | cut -d' ' -f1-2)"
# Print the calculation result
- printf '%s %s %s\n' "${date}" "${Tsrf}" "${mc}"
+ printf '%s %s %s %s\n' "${date}" "${Tsrf}" "${mc0}" "${mc1}"
# Print progression bar
isimu="$((isimu+1))"
