stardis-solver

commit 0f451ff9ec73e3e4578023e45b08f9cb709e1de1
parent 1ab9726c0f55a151756d1b0512dce35dc1d32e24
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 24 Apr 2019 08:58:12 +0200

Remove the volumetric power corrective term

Actually, this term should be still handled to correct the bias introduced
by the random walk steps that progress near the boundaries. However, the new
conductive random walk algorithm that uses several directions to
estimate the nearest distance to the boundary drastically reduces the
bias that actually is now really hard to observe in practice.

Diffstat:
M
src/sdis_heat_path_conductive_Xd.h
 | 
71
++++++++++++++++++++++-------------------------------------------------

1 file changed, 22 insertions(+), 49 deletions(-)
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -59,10 +59,15 @@ XD(sample_next_step)
   ssp_ran_circle_uniform_float(rng, dirs[0], NULL);
 
   /* Compute in dirs[2] a direction orthogonal to dirs[0] */
+#if 0
+  dirs[2][0] =  dirs[0][0];
+  dirs[2][1] =  dirs[0][1];
+#else
   dirs[2][0] = -dirs[0][1];
   dirs[2][1] =  dirs[0][0];
   ASSERT(f2_is_normalized(dirs[2]));
   ASSERT(eq_epsf(f2_dot(dirs[0], dirs[2]), 0, 1.e-6f));
+#endif
 
   /* Negate the orthornormal frame */
   f2_minus(dirs[1], dirs[0]);
@@ -240,55 +245,9 @@ XD(conductive_path)
 
     /* Add the volumic power density to the measured temperature */
     if(power != SDIS_VOLUMIC_POWER_NONE) {
-      if((S3D_HIT_NONE(&hit0) && S3D_HIT_NONE(&hit1))) { /* Hit nothing */
-        const double delta_in_meter = delta * fp_to_meter;
-        power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
-        T->value += power * power_factor;
-      } else {
-        const double delta_s_adjusted = delta_solid * RAY_RANGE_MAX_SCALE;
-        const double delta_s_in_meter = delta_solid * fp_to_meter;
-        double h;
-        double h_in_meter;
-        double cos_U_N;
-        float N[DIM];
-
-        if(delta == hit0.distance) {
-          fX(normalize)(N, hit0.normal);
-          cos_U_N = fX(dot)(dir0, N);
-        } else {
-          ASSERT(delta == hit1.distance);
-          fX(normalize)(N, hit1.normal);
-          cos_U_N = fX(dot)(dir1, N);
-        }
-
-        h = delta * fabs(cos_U_N);
-        h_in_meter = h * fp_to_meter;
-
-        /* The regular power term at wall */
-        tmp = h_in_meter * h_in_meter / (2.0 * lambda);
-
-        /* Add the power corrective term. Be careful to use the adjusted
-         * delta_solid to correctly handle the RAY_RANGE_MAX_SCALE factor in
-         * the computation of the limit angle. But keep going with the
-         * unmodified delta_solid in the corrective term since it was the one
-         * that was "wrongly" used in the previous step and that must be
-         * corrected. */
-        if(h == delta_s_adjusted) {
-          tmp += -(delta_s_in_meter * delta_s_in_meter)/(2.0*DIM*lambda);
-        } else if(h < delta_s_adjusted) {
-          const double sin_a = h / delta_s_adjusted;
-#if DIM==2
-          /* tmp1 = sin(2a) / (PI - 2*a) */
-          const double tmp1 = sin_a * sqrt(1 - sin_a*sin_a)/acos(sin_a);
-          tmp += -(delta_s_in_meter * delta_s_in_meter)/(4.0*lambda) * tmp1;
-#else
-          const double tmp1 = (sin_a*sin_a*sin_a - sin_a)/ (1-sin_a);
-          tmp += (delta_s_in_meter * delta_s_in_meter)/(6.0*lambda) * tmp1;
-#endif
-        }
-        power_factor = tmp;
-        T->value += power * power_factor;
-      }
+      const double delta_in_meter = delta * fp_to_meter;
+      power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
+      T->value += power * power_factor;
     }
 
     /* Register the power term for the green function. Delay its registration
@@ -387,6 +346,20 @@ XD(conductive_path)
 
     ++istep;
 
+#if 0
+   {
+#define Tf 100.0
+#define Power 10000.0
+#define H 50.0
+#define LAMBDA 100.0
+#define DELTA 0.4/*(1.0/2.0)*/
+
+      T->value += -Power / (2*LAMBDA) * rwalk->vtx.P[1]*rwalk->vtx.P[1] + Tf + Power/(2*H) + Power/(8*LAMBDA);
+      T->done = 1;
+      break;
+    }
+#endif
+
   /* Keep going while the solid random walk does not hit an interface */
   } while(SXD_HIT_NONE(&rwalk->hit));