commit 1dfc9e764907a00b25cfc4eabf909948d8161eed
parent e6e567e2a6c398c6e5b65f9fb2c1356b78949d8d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 9 Feb 2024 11:19:05 +0100
Making the power of external sources time dependent
This makes it possible to simulate, for example, how solar radiation is
attenuated by the atmosphere (it is more attenuated on grazing angles).
Diffstat:
9 files changed, 90 insertions(+), 40 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -140,19 +140,25 @@ struct sdis_info {
#define SDIS_INFO_NULL__ {0}
static const struct sdis_info SDIS_INFO_NULL = SDIS_INFO_NULL__;
-/* Type of functor used to retrieve the source's position relative to time. */
+/* Type of functor used to retrieve the source's position relative to time */
typedef void
(*sdis_get_position_T)
(const double time,
double pos[3],
struct sdis_data* data);
+/* Type of functor used to retrieve the source's power relative to time */
+typedef double
+(*sdis_get_power_T)
+ (const double time,
+ struct sdis_data* data);
+
/* Input arguments of the sdis_spherical_source_create function */
struct sdis_spherical_source_create_args {
sdis_get_position_T position; /* [m] */
+ sdis_get_power_T power; /* Total power [W] */
struct sdis_data* data; /* Data sent to the position functor */
double radius; /* [m] */
- double power; /* Total power [W] */
};
#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, 0, 0}
static const struct sdis_spherical_source_create_args
@@ -999,10 +1005,10 @@ SDIS_API res_T
sdis_source_ref_put
(struct sdis_source* source);
-SDIS_API res_T
+SDIS_API double
sdis_source_get_power
- (const struct sdis_source* src,
- double* power);/* [W] */
+ (struct sdis_source* source,
+ const double time); /* [s] */
/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
@@ -1526,7 +1532,7 @@ sdis_solve_probe_boundary_list
*
* Also note that the green solvers assume that the interface fluxes are
* constant in time and space. The same applies to the volumic power of the
- * solid media.
+ * solid media and the power of external sources.
*
* If these assumptions are not ensured by the caller, the behavior of the
* estimated green function is undefined.
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -453,8 +453,10 @@ green_function_solve_path
/* Compute external flux */
external_flux = 0;
if(green->scn->source) {
+ /* NOTE: The power of the source is assumed to be constant over time and is
+ * therefore recovered at steady state */
external_flux =
- path->external_flux_term * source_get_power(green->scn->source);
+ path->external_flux_term * source_get_power(green->scn->source, INF);
}
/* Compute path's end temperature */
diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h
@@ -505,7 +505,7 @@ XD(handle_external_net_flux)
external_flux_term = net_flux / (args->h_radi + args->h_conv + args->h_cond);
/* Update the Monte Carlo weight */
- T->value += external_flux_term * source_get_power(scn->source);
+ T->value += external_flux_term * source_get_power(scn->source, frag.time);
/* Register the external net flux term */
if(args->green_path) {
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -47,15 +47,14 @@ check_spherical_source_create_args
return RES_BAD_ARG;
}
- if(args->radius < 0) {
- log_err(dev, "%s: invalid source radius '%g' m. It cannot be negative.\n",
- func_name, args->radius);
+ if(!args->power) {
+ log_err(dev, "%s: the power functor is missing.\n", func_name);
return RES_BAD_ARG;
}
- if(args->power < 0) {
- log_err(dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
- func_name, args->power);
+ if(args->radius < 0) {
+ log_err(dev, "%s: invalid source radius '%g' m. It cannot be negative.\n",
+ func_name, args->radius);
return RES_BAD_ARG;
}
@@ -126,12 +125,11 @@ sdis_source_ref_put(struct sdis_source* src)
return RES_OK;
}
-res_T
-sdis_source_get_power(const struct sdis_source* src, double* power)
+double
+sdis_source_get_power(struct sdis_source* src, const double time /* [s] */)
{
- if(!src || !power) return RES_BAD_ARG;
- *power = source_get_power(src);
- return RES_OK;
+ ASSERT(src);
+ return source_get_power(src, time);
}
/*******************************************************************************
@@ -151,14 +149,23 @@ source_sample
double cos_half_angle; /* [radians] */
double dst; /* [m] */
double radius; /* Source radius [m] */
+ double power; /* Source power [W] */
double area; /* Source area [m^2] */
res_T res = RES_OK;
ASSERT(src && rng && pos && sample);
- /* Retrieve current source position and radius */
+ /* Retrieve current source position, radius and power */
src->spherical.position(time, src_pos, src->spherical.data);
+ power = src->spherical.power(time, src->spherical.data);
radius = src->spherical.radius;
+ if(power < 0) {
+ log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ FUNC_NAME, power);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
area = 4*PI*radius*radius; /* [m^2] */
/* compute the direction of `pos' toward the center of the source */
@@ -186,7 +193,7 @@ source_sample
d3_set(sample->dir, main_dir);
sample->pdf = 1;
sample->dst = dst;
- sample->power = src->spherical.power; /* [W] */
+ sample->power = power; /* [W] */
sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */
sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
@@ -201,7 +208,7 @@ source_sample
/* Set other sample variables */
sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->power = src->spherical.power; /* [W] */
+ sample->power = power; /* [W] */
sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
@@ -223,6 +230,7 @@ source_trace_to
double src_pos[3]; /* [m] */
double main_dir[3];
double radius; /* [m] */
+ double power; /* [W] */
double dst; /* Distance from pos to the source center [m] */
double half_angle; /* [radian] */
res_T res = RES_OK;
@@ -237,8 +245,16 @@ source_trace_to
goto exit;
}
- /* Retrieve current source position and radius */
+ /* Retrieve current source position and power */
src->spherical.position(time, src_pos, src->spherical.data);
+ power = src->spherical.power(time, src->spherical.data);
+
+ if(power < 0) {
+ log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ FUNC_NAME, power);
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* compute the direction of `pos' toward the center of the source */
d3_sub(main_dir, src_pos, pos);
@@ -274,7 +290,7 @@ source_trace_to
d3_set(sample->dir, dir);
sample->pdf = 1;
sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->power = src->spherical.power; /* [W] */
+ sample->power = power; /* [W] */
sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
@@ -286,11 +302,11 @@ error:
goto exit;
}
-double
-source_get_power(const struct sdis_source* src)
+double /* [W] */
+source_get_power(const struct sdis_source* src, const double time /* [s] */)
{
ASSERT(src);
- return src->spherical.power;
+ return src->spherical.power(time, src->spherical.data);
}
void
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -63,7 +63,8 @@ source_trace_to
extern LOCAL_SYM double /* [W] */
source_get_power
- (const struct sdis_source* source);
+ (const struct sdis_source* source,
+ const double time); /* [s] */
extern LOCAL_SYM void
source_compute_signature
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -182,6 +182,13 @@ source_get_position
pos[2] = 5.0;
}
+static double
+source_get_power(const double time, struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ return SOURCE_POWER; /* [W] */
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -189,9 +196,9 @@ create_source(struct sdis_device* sdis)
struct sdis_source* source = NULL;
args.position = source_get_position;
+ args.power = source_get_power;
args.data = NULL;
args.radius = 3e-1; /* [m] */
- args.power = SOURCE_POWER; /* [W] */
OK(sdis_spherical_source_create(sdis, &args, &source));
return source;
}
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -261,6 +261,13 @@ source_get_position
pos[2] = 0;
}
+static double
+source_get_power(const double time, struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ return 3.845e26; /* [W] */
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -268,9 +275,9 @@ create_source(struct sdis_device* sdis)
struct sdis_source* src = NULL;
args.position = source_get_position;
+ args.power = source_get_power;
args.data = NULL;
args.radius = 6.5991756e8; /* [m] */
- args.power = 3.845e26; /* [W] */
OK(sdis_spherical_source_create(sdis, &args, &src));
return src;
}
diff --git a/src/test_sdis_source.c b/src/test_sdis_source.c
@@ -26,7 +26,16 @@ spherical_source_get_position
struct sdis_data* data)
{
(void)time, (void)data;
- pos[0] = pos[1] = pos[2] = 1.234;
+ pos[0] = pos[1] = pos[2] = 1.234; /* [m] */
+}
+
+static double
+spherical_source_get_power
+ (const double time,
+ struct sdis_data* data)
+{
+ (void)time, (void)data;
+ return 10; /* [W] */
}
static void
@@ -36,25 +45,21 @@ check_spherical_source(struct sdis_device* dev)
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
struct sdis_source* src = NULL;
struct sdis_data* data = NULL;
- double power = 0;
/* Create a data to check its memory management */
OK(sdis_data_create(dev, sizeof(double[3]), ALIGNOF(double[3]), NULL, &data));
args.position = spherical_source_get_position;
+ args.power = spherical_source_get_power;
args.data = data;
args.radius = 1;
- args.power = 10;
BA(sdis_spherical_source_create(NULL, &args, &src));
BA(sdis_spherical_source_create(dev, NULL, &src));
BA(sdis_spherical_source_create(dev, &args, NULL));
OK(sdis_spherical_source_create(dev, &args, &src));
- BA(sdis_source_get_power(NULL, &power));
- BA(sdis_source_get_power(src, NULL));
- OK(sdis_source_get_power(src, &power));
- CHK(power == args.power);
+ CHK(sdis_source_get_power(src, INF) == 10);
BA(sdis_source_ref_get(NULL));
OK(sdis_source_ref_get(src));
@@ -67,6 +72,12 @@ check_spherical_source(struct sdis_device* dev)
args.data = NULL;
OK(sdis_spherical_source_create(dev, &args, &src));
OK(sdis_source_ref_put(src));
+
+ args.position = NULL;
+ BA(sdis_spherical_source_create(dev, &args, &src));
+ args.position = spherical_source_get_position;
+ args.power = NULL;
+ BA(sdis_spherical_source_create(dev, &args, &src));
}
/*******************************************************************************
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -99,7 +99,6 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
struct sdis_data* data = NULL;
enum sdis_medium_type type;
enum sdis_green_path_end_type end_type;
- double source_power = 0; /* [W] */
double power = 0;
double flux = 0;
double external_flux = 0; /* [W/m^2] */
@@ -139,8 +138,9 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
if(source == NULL) {
CHK(external_flux == 0);
} else {
- OK(sdis_source_get_power(source, &source_power));
- external_flux *= source_power;
+ /* NOTE: source power is assumed constant in time and is therefore retrieved
+ * at steady state*/
+ external_flux *= sdis_source_get_power(source, INF); /* [W] */
}
BA(sdis_green_path_get_end_type(NULL, NULL));
