stardis-solver

commit 3ee7be850308d343ea9d8028e63f12e40bd51e6f
parent 5db76da4ea9d5580842ec9d4cd407176dc915c50
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  5 Nov 2021 16:45:39 +0100

Fix issues in the picardN algorithm

Diffstat:
M
src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
 | 
15
++++++++-------
M
src/sdis_realisation_Xd.h
 | 
22
++++++++++++++++++----

2 files changed, 26 insertions(+), 11 deletions(-)
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -40,7 +40,7 @@ restart_heat_path
   size_t nbreaks = 0;
   res_T res = RES_OK;
 
-  if(path) goto exit;
+  if(!path) goto exit;
   ASSERT(vtx);
 
   nbreaks = darray_size_t_size_get(&path->breaks);
@@ -78,6 +78,7 @@ XD(sample_path)
 
   /* Clean-up the output variable */
   *T = XD(TEMPERATURE_NULL);
+  T->func = XD(boundary_path);
 
   /* Init the random walk */
   rwalk.vtx = rwalk_from->vtx;
@@ -286,7 +287,7 @@ XD(solid_fluid_boundary_picardN_path)
     if(ctx->heat_path) hvtx_s = *heat_path_get_last_vertex(ctx->heat_path);
 
     h_radi_min = 4.0 * BOLTZMANN_CONSTANT * Tmin3 * epsilon;
-    p_radi_min = h_radi_min / h_radi_hat;
+    p_radi_min = h_radi_min / h_hat;
 
     /* Switch in radiative path */
     if(r < p_conv + p_cond + p_radi_min) { *rwalk = rwalk_s; *T = T_s; break; }
@@ -311,16 +312,16 @@ XD(solid_fluid_boundary_picardN_path)
     } (void)0
 
     #define CHECK_PMIN_PMAX {                                                  \
-      p_radi_min = h_radi_min*epsilon / h_radi_hat;                            \
-      p_radi_max = h_radi_max*epsilon / h_radi_hat;                            \
+      p_radi_min = h_radi_min*epsilon / h_hat;                                 \
+      p_radi_max = h_radi_max*epsilon / h_hat;                                 \
       if(r < p_conv + p_cond + p_radi_min) { SWITCH_IN_RADIATIVE; break; }     \
       if(r > p_conv + p_cond + p_radi_max) { NULL_COLLISION; continue; }       \
     } (void)0
 
     /* Sample a 1st heat path at the end of the radiative path */
     COMPUTE_TEMPERATURE(T0, &rwalk_s);
-    h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 * 3*Tmin2*T0);
-    h_radi_max = BOLTZMANN_CONSTANT*(That3 * 3*That2*T0);
+    h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + 3*Tmin2*T0);
+    h_radi_max = BOLTZMANN_CONSTANT*(That3 + 3*That2*T0);
     CHECK_PMIN_PMAX;
 
     /* Sample a 2nd heat path at the end of the radiative path */
@@ -350,7 +351,7 @@ XD(solid_fluid_boundary_picardN_path)
     /* Sample a 3rd heat path at the current position onto the interface */
     COMPUTE_TEMPERATURE(T5, rwalk);
     h_radi = BOLTZMANN_CONSTANT*(T3*T4*T5 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2);
-    p_radi = h_radi * epsilon / h_radi_hat;
+    p_radi = h_radi * epsilon / h_hat;
 
     /* Switch in radiative path */
     if(r < p_cond + p_conv + p_radi) { SWITCH_IN_RADIATIVE; break; }
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -181,8 +181,12 @@ XD(probe_realisation)
 
   ctx.green_path = green_path;
   ctx.heat_path = heat_path;
-  ctx.That2 = scn->tmax * scn->tmax;
-  ctx.That3 = scn->tmax * ctx.That2;
+  ctx.Tmin  = scn->tmin;
+  ctx.Tmin2 = ctx.Tmin * ctx.Tmin;
+  ctx.Tmin3 = ctx.Tmin * ctx.Tmin2;
+  ctx.That  = scn->tmax;
+  ctx.That2 = ctx.That * ctx.That;
+  ctx.That3 = ctx.That * ctx.That2;
 
   res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
   if(res != RES_OK) goto error;
@@ -256,8 +260,12 @@ XD(boundary_realisation)
 
   ctx.green_path = green_path;
   ctx.heat_path = heat_path;
-  ctx.That2 = scn->tmax * scn->tmax;
-  ctx.That3 = scn->tmax * ctx.That2;
+  ctx.Tmin  = scn->tmin;
+  ctx.Tmin2 = ctx.Tmin * ctx.Tmin;
+  ctx.Tmin3 = ctx.Tmin * ctx.Tmin2;
+  ctx.That  = scn->tmax;
+  ctx.That2 = ctx.That * ctx.That;
+  ctx.That3 = ctx.That * ctx.That2;
 
   res = XD(compute_temperature)(scn, &ctx, &rwalk, rng, &T);
   if(res != RES_OK) goto error;
@@ -296,6 +304,9 @@ XD(boundary_flux_realisation)
 #endif
   double P[SDIS_XD_DIMENSION];
   float N[SDIS_XD_DIMENSION];
+  const double Tmin = scn->tmin;
+  const double Tmin2 = Tmin * Tmin;
+  const double Tmin3 = Tmin * Tmin2;
   const double That = scn->tmax;
   const double That2 = That * That;
   const double That3 = That * That2;
@@ -333,6 +344,9 @@ XD(boundary_flux_realisation)
     rwalk.mdm = (Mdm);                                                         \
     rwalk.hit.prim = prim;                                                     \
     SET_PARAM(rwalk.hit, st);                                                  \
+    ctx.Tmin  = Tmin;                                                          \
+    ctx.Tmin3 = Tmin3;                                                         \
+    ctx.That  = That;                                                          \
     ctx.That2 = That2;                                                         \
     ctx.That3 = That3;                                                         \
     dX(set)(rwalk.vtx.P, P);                                                   \