Merge pull request #677 from MFraters/fix_nan_add_asserts

Tracing and fixing NaN in distance_point_from_curved_planes
GeodynamicWorldBuilder · Feb 27, 2024 · 601b6d2 · 601b6d2
2 parents b6ae1ad + 9efa94f
commit 601b6d2
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 1 deletion.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -33,7 +33,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - Fixed an issue where the ridge feature in spherical geometries for both the half space cooling and plate cooling models gave a discontinuous spreading center when crossing longitudes at intervals of 90 degrees. \[Daniel Douglas; 2024-01-22; [#520](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/520),[#518](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/518)\]
 - In some cases the bezier curve closest_point_on_curve_segment function would include a half circle around the end point(s) as part of a slab/fault. \[Menno Fraters, reported by Daniel Douglas;2023-12-07; [#522](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/522) and [#523](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/523)\]
 - In some cases the fault feature would not recognize points as inside the fault if they were exactly on the fault line. This is fixed now. \[Rene Gassmoeller; 2024-02-16; [#640](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/640)\]
-- The input name for the annulus was misspelled as `anullus`. This is now fixed. \[Menno Fraters and Juliane Dannberg; 2024-02-22; [#539](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/539) and [#666](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/666)]
+- The input name for the annulus was misspelled as `anullus`. This is now fixed. \[Menno Fraters and Juliane Dannberg; 2024-02-22; [#539](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/539) and [#666](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/666)\]
+- Fixed bug in distance_point_from_curved_planes function. It was not setting the depth_reference_surface variable, leaving it NaN, when it the check point is exactly on a trench/fault point. \[Menno Fraters; 2024-02-26; [#677](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/677)\]
 
 ## [0.5.0]
 ### Added

diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc
@@ -17,6 +17,7 @@
    along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
+#include "world_builder/assert.h"
 #include <algorithm>
 #include <cmath>
 #include <fstream>
@@ -586,6 +587,7 @@ namespace WorldBuilder
                   return_values.section = i_section_min_distance;
                   return_values.segment = 0;
                   return_values.average_angle = total_average_angle;
+                  return_values.depth_reference_surface = 0.0;
                   return_values.closest_trench_point = closest_point_on_line_cartesian;
                   return return_values;
                 }
@@ -837,6 +839,10 @@ namespace WorldBuilder
                           new_distance = side_of_line * (check_point_2d - Pb).norm();
                           new_along_plane_distance = (begin_segment - Pb).norm();
                           new_depth_reference_surface = start_radius - Pb[1];
+
+                          WBAssert(!std::isnan(new_depth_reference_surface),
+                                   "new_depth_reference_surface is not a number: " << new_depth_reference_surface << ". "
+                                   << "start_radius = " << start_radius << ",Pb[1] = " << Pb[1] << ".");
                         }
                     }
                 }
@@ -971,6 +977,11 @@ namespace WorldBuilder
                       new_along_plane_distance = (radius_angle_circle * check_point_angle - radius_angle_circle * interpolated_angle_top) * (difference_in_angle_along_segment < 0 ? 1 : -1);
                       // compute the new depth by rotating the begin point to the check point location.
                       new_depth_reference_surface = start_radius-(sin(check_point_angle + interpolated_angle_top) * BSPC[0] + cos(check_point_angle + interpolated_angle_top) * BSPC[1] + center_circle[1]);
+
+                      WBAssert(!std::isnan(new_depth_reference_surface),
+                               "new_depth_reference_surface is not a number: " << new_depth_reference_surface << ". "
+                               << "start_radius = " << start_radius << ", check_point_angle = " << check_point_angle << ", interpolated_angle_top = " << interpolated_angle_top
+                               << ", BSPC[0] = " << BSPC[0] << ".");
                     }
 
                 }
@@ -1011,6 +1022,8 @@ namespace WorldBuilder
             }
         }
 
+      WBAssert(!std::isnan(depth_reference_surface), "depth_reference_surface is not a number: " << depth_reference_surface << ".");
+
       PointDistanceFromCurvedPlanes return_values(natural_coordinate.get_coordinate_system());
       return_values.distance_from_plane = distance;
       return_values.distance_along_plane = along_plane_distance;