stardis-solver

commit f8ac59923c13813a8d043fad679390e6a9d07548
parent fa23f7176bd7cca747ee3ff83ec26d1fe7c1e249
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 23 Nov 2021 16:31:21 +0100

Check the picard order with a flux > 0

Picard N cannot be used with a flux > 0

Diffstat:
M
src/sdis_heat_path_boundary_Xd_fixed_flux.h
 | 
16
+++++++++++++++-
M
src/test_sdis_flux.c
 | 
10
++++++++++

2 files changed, 25 insertions(+), 1 deletion(-)
diff --git a/src/sdis_heat_path_boundary_Xd_fixed_flux.h b/src/sdis_heat_path_boundary_Xd_fixed_flux.h
@@ -51,12 +51,13 @@ XD(solid_boundary_with_flux_path)
   struct sdis_interface* interf = NULL;
   struct sdis_medium* solid = NULL;
 
-  /* Miscellaneous terms */
+  /* Miscellaneous variables */
   double lambda; /* Solid conductivity */
   double delta_boundary; /* Orthogonal reinjection dst at the boundary */
   double delta; /* Orthogonal fitted reinjection dst at the boundary */
   double delta_m; /* Delta in meters */
   double flux_term;
+  size_t picard_order;
   enum sdis_side solid_side = SDIS_SIDE_NULL__;
   res_T res = RES_OK;
 
@@ -64,6 +65,19 @@ XD(solid_boundary_with_flux_path)
   ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
   (void)ctx;
 
+  /* Currently, the flux can be correctly taken into account only when the
+   * radiative temperature is linearized, i.e. when the picard order is equal
+   * to 1 */
+  picard_order = get_picard_order(ctx);
+  if(picard_order > 1 && phi > 0) {
+    log_err(scn->dev,
+      "%s: invalid flux '%g' W/m^2. Could not manage a flux != 0 when the "
+      "picard order is not equal to 1; Picard order is currently set to %lu.\n",
+      FUNC_NAME, phi, (unsigned long)picard_order);
+    res = RES_BAD_ARG;
+    goto error;
+  }
+
   /* Retrieve the solid split by the interface */
   interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
   solid = interface_get_medium(interf, frag->side);
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -153,6 +153,7 @@ solve
   ASSERT(scn && rng && interf);
 
   OK(sdis_scene_get_dimension(scn, &dim));
+
   FOR_EACH(isimul, 0, nsimuls) {
     int steady = (isimul % 2) == 0;
 
@@ -335,6 +336,15 @@ solve
 
     printf("\n\n");
   }
+
+  /* Picard N is not supported with a flux != 0 */
+  solve_args.position[0] = 0.1;
+  solve_args.position[1] = 0.1;
+  solve_args.position[2] = dim == SDIS_SCENE_2D ? 0 : 0.1;
+  solve_args.time_range[0] = INF;
+  solve_args.time_range[1] = INF;
+  solve_args.picard_order = 2;
+  BA(sdis_solve_probe(scn, &solve_args, &estimator));
 }
 
 /*******************************************************************************