Add spherical test and fix broken variable spreading in spherical

include/world_builder/features/oceanic_plate_models/temperature/half_space_model.h

@@ -89,7 +89,7 @@ namespace WorldBuilder
             Objects::Surface max_depth_surface;
             double top_temperature;
             double bottom_temperature;
-            std::pair<std::vector<double>,std::vector<double>> spreading_velocities;
+            std::vector<std::vector<double>> spreading_velocities;
             std::vector<std::vector<Point<2> > > mid_oceanic_ridges;
             std::vector<std::vector<double>> spreading_velocities_at_each_ridge_point;
             Operations operation;

include/world_builder/types/value_at_points.h

@@ -40,7 +40,7 @@ namespace WorldBuilder
          * A constructor
          */
         ValueAtPoints(const double default_value,
-                      const double max_values_in_array,
+                      size_t max_values_in_array,
                       std::vector<Point<2>> default_points_ = std::vector<Point<2>>());
 
         /**

include/world_builder/utilities.h

@@ -426,6 +426,24 @@ namespace WorldBuilder
     */
     std::string
     read_and_distribute_file_content(const std::string &filename);
+
+    /**
+     * Calculate the distance of a point from a mid oceanic ridge, and also calulate
+     * the spreading velocity of the ridge at this point.
+     * TODO: make the spreading velocity spatially/temporally variable
+     *
+     * @param mid_oceanic_ridges The coordinates of the mid oceanic ridges
+     * @param spreading_velocity The spreading rate of the mid oceanic ridges
+     * @param coordinate_system The coordinate system
+     * @param position_in_natural_coordinates_at_min_depth the current position in natural_coordinates
+     * @return The content of the file.
+    */
+    std::pair<double, double>
+    calculate_ridge_distance_and_spreading(std::vector<std::vector<Point<2>>> mid_oceanic_ridges,
+                                           std::vector<std::vector<double>> mid_oceanic_spreading_velocities,
+                                           const std::unique_ptr<WorldBuilder::CoordinateSystems::Interface> &coordinate_system,
+                                           const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth);
+
   } // namespace Utilities
 } // namespace WorldBuilder
 

source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc

@@ -83,7 +83,7 @@ namespace WorldBuilder
                             "in degree Kelvin for this feature. If the model has an adiabatic gradient"
                             "this should be the mantle potential temperature, and T = Tad + Thalf. ");
 
-          prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., std::numeric_limits<double>::max()))),
+          prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::Array(Types::Double(0.01)))),
                             "The spreading velocity of the plate in meter per year. "
                             "This is the velocity with which one side moves away from the ridge.");
 
@@ -105,27 +105,32 @@ namespace WorldBuilder
           operation = string_operations_to_enum(prm.get<std::string>("operation"));
           top_temperature = prm.get<double>("top temperature");
           bottom_temperature = prm.get<double>("bottom temperature");
-          spreading_velocities = prm.get_value_at_array("spreading velocity");
+          spreading_velocities = prm.get_vector<std::vector<double>>("spreading velocity");
 
           mid_oceanic_ridges = prm.get_vector<std::vector<Point<2>>>("ridge coordinates");
           const double dtr = prm.coordinate_system->natural_coordinate_system() == spherical ? Consts::PI / 180.0 : 1.0;
-
-          unsigned int index_x = 0;
-          unsigned int index_y = 0;
-          unsigned int test_ind = 0;
-          for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
-            {
-              std::vector<double> spreading_rates_for_ridge;
-              for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
-                {
-                  if (spreading_velocities.second.size() <= 1)
-                    spreading_rates_for_ridge.push_back(spreading_velocities.first[0]);
-                  else
-                    spreading_rates_for_ridge.push_back(spreading_velocities.second[test_ind]);
-                  test_ind += 1;
-                }
-              spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge);
-            }
+          for (auto &ridge_coordinates : mid_oceanic_ridges)
+            for (auto &ridge_coordinate : ridge_coordinates)
+              {
+                ridge_coordinate *= dtr;
+              }
+
+        //   unsigned int index_x = 0;
+        //   unsigned int index_y = 0;
+        //   unsigned int ridge_point_index = 0;
+        //   for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
+        //     {
+        //       std::vector<double> spreading_rates_for_ridge;
+        //       for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
+        //         {
+        //           if (spreading_velocities.second.size() <= 1)
+        //             spreading_rates_for_ridge.push_back(spreading_velocities.first[0] / 31557600); // m/s
+        //           else
+        //             spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index] / 31557600); // m/s
+        //           ridge_point_index += 1;
+        //         }
+        //       spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge);
+        //     }
         }
 
         double
@@ -156,7 +161,10 @@ namespace WorldBuilder
                                                             this->world->specific_heat) * depth);
                     }
 
-                  double distance_ridge = std::numeric_limits<double>::max();
+                  std::pair<double, double> ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges,
+                                                               spreading_velocities,
+                                                               world->parameters.coordinate_system,
+                                                               position_in_natural_coordinates_at_min_depth);
 
                   const CoordinateSystem coordinate_system = world->parameters.coordinate_system->natural_coordinate_system();
 
@@ -272,12 +280,12 @@ namespace WorldBuilder
                            << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
                            << ", top_temperature = " << top_temperature
                            << ", max_depth = " << max_depth
-                           << ", spreading_velocity = " << spreading_velocity
+                           << ", spreading_velocity = " << ridge_parameters.first
                            << ", thermal_diffusivity = " << thermal_diffusivity
                            << ", age = " << age << '.');
                   WBAssert(std::isfinite(temperature), "Temperature inside half-space cooling model is not a finite: " << temperature                           << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
                            << ", top_temperature = " << top_temperature
-                           << ", spreading_velocity = " << spreading_velocity
+                           << ", spreading_velocity = " << ridge_parameters.first
                            << ", thermal_diffusivity = " << thermal_diffusivity
                            << ", age = " << age << '.');
 

source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc

@@ -79,7 +79,7 @@ namespace WorldBuilder
           prm.declare_entry("bottom temperature", Types::Double(-1),
                             "The temperature in degree Kelvin which this feature should have");
 
-          prm.declare_entry("spreading velocity", Types::Double(-1),
+          prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::Array(Types::Double(0.01)))),
                             "The spreading velocity of the plate in meter per year. "
                             "This is the velocity with which one side moves away from the ridge.");
 

source/world_builder/types/value_at_points.cc

@@ -29,7 +29,7 @@ namespace WorldBuilder
   {
 
     ValueAtPoints::ValueAtPoints(const double default_value_,
-                                 const double max_values_in_array_,
+                                 size_t max_values_in_array_,
                                  std::vector<Point<2>> default_points_)
       :
       default_value(default_value_),
@@ -69,7 +69,7 @@ namespace WorldBuilder
         Pointer((base + "/type").c_str()).Set(declarations,"array");
         Pointer((base + "/additionalProperties").c_str()).Set(declarations,false);
         Pointer((base + "/minItems").c_str()).Set(declarations,1);
-        Pointer((base + "/maxItems").c_str()).Set(declarations,2);
+        Pointer((base + "/maxItems").c_str()).Set(declarations, max_values_in_array);
         Pointer((base + "/documentation").c_str()).Set(declarations,documentation.c_str());
 
         {

source/world_builder/utilities.cc

@@ -1219,6 +1219,143 @@ namespace WorldBuilder
 
     template std::array<double,2> convert_point_to_array<2>(const Point<2> &point_);
     template std::array<double,3> convert_point_to_array<3>(const Point<3> &point_);
+
+
+    std::pair<double, double>
+    calculate_ridge_distance_and_spreading(std::vector<std::vector<Point<2>>> mid_oceanic_ridges,
+                                           std::vector<std::vector<double>> mid_oceanic_spreading_velocities,
+                                           const std::unique_ptr<WorldBuilder::CoordinateSystems::Interface> &coordinate_system,
+                                           const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth)
+    {
+
+      double distance_ridge = std::numeric_limits<double>::max();
+      double spreading_velocity_at_ridge = 0;
+
+      // first find if the coordinate is on this side of a ridge
+      unsigned int relevant_ridge = 0;
+      const Point<2> check_point(position_in_natural_coordinates_at_min_depth.get_surface_coordinates(),
+                                 position_in_natural_coordinates_at_min_depth.get_coordinate_system());
+
+      Point<2> other_check_point = check_point;
+      if (check_point.get_coordinate_system() == CoordinateSystem::spherical)
+        {
+          other_check_point[0] += check_point[0] < 0 ? 2.0 * WorldBuilder::Consts::PI : -2.0 * WorldBuilder::Consts::PI;
+        }
+
+      // if there is only one ridge, there is no transform
+      if (mid_oceanic_ridges.size() > 1)
+        {
+          // There are more than one ridge, so there are transform faults
+          // Find the first which is on the same side
+          for (relevant_ridge = 0; relevant_ridge < mid_oceanic_ridges.size()-1; relevant_ridge++)
+            {
+              const Point<2> transform_point_0 = mid_oceanic_ridges[relevant_ridge+1][0];
+              const Point<2> transform_point_1 = mid_oceanic_ridges[relevant_ridge][mid_oceanic_ridges[relevant_ridge].size()-1];
+              const Point<2> reference_point   = mid_oceanic_ridges[relevant_ridge][0];
+
+              const bool reference_on_side_of_line = (transform_point_1[0] - transform_point_0[0])
+                                                     * (reference_point[1] - transform_point_0[1])
+                                                     - (transform_point_1[1] - transform_point_0[1])
+                                                     * (reference_point[0] - transform_point_0[0])
+                                                     < 0;
+              const bool checkpoint_on_side_of_line = (transform_point_1[0] - transform_point_0[0])
+                                                      * (check_point[1] - transform_point_0[1])
+                                                      - (transform_point_1[1] - transform_point_0[1])
+                                                      * (check_point[0] - transform_point_0[0])
+                                                      < 0;
+
+
+              if (reference_on_side_of_line == checkpoint_on_side_of_line)
+                {
+                  break;
+                }
+
+            }
+        }
+
+      for (unsigned int i_coordinate = 0; i_coordinate < mid_oceanic_ridges[relevant_ridge].size() - 1; i_coordinate++)
+        {
+          const Point<2> segment_point0 = mid_oceanic_ridges[relevant_ridge][i_coordinate];
+          const Point<2> segment_point1 = mid_oceanic_ridges[relevant_ridge][i_coordinate + 1];
+
+          const double spreading_velocity_point0 = mid_oceanic_spreading_velocities[relevant_ridge][i_coordinate];
+          const double spreading_velocity_point1 = mid_oceanic_spreading_velocities[relevant_ridge][i_coordinate + 1];
+
+          {
+            // based on http://geomalgorithms.com/a02-_lines.html
+            const Point<2> v = segment_point1 - segment_point0;
+            const Point<2> w1 = check_point - segment_point0;
+            const Point<2> w2 = other_check_point - segment_point0;
+
+            const double c1 = (w1[0] * v[0] + w1[1] * v[1]);
+            const double c = (v[0] * v[0] + v[1] * v[1]);
+            const double c2 = (w2[0] * v[0] + w2[1] * v[1]);
+
+
+            Point<2> Pb1(coordinate_system->natural_coordinate_system());
+            // This part is needed when we want to consider segments instead of lines
+            // If you want to have infinite lines, use only the else statement.
+
+            if (c1 <= 0)
+              {
+                Pb1=segment_point0;
+                spreading_velocity_at_ridge = spreading_velocity_point0;
+              }
+            else if (c <= c1)
+              {
+                Pb1=segment_point1;
+                spreading_velocity_at_ridge = spreading_velocity_point1;
+              }
+            else
+              {
+                Pb1=segment_point0 + (c1 / c) * v;
+                spreading_velocity_at_ridge = spreading_velocity_point0 + (spreading_velocity_point1 - spreading_velocity_point0) * (c1 / c) * 1;
+              }
+
+            Point<2> Pb2(coordinate_system->natural_coordinate_system());
+            if (c2 <= 0)
+              {
+                Pb2=segment_point0;
+                spreading_velocity_at_ridge = spreading_velocity_point0;
+              }
+            else if (c <= c2)
+              {
+                Pb2=segment_point1;
+                spreading_velocity_at_ridge = spreading_velocity_point1;
+              }
+            else
+              {
+                Pb2=segment_point0 + (c1 / c) * v;
+                spreading_velocity_at_ridge = spreading_velocity_point0 + (spreading_velocity_point1 - spreading_velocity_point0) * (c1 / c) * 1;
+              }
+
+            Point<3> compare_point1(coordinate_system->natural_coordinate_system());
+            Point<3> compare_point2(coordinate_system->natural_coordinate_system());
+
+            compare_point1[0] = coordinate_system->natural_coordinate_system() == cartesian ? Pb1[0] :  position_in_natural_coordinates_at_min_depth.get_depth_coordinate();
+            compare_point1[1] = coordinate_system->natural_coordinate_system() == cartesian ? Pb1[1] : Pb1[0];
+            compare_point1[2] = coordinate_system->natural_coordinate_system() == cartesian ? position_in_natural_coordinates_at_min_depth.get_depth_coordinate() : Pb1[1];
+
+            compare_point2[0] = coordinate_system->natural_coordinate_system() == cartesian ? Pb2[0] :  position_in_natural_coordinates_at_min_depth.get_depth_coordinate();
+            compare_point2[1] = coordinate_system->natural_coordinate_system() == cartesian ? Pb2[1] : Pb2[0];
+            compare_point2[2] = coordinate_system->natural_coordinate_system() == cartesian ? position_in_natural_coordinates_at_min_depth.get_depth_coordinate() : Pb2[1];
+
+            distance_ridge = std::min(distance_ridge,
+                                      coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),
+                                          position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
+                                          compare_point1));
+
+            distance_ridge = std::min(distance_ridge,
+                                      coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),
+                                          position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
+                                          compare_point2));
+          }
+        }
+      std::pair<double, double> result;
+      result.first = spreading_velocity_at_ridge;
+      result.second = distance_ridge;
+      return result;
+    }
   } // namespace Utilities
 } // namespace WorldBuilder
 

tests/gwb-dat/spherical_variable_spreading.dat

@@ -0,0 +1,15 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 3
+# compositions = 0
+# convert spherical = true
+# R long lat d T
+6321e3 22.5 16 50e3
+6321e3 22.5 18 50e3
+6321e3 22.5 20 50e3
+6321e3 22.5 22 50e3
+6321e3 22.5 23 50e3
+6321e3 22.5 25 50e3
+6321e3 22.5 27 50e3
+6321e3 22.5 29 50e3

tests/gwb-dat/spherical_variable_spreading.wb

@@ -0,0 +1,17 @@
+{
+  "version": "0.6",
+  "coordinate system":{"model":"spherical", "depth method":"starting point"},
+  "gravity model":{"model":"uniform", "magnitude":10},
+  "surface temperature":273, "force surface temperature":true,
+  "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, 
+  "specific heat":1250, "thermal diffusivity":1.0e-6,
+  "features":
+  [
+    { "model":"oceanic plate", "name":"Subducting", "max depth":250e3,"min depth":0,
+      "coordinates" :[[15, 15],[15, 30],[30, 30],[30, 15]],
+      "temperature models":[
+        {"model":"half space model", "min depth":0,  "max depth":100e3, "spreading velocity":[ [0,[[0.005, 0.0001, 0.005]]], [0, [[0.001, 0.00005]]] ],
+        "top temperature":273,
+         "ridge coordinates": [[[20,10], [21, 16], [20, 22.5]], [[25, 22.5], [25, 31]]]}]}      
+  ]
+}

tests/gwb-dat/spherical_variable_spreading/screen-output.log

@@ -0,0 +1,10 @@
+# x y z d g T 
+6321e3 22.5 16 50e3 1315.37 
+6321e3 22.5 18 50e3 1258.52 
+6321e3 22.5 20 50e3 1208.46 
+6321e3 22.5 22 50e3 1166.44 
+6321e3 22.5 23 50e3 696.071 
+6321e3 22.5 25 50e3 647.504 
+6321e3 22.5 27 50e3 588.45 
+6321e3 22.5 29 50e3 510.556 
+