stardis-solver

commit 27c91d56f40bfdda3e6af27cd8c0d208abe3afc9
parent af831a26c8b2a3f72c082aab2cb4e726348fcf39
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  1 Mar 2019 10:14:28 +0100

Fix the volumic power corrective term.

Adjust how the RAY_RANGE_SCALE_MAX factor was handle : scale the origin
"delta_solid" by it only to find the limit angle, not in the term to
correct.

Diffstat:
M
src/sdis_heat_path_conductive_Xd.h
 | 
19
++++++++++++-------
M
src/sdis_misc.h
 | 
2
+-

2 files changed, 13 insertions(+), 8 deletions(-)
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -142,7 +142,8 @@ XD(conductive_path)
 
       /* Add the volumic power density to the measured temperature */
       if(power != SDIS_VOLUMIC_POWER_NONE) {
-        const double delta_s_in_meter = delta_solid * RAY_RANGE_MAX_SCALE * fp_to_meter;
+        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;
@@ -163,9 +164,16 @@ XD(conductive_path)
         /* The regular power term at wall */
         tmp = h_in_meter * h_in_meter / (2.0 * lambda);
 
-        /* Add the power corrective term */
-        if(h < delta_solid) {
-          const double sin_a = h / delta_solid;
+        /* 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);
@@ -174,9 +182,6 @@ XD(conductive_path)
           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
-
-        } else if(h == delta_solid) {
-          tmp += -(delta_s_in_meter * delta_s_in_meter)/(2.0*DIM*lambda);
         }
         power_factor = tmp;
         T->value += power * power_factor;
diff --git a/src/sdis_misc.h b/src/sdis_misc.h
@@ -89,7 +89,7 @@ sample_time(struct ssp_rng* rng, const double time_range[2])
 
 static INLINE res_T
 register_heat_vertex
-  (struct sdis_heat_path* path, 
+  (struct sdis_heat_path* path,
    const struct sdis_rwalk_vertex* vtx,
    const double weight,
    const enum sdis_heat_vertex_type type)