along interface velocity for subducting slab

include/world_builder/features/subducting_plate_models/velocity/along_surface.h

@@ -0,0 +1,95 @@
+/*
+  Copyright (C) 2018-2026 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef WORLD_BUILDER_FEATURES_SUBDUCTING_PLATE_MODELS_VELOCITY_ALONG_SURFACE_H
+#define WORLD_BUILDER_FEATURES_SUBDUCTING_PLATE_MODELS_VELOCITY_ALONG_SURFACE_H
+
+
+#include "world_builder/features/subducting_plate_models/velocity/interface.h"
+#include "world_builder/features/feature_utilities.h"
+
+
+namespace WorldBuilder
+{
+
+  namespace Features
+  {
+    namespace SubductingPlateModels
+    {
+      namespace Velocity
+      {
+        /**
+         * This class represents a subducting plate and can implement submodules
+         * for velocity and composition. These submodules determine what
+         * the returned velocity or composition of the velocity and composition
+         * functions of this class will be.
+         */
+        class AlongSurface final: public Interface
+        {
+          public:
+            /**
+             * constructor
+             */
+            AlongSurface(WorldBuilder::World *world);
+
+            /**
+             * Destructor
+             */
+            ~AlongSurface() override final;
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            static
+            void declare_entries(Parameters &prm, const std::string &parent_name = "");
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            void parse_entries(Parameters &prm) override final;
+
+
+            /**
+             * Returns a velocity based on the given position, depth in the model,
+             * gravity and current velocity.
+             */
+            std::array<double,3> get_velocity(const Point<3> &position,
+                                              const double depth,
+                                              const double gravity,
+                                              std::array<double,3> velocity,
+                                              const double feature_min_depth,
+                                              const double feature_max_depth,
+                                              const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_planes,
+                                              const AdditionalParameters &additional_parameters) const override final;
+
+
+          private:
+            // uniform raw velocity submodule parameters
+            double min_depth;
+            double max_depth;
+            double velocity_magnitude;
+            Operations operation;
+
+        };
+      } // namespace Velocity
+    } // namespace SubductingPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+
+#endif

include/world_builder/utilities.h

@@ -308,6 +308,7 @@ namespace WorldBuilder
         fraction_of_segment(NaN::DSNAN),
         section(NaN::ISNAN),
         segment(NaN::ISNAN),
+        angle(NaN::DSNAN),
         average_angle(NaN::DSNAN),
         depth_reference_surface(NaN::DSNAN),
         closest_trench_point(Point<3>(coordinate_system))
@@ -346,6 +347,13 @@ namespace WorldBuilder
        */
       size_t segment;
 
+      // todo_velo
+      /**
+       * The dip angle of the plane at the location where the
+       * point is projected onto the plane.
+       */
+      double angle;
+
       /**
        * The average dip angle of the plane at the location where the
        * point is projected onto the plane.

source/world_builder/features/subducting_plate_models/velocity/along_surface.cc

@@ -0,0 +1,122 @@
+/*
+  Copyright (C) 2018-2026 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "world_builder/features/subducting_plate_models/velocity/along_surface.h"
+
+
+#include "world_builder/nan.h"
+#include "world_builder/types/array.h"
+#include "world_builder/types/double.h"
+#include "world_builder/types/object.h"
+#include "world_builder/utilities.h"
+
+
+namespace WorldBuilder
+{
+
+  using namespace Utilities;
+
+  namespace Features
+  {
+    namespace SubductingPlateModels
+    {
+      namespace Velocity
+      {
+        AlongSurface::AlongSurface(WorldBuilder::World *world_)
+          :
+          min_depth(NaN::DSNAN),
+          max_depth(NaN::DSNAN),
+          velocity_magnitude(NaN::DSNAN),
+          operation(Operations::REPLACE)
+        {
+          this->world = world_;
+          this->name = "along surface";
+        }
+
+        AlongSurface::~AlongSurface()
+          = default;
+
+        void
+        AlongSurface::declare_entries(Parameters &prm, const std::string & /*unused*/)
+        {
+          // Document plugin and require entries if needed.
+          // Add `velocity magnitude` and to the required parameters.
+          prm.declare_entry("", Types::Object({"velocity magnitude"}),
+                            "Uniform velocity model. Set the velocity to a constant value.");
+
+          // Declare entries of this plugin
+          prm.declare_entry("min distance slab top", Types::Double(0),
+                            "todo The depth in meters from which the composition of this feature is present.");
+
+          prm.declare_entry("max distance slab top", Types::Double(std::numeric_limits<double>::max()),
+                            "todo The depth in meters to which the composition of this feature is present.");
+
+          prm.declare_entry("velocity magnitude", Types::Double(0),
+                            "The velocity in meter per year");
+
+        }
+
+        void
+        AlongSurface::parse_entries(Parameters &prm)
+        {
+
+          min_depth = prm.get<double>("min distance slab top");
+          max_depth = prm.get<double>("max distance slab top");
+          operation = string_operations_to_enum(prm.get<std::string>("operation"));
+          velocity_magnitude = prm.get<double>("velocity magnitude");
+        }
+
+
+        std::array<double,3>
+        AlongSurface::get_velocity(const Point<3> & /*position_in_cartesian_coordinates*/,
+                                   const double  /*depth*/,
+                                   const double  /*gravity*/,
+                                   std::array<double,3> velocity_,
+                                   const double  /*feature_min_depth*/,
+                                   const double  /*feature_max_depth*/,
+                                   const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_plane,
+                                   const AdditionalParameters & /*additional_parameters*/) const
+        {
+
+          if (distance_from_plane.distance_from_plane <= max_depth && distance_from_plane.distance_from_plane >= min_depth)
+            {
+              const double angle = distance_from_plane.angle;
+
+              const double sign_vx = -1.0; // place holder for modify vx direction
+              const double vx = sign_vx * velocity_magnitude * std::cos(angle);
+              const double vy = 0.0;
+              const double vz = -velocity_magnitude * std::sin(angle);
+
+              return {{
+                  apply_operation(operation,velocity_[0],vx),
+                  apply_operation(operation,velocity_[1],vy),
+                  apply_operation(operation,velocity_[2],vz)
+                }
+              };
+            }
+
+          return velocity_;
+        }
+
+        WB_REGISTER_FEATURE_SUBDUCTING_PLATE_VELOCITY_MODEL(AlongSurface, along surface)
+      } // namespace Velocity
+    } // namespace SubductingPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+

source/world_builder/utilities.cc

@@ -444,8 +444,10 @@ namespace WorldBuilder
       // The 'vertical' fraction, indicates how far in the current segment the
       // point on the line is.
       double segment_fraction = 0.0;
+      double current_angle = 0.0;
       double total_average_angle = 0.0;
       double depth_reference_surface = 0.0;
+      double new_check_point_angle = 0.0;
 
       const DepthMethod depth_method = coordinate_system->depth_method();
 
@@ -856,6 +858,7 @@ namespace WorldBuilder
                           new_distance = std::numeric_limits<double>::infinity();
                           new_along_plane_distance = std::numeric_limits<double>::infinity();
                           new_depth_reference_surface = std::numeric_limits<double>::infinity();
+                          // new_check_point_angle = std::numeric_limits<double>::infinity();
                         }
                       else
                         {
@@ -867,6 +870,7 @@ 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];
+                          new_check_point_angle = interpolated_angle_top;
 
                           WBAssert(!std::isnan(new_depth_reference_surface),
                                    "new_depth_reference_surface is not a number: " << new_depth_reference_surface << ". "
@@ -1005,6 +1009,7 @@ 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]);
+                      new_check_point_angle = check_point_angle;
 
                       WBAssert(!std::isnan(new_depth_reference_surface),
                                "new_depth_reference_surface is not a number: " << new_depth_reference_surface << ". "
@@ -1033,8 +1038,9 @@ namespace WorldBuilder
                   section_fraction = fraction_CPL_P1P2;
                   segment = i_segment;
                   segment_fraction = new_along_plane_distance / interpolated_segment_length;
+                  current_angle = 0.5 * (interpolated_angle_top + interpolated_angle_bottom  - 2 * add_angle);
                   total_average_angle = (average_angle * total_length
-                                         + 0.5 * (interpolated_angle_top + interpolated_angle_bottom  - 2 * add_angle) * new_along_plane_distance);
+                                         +  current_angle * new_along_plane_distance);
                   total_average_angle = (std::fabs(total_average_angle) < std::numeric_limits<double>::epsilon() ? 0 : total_average_angle /
                                          (total_length + new_along_plane_distance));
                   depth_reference_surface = new_depth_reference_surface;
@@ -1059,6 +1065,7 @@ namespace WorldBuilder
       return_values.fraction_of_segment = segment_fraction;
       return_values.section = section;
       return_values.segment = segment;
+      return_values.angle = new_check_point_angle; // current_angle; // add_angle; // total_average_angle; // interpolated_angle_top;
       return_values.average_angle = total_average_angle;
       return_values.depth_reference_surface = depth_reference_surface;
       return_values.closest_trench_point = closest_point_on_line_cartesian;

tests/gwb-grid/oceanic_subduction_along_surface_velocity.grid

@@ -0,0 +1,14 @@
+# Grid File for creating vtk files to view in Paraview or Visit
+grid_type = cartesian
+dim = 2
+compositions = 1
+
+# domain of the grid (x and y are the surface, z is vertical with 0 at the bottom)
+x_min = 0e3
+x_max = 3000e3
+z_min = 0e3
+z_max = 1000e3
+
+# grid properties
+n_cell_x = 30
+n_cell_z = 10