commit 393e031f17f861869248a325bbca9ff66f3ed693
parent a0b4b834351f769283a221fca9be347972e03c03
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 28 Mar 2024 18:21:18 +0100
Add support for green function estimation to Wos
When estimating the green function, the WoS conduction algorithm now
records the power term if there is one, and saves the path position when
the initial condition is reached.
The test_sdis_volumic_power was then adjusted to test the WoS algorithm
in addition to the delta sphere algorithm, in both 2D and 3D with and
without green. 4 tests are performed in 2D and 3D. 2 in steady state and
2 in unsteady state. Half of these tests are now carried out with WoS to
test this algorithm without increasing the overall execution time: the
total number of tests remains the same.
Diffstat:
2 files changed, 128 insertions(+), 96 deletions(-)
diff --git a/src/sdis_heat_path_conductive_wos_Xd.h b/src/sdis_heat_path_conductive_wos_Xd.h
@@ -483,23 +483,62 @@ XD(handle_volumic_power_wos)
(struct sdis_scene* scn,
const struct solid_props* props,
const double distance, /* [m/fp_to_meter] */
+ double* power_term,
struct XD(temperature)* T)
{
double dst = distance * scn->fp_to_meter; /* [m] */
+ double term = 0;
res_T res = RES_OK;
- ASSERT(scn && props && distance >= 0 && T);
+ ASSERT(scn && props && distance >= 0 && power_term && T);
if(props->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
/* No displacement => no power density */
if(distance == 0) goto exit;
- T->value += props->power * dst*dst / (2*DIM*props->lambda);
+ term = dst*dst / (2*DIM*props->lambda);
+ T->value += props->power * term;
exit:
+ *power_term = term;
return res;
}
+static res_T
+XD(update_green_path)
+ (struct green_path_handle* green_path,
+ struct XD(rwalk)* rwalk,
+ struct sdis_medium* mdm,
+ const struct solid_props* props,
+ const double power_term,
+ const struct XD(temperature)* T)
+{
+ res_T res = RES_OK;
+ ASSERT(mdm && props && T);
+
+ /* Is the green function estimated? */
+ if(!green_path) goto exit;
+
+ /* Save power term for green function if any */
+ if(props->power != SDIS_VOLUMIC_POWER_NONE) {
+ res = green_path_add_power_term(green_path, mdm, &rwalk->vtx, power_term);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Set the green path limit to the current position if the initial condition
+ * has been reached */
+ if(T->done) {
+ res = green_path_set_limit_vertex
+ (green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local function
******************************************************************************/
@@ -512,13 +551,15 @@ XD(conductive_path_wos)
struct XD(temperature)* T)
{
/* Properties */
+ struct sdis_medium* mdm = NULL;
struct solid_props props_ref = SOLID_PROPS_NULL;
struct solid_props props = SOLID_PROPS_NULL;
double alpha = 0; /* diffusivity, i.e. lambda/(rho*cp) */
/* Miscellaneous */
size_t ndiffusion_steps = 0; /* For debug */
- const int green = 0;
+ double green_power_term = 0;
+ int green = 0;
const int wos = 1;
res_T res = RES_OK;
(void)ctx; /* Avoid the "unused variable" warning */
@@ -527,6 +568,14 @@ XD(conductive_path_wos)
ASSERT(scn && ctx && rwalk && rng && T);
ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_SOLID);
+ /* Is green evaluated evaluated */
+ green = ctx->green_path != NULL;
+
+ /* Keep track of the solid. After conduction, a boundary may have been
+ * reached, so the random walk medium is NULL. However, this medium is still
+ * needed to update the green path. Hence this backup */
+ mdm = rwalk->mdm;
+
res = XD(check_medium_consistency)(scn, rwalk);
if(res != RES_OK) goto error;
@@ -548,6 +597,7 @@ XD(conductive_path_wos)
/* Sample a diffusive path */
for(;;) {
+ double power_term = 0; /* */
double dst = 0; /* [m/fp_to_meter] */
/* Register the new vertex against the heat path */
@@ -574,9 +624,12 @@ XD(conductive_path_wos)
if(res != RES_OK) goto error;
/* Add the volumic power density */
- res = XD(handle_volumic_power_wos)(scn, &props, dst, T);
+ res = XD(handle_volumic_power_wos)(scn, &props, dst, &power_term, T);
if(res != RES_OK) goto error;
+ /* Accumulate the power term */
+ if(green) green_power_term += power_term;
+
/* The path reaches the initial condition */
if(T->done) {
REGISTER_HEAT_VERTEX;
@@ -603,6 +656,10 @@ XD(conductive_path_wos)
++ndiffusion_steps; /* For debug */
}
+ /* Save green function data */
+ res = XD(update_green_path)
+ (ctx->green_path, rwalk, mdm, &props_ref, green_power_term, T);
+
exit:
return res;
error:
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -53,6 +53,22 @@
#define DELTA 1.0/55.0
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static const char*
+algo_cstr(const enum sdis_diffusion_algorithm diff_algo)
+{
+ const char* cstr = "none";
+
+ switch(diff_algo) {
+ case SDIS_DIFFUSION_DELTA_SPHERE: cstr = "delta sphere"; break;
+ case SDIS_DIFFUSION_WOS: cstr = "WoS"; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return cstr;
+}
+
+/*******************************************************************************
* Media
******************************************************************************/
struct solid {
@@ -173,7 +189,6 @@ solve
struct sdis_mc time = SDIS_MC_NULL;
size_t nreals;
size_t nfails;
- double ref = SDIS_TEMPERATURE_NONE;
enum sdis_scene_dimension dim;
const int nsimuls = 4;
int isimul;
@@ -181,11 +196,16 @@ solve
OK(sdis_scene_get_dimension(scn, &dim));
FOR_EACH(isimul, 0, nsimuls) {
- int steady = (isimul % 2) == 0;
+ const enum sdis_diffusion_algorithm algo = (isimul / 2) % 2
+ ? SDIS_DIFFUSION_WOS
+ : SDIS_DIFFUSION_DELTA_SPHERE;
+ const int steady = (isimul % 2) == 0;
+ double power = P0 == SDIS_VOLUMIC_POWER_NONE ? 0 : P0;
/* Restore power value */
solid->vpower = P0;
+ solve_args.diff_algo = algo;
solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9);
solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9);
solve_args.position[2] =
@@ -209,42 +229,28 @@ solve
OK(sdis_estimator_get_temperature(estimator, &T));
OK(sdis_estimator_get_realisation_time(estimator, &time));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s "
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ } else {
+ printf(
+ "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
printf("Elapsed time = %s\n", dump);
printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- if(steady) CHK(eq_eps(T.E, ref, T.SE * 3));
/* Check green function */
time_current(&t0);
@@ -257,37 +263,21 @@ solve
OK(sdis_estimator_get_failure_count(estimator2, &nfails));
OK(sdis_estimator_get_temperature(estimator2, &T));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Green Steady temperature at (%g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Green Mean temperature at (%g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Green Steady temperature at (%g, %g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Green Mean temperature at (%g, %g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Green steady temperature - pos: %g, %g, %g; Power: %g algo: %s"
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ } else {
+ printf(
+ "Green mean temperature - pos: %g, %g, %g; t: [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
@@ -305,7 +295,7 @@ solve
printf("\n");
/* Check same green used at a different power level */
- solid->vpower = 3 * P0;
+ solid->vpower = power = 3 * P0;
time_current(&t0);
OK(sdis_solve_probe(scn, &solve_args, &estimator));
@@ -317,43 +307,28 @@ solve
OK(sdis_estimator_get_temperature(estimator, &T));
OK(sdis_estimator_get_realisation_time(estimator, &time));
- switch(dim) {
- case SDIS_SCENE_2D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g) with Power=%g = %g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g) with t in [%g %g] and Power=%g"
- " ~ %g +/- %g\n",
- SPLIT2(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- case SDIS_SCENE_3D:
- if(steady) {
- double x = solve_args.position[0] - 0.5;
- ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
- printf("Steady temperature at (%g, %g, %g) with Power=%g = %g"
- " ~ %g +/- %g\n",
- SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE);
- } else {
- printf("Mean temperature at (%g, %g, %g) with t in [%g %g] and"
- " Power=%g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
- solid->vpower, T.E, T.SE);
- }
- break;
- default: FATAL("Unreachable code.\n"); break;
+ if(steady) {
+ const double x = solve_args.position[0] - 0.5;
+ const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+ printf
+ ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s "
+ "= %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ } else {
+ printf(
+ "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s "
+ "~ %g +/- %g\n",
+ SPLIT3(solve_args.position), SPLIT2(solve_args.time_range),
+ power, algo_cstr(algo), T.E, T.SE);
}
+
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
printf("Elapsed time = %s\n", dump);
printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
CHK(nfails + nreals == N);
CHK(nfails <= N/1000);
- if(steady) CHK(eq_eps(T.E, ref, T.SE * 3));
time_current(&t0);
OK(sdis_green_function_solve(green, &estimator2));
