stardis-solver

commit 4f8944d8bb8fe3415b594d59b68b4f1ac0060114
parent 7647c77bb0c6f6faa2f9b36d69ba539a38f247b1
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  3 Mar 2023 12:05:23 +0100

Change the value of the constant SDIS_VOLUMIC_POWER_NONE

Let be a solid medium with a 'volumic_power' functor defined and whose
value returned is zero. The library assumed that there was no volumic
power and therefore simply ignored it. This is fine except for the
propagator estimation: ignoring the medium's volumic power means that
the medium is not registered in the list of mediums with a volumic
power. A then non-zero volumic power was therefore not taken into
account when the propagator was re-evaluated.

This commit corrects this problem by differentiating between a medium
with a volumic power of 0, and a medium with no volumic power at all.
SDIS_VOLUMIC_POWER_NONE is now set to DBL_MAX which means that the
medium has no volumic power, while a value of 0 is now treated as any
valid volumic power term. In fact, volumic powers and imposed fluxes
are now treated in the same way.

Diffstat:
M
src/sdis.h
 | 
2
+-
M
src/test_sdis_volumic_power.c
 | 
43
++++++++++++++++++++++++++++++++++++++++++-

2 files changed, 43 insertions(+), 2 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -43,7 +43,7 @@
  * as CPU cores */
 #define SDIS_NTHREADS_DEFAULT (~0u)
 
-#define SDIS_VOLUMIC_POWER_NONE 0 /* <=> No volumic power */
+#define SDIS_VOLUMIC_POWER_NONE DBL_MAX /* <=> No volumic power */
 #define SDIS_FLUX_NONE DBL_MAX /* <=> No flux */
 #define SDIS_PRIMITIVE_NONE SIZE_MAX /* Invalid primitive */
 
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -373,6 +373,44 @@ solve
   }
 }
 
+static void
+check_null_power_term_with_green
+  (struct sdis_scene* scn,
+   struct ssp_rng* rng,
+   struct solid* solid)
+{
+  struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+  struct sdis_mc T = SDIS_MC_NULL;
+  struct sdis_estimator* estimator = NULL;
+  struct sdis_green_function* green = NULL;
+  double x = 0;
+  double ref = 0;
+  ASSERT(scn && rng && solid);
+
+  solve_args.position[0] = 0.5;
+  solve_args.position[1] = 0.5;
+  solve_args.position[2] = 0;
+  solve_args.nrealisations = N;
+  solve_args.time_range[0] =
+  solve_args.time_range[1] = INF;
+
+  solid->vpower = 0;
+  OK(sdis_solve_probe_green_function(scn, &solve_args, &green));
+
+  solid->vpower = P0;
+  OK(sdis_green_function_solve(green, &estimator));
+  OK(sdis_estimator_get_temperature(estimator, &T));
+
+  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);
+  CHK(eq_eps(ref, T.E, 3*T.SE));
+
+  OK(sdis_estimator_ref_put(estimator));
+  OK(sdis_green_function_ref_put(green));
+}
+
 /*******************************************************************************
  * Test
  ******************************************************************************/
@@ -485,12 +523,15 @@ main(int argc, char** argv)
   OK(sdis_interface_ref_put(interf_T0));
 
   /* Solve */
-  OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng));
+  OK(ssp_rng_create(NULL, SSP_RNG_MT19937_64, &rng));
   printf(">> Box scene\n");
   solve(box_scn, rng, solid_props);
   printf(">> Square scene\n");
   solve(square_scn, rng, solid_props);
 
+  /* Check green registration with a null power term */
+  check_null_power_term_with_green(box_scn, rng, solid_props);
+
   OK(sdis_scene_ref_put(box_scn));
   OK(sdis_scene_ref_put(square_scn));
   OK(sdis_device_ref_put(dev));