Add gwb-dat tests

GeodynamicWorldBuilder · Mar 2, 2024 · 16a77d8 · 16a77d8
1 parent fe32a03
commit 16a77d8
Show file tree

Hide file tree

Showing 13 changed files with 269 additions and 32 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -11,7 +11,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 
 ## [Unreleased]
 ### Added
-- Implemented a method for modifying the mass conserving temperature model based on the movement of a spreading center through time. \[Daniel Douglas; Haoyuan Li; 2024-02-29; [#654](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/654)\]
+- Implemented the framework that will allow the mass conserving temperature model to account for the the movement of a spreading center through time. \[Daniel Douglas; Haoyuan Li; 2024-02-29; [#654](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/654)\]
 - Added an option to apply a spline in the mass conserving temperature of the slab. \[Haoyuan Li; 2024-02-27; [#659](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/659)\]
 - Added a cookbook tutorial for a simple subduction model in 2D Chunk geometry. \[Magali Billen; 2024-02-16; [#535](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/535)\]
 - Added a cookbook tutorial for a simple subduction model in 2D Cartesian geometry. \[Magali Billen; 2024-02-14; [#535](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/535)\]
@@ -109,7 +109,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - Changed the underlying function which determines when you are inside a slab or a fault. This may slightly change the result, but the new result is more realistic I think. \[Menno Fraters; 2021-05-01; [#229](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/229)\]
 - The Vizualizer now uses compressed output by default. This decreases the file size and increases performance. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\]
 - The Vizualizer buffers output before it writes it to a file. This increases performance. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\]
-- In the fault temperature model linear, the entiry top temperature is changed to center temperature and the entry bottom temperature is changed to side temperature, since this represents the actual sides more accurately. \[Menno Fraters; 2021-07-09; [#260](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/260)\]
+- In the fault temperature model linear, the entry top temperature is changed to center temperature and the entry bottom temperature is changed to side temperature, since this represents the actual sides more accurately. \[Menno Fraters; 2021-07-09; [#260](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/260)\]
 - A large overhoal of the distance_point_from_curved_planes function improving the accuracy in spherical coordinates. This may slightly change the results, but it should be an improvement both in accuracy and performance. Also makes available a coordinate system independent distrance function and some fast trigonometric functions. \[Menno Fraters; 2021-06-11; [#255](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/255)\]
 
 ### Fixed

diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h
@@ -455,9 +455,17 @@ namespace WorldBuilder
      * 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 mid_oceanic_spreading_velocities The spreading rate of the mid oceanic ridges at each ridge coordinate
      * @param coordinate_system The coordinate system
      * @param position_in_natural_coordinates_at_min_depth the current position in natural_coordinates
+     * @param subducting_plate_velocities the subducting plate velocities, currently this is only an optional parameter
+     * that can be used in the mass conserving subducting plate temperature model. This parameter allows the user to
+     * track the effect of ridge migration on the slab thermal structure
+     * @param ridge_migration_times the times that the corresponding section of the ridge has been moving, in years. This
+     * is used in combination with subducting_plate_velocities, and mid_oceanic_spreading_velocities to compute the distance
+     * that the spreading center has migrated. This vector is obtained from the input parameter "spreading velocity" in the
+     * mass conserving model when "spreading velocity" has the form: [ [t1,[[v11, v12, ...]], [t2,[[v21, v22, ...]], ... ].
+     * where tn is the time that ridge section n has been moving.
      * @return The content of the file.
     */
     std::vector<double>

diff --git a/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc b/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc
@@ -258,6 +258,24 @@ namespace WorldBuilder
 
           apply_spline = prm.get<bool>("apply spline");
           spline_n_points = prm.get<int>("number of points in spline");
+
+          if (subducting_velocities[0].size() > 1)
+            {
+              for (unsigned int ridge_index = 0; ridge_index < mid_oceanic_ridges.size(); ridge_index++)
+                {
+                  WBAssertThrow(subducting_velocities.size() == mid_oceanic_ridges.size() &&
+                                subducting_velocities[ridge_index].size() == mid_oceanic_ridges[ridge_index].size(),
+                                "subducting velocity must have the same dimension as the ridge coordinates parameter.");
+                  for (unsigned int point_index = 0; point_index < mid_oceanic_ridges[ridge_index].size(); point_index++)
+                    {
+                      WBAssertThrow(subducting_velocities[ridge_index][point_index] ==
+                                    ridge_spreading_velocities_at_each_ridge_point[ridge_index][point_index],
+                                    "Currently, subducting velocity must equal the spreading velocity to satisfy conservation of mass. "
+                                    "This will be changed in the future to allow for the spreading center to move depending on the relative "
+                                    "difference between the spreading velocity and the subducting velocity, thereby conserving mass.");
+                    }
+                }
+            }
         }
 
         double

diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc
@@ -1358,12 +1358,12 @@ namespace WorldBuilder
           // When subducting_velocities is input as an array, spatial variation
           if (subducting_plate_velocities[0].size() > 1)
             {
-              WBAssertThrow(subducting_plate_velocities.size() == mid_oceanic_ridges.size() && \
-                            subducting_plate_velocities[relevant_ridge].size() == mid_oceanic_ridges[relevant_ridge].size(),
-                            "subducting velocity and ridge coordinates must be the same dimension");
-              WBAssertThrow(ridge_migration_times.size() == mid_oceanic_ridges.size(),
-                            "the times for ridge migration specified in 'spreading velocity' must be the same dimension "
-                            "as ridge coordinates.");
+              WBAssert(subducting_plate_velocities.size() == mid_oceanic_ridges.size() && \
+                       subducting_plate_velocities[relevant_ridge].size() == mid_oceanic_ridges[relevant_ridge].size(),
+                       "subducting velocity and ridge coordinates must be the same dimension");
+              WBAssert(ridge_migration_times.size() == mid_oceanic_ridges.size(),
+                       "the times for ridge migration specified in 'spreading velocity' must be the same dimension "
+                       "as ridge coordinates.");
               subducting_velocity_point0 = subducting_plate_velocities[relevant_ridge][i_coordinate];
               subducting_velocity_point1 = subducting_plate_velocities[relevant_ridge][i_coordinate + 1];
               ridge_migration_time = ridge_migration_times[relevant_ridge];
@@ -1452,6 +1452,8 @@ namespace WorldBuilder
             double spreading_velocity_at_ridge_pt = spreading_velocity_at_ridge_pt1;
             double subducting_velocity_at_trench_pt = subducting_velocity_at_trench_pt1;
 
+            // This is required in spherical coordinates to ensure that the distance
+            // returned is the shortest distance around the sphere.
             if (compare_distance2 < compare_distance1)
               {
                 compare_distance = compare_distance2;
@@ -1476,41 +1478,23 @@ namespace WorldBuilder
       return result;
     }
 
-    // todo_effective
+    // TODO: implement method for modifying the age of the slab based on ridge/trench migration.
     std::vector<double>
     calculate_effective_trench_and_plate_ages(std::vector<double> ridge_parameters, double distance_along_plane)
     {
       WBAssert(ridge_parameters.size() == 4, "Internal error: ridge_parameters have the wrong size: " << ridge_parameters.size() << " instead of 4.");
       const double seconds_in_year = 60.0 * 60.0 * 24.0 * 365.25;  // sec/y
       const double spreading_velocity = ridge_parameters[0] * seconds_in_year; // m/yr
       double subducting_velocity = ridge_parameters[2] * seconds_in_year; // m/yr
-      const double ridge_velocity = spreading_velocity - subducting_velocity; // m/yr
-      const int sign_v = (ridge_velocity > 0) ? 1 : ((ridge_velocity < 0) ? -1 : 0);
 
-      const double time_of_ridge_migration = ridge_parameters[3]; //yr
-      const double ridge_drift_distance = abs(ridge_velocity) * time_of_ridge_migration; // m
-
-      double effective_age_shift = 0;
-      double trench_age_shift = 0;
       if (subducting_velocity <= 0)
         subducting_velocity = spreading_velocity;
 
-      else
-        {
-          if (distance_along_plane < ridge_drift_distance)
-            {
-              effective_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * \
-                                    distance_along_plane / subducting_velocity; // yr
-              trench_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * \
-                                 distance_along_plane / spreading_velocity; // yr
-            }
-        }
-
-      const double age_at_trench = ridge_parameters[1] / spreading_velocity + trench_age_shift * 0; // m/(m/y) = yr
+      const double age_at_trench = ridge_parameters[1] / spreading_velocity; // m/(m/y) = yr
       const double plate_age_sec = age_at_trench * seconds_in_year; // y --> seconds
 
       // Plate age increases with distance along the slab in the mantle
-      double effective_plate_age = plate_age_sec + (distance_along_plane / subducting_velocity + effective_age_shift) * seconds_in_year; // m/(m/y) = y(seconds_in_year)
+      double effective_plate_age = plate_age_sec + (distance_along_plane / subducting_velocity) * seconds_in_year; // m/(m/y) = y(seconds_in_year)
       WBAssertThrow(effective_plate_age >= 0, "The age of the subducting plate is less than or equal to 0. "
                     "Effective plate age: " << effective_plate_age);
       std::vector<double> result;
@@ -1519,7 +1503,6 @@ namespace WorldBuilder
       return result;
 
     }
-
   } // namespace Utilities
 } // namespace WorldBuilder
 

diff --git a/tests/data/type_data.json b/tests/data/type_data.json
@@ -13,6 +13,8 @@
 "value at points set ap val string": [[100,[[1,2],[3,4]]],[105],["test3",[[5,6]]],[300,[[-10,10]]]],
 "value at points set ap p1 string": [[100,[[1,2],[3,4]]],[105],[200,[[5,6]]],[300,[["test4",10]]]],
 "value at points set ap p2 string": [[100,[[1,2],[3,4]]],[105],[200,[[5,6]]],[300,[[-10,"test5"]]]],
+"value at array full": [[10,[[1, 2, 3]]], [20,[[4, 5]]]],
+"vector for vector or double": [[10, 20, 30, 40], [100, 200]],
 "unsigned int array" :[25,26,27],
 "bool array" :[true,false,true,false,true,false],
 "bool array nob" :[true,false,true,false,"True","false"],

diff --git a/tests/gwb-dat/cartesian_no_ridge_drift.dat b/tests/gwb-dat/cartesian_no_ridge_drift.dat
@@ -0,0 +1,20 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 2
+# compositions = 0
+# x y d T C0
+500e3 380e3 20e3
+525e3 380e3 20e3
+550e3 380e3 20e3
+575e3 380e3 20e3
+600e3 380e3 20e3
+625e3 380e3 20e3
+650e3 380e3 20e3
+675e3 380e3 20e3
+700e3 380e3 20e3
+600e3 360e3 40e3
+650e3 340e3 60e3
+700e3 320e3 80e3
+750e3 280e3 120e3
+800e3 240e3 160e3
diff --git a/tests/gwb-dat/cartesian_no_ridge_drift.wb b/tests/gwb-dat/cartesian_no_ridge_drift.wb
@@ -0,0 +1,87 @@
+{
+  "version":"0.6",
+  "coordinate system":{"model":"cartesian"},
+  "gravity model":{"model":"uniform", "magnitude":10},
+  "cross section":[[0,0],[10000e3,0]],"surface temperature":273, "force surface temperature":true,
+  "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, 
+  "specific heat":1000, "thermal diffusivity":1.0e-6,
+  "features":
+  [
+    {
+      "model": "oceanic plate",
+      "name": "Subducting Plate plate",
+      "max depth": 100e3,
+      "min depth": 0.0,
+      "coordinates": [ [0, 1000], [0, -1000], [500e3, -1000], [500e3, 1000] ],
+      "temperature models": [
+        {
+          "model": "half space model",
+          "min depth": 0.0,
+          "max depth": 100e3,
+          "spreading velocity": 0.05,
+          "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]]
+        }
+      ],
+      "composition models": 
+      [
+        {"model": "uniform", "min depth": 0.0, "max depth": 100e3, "compositions": [0]}
+      ]
+    },
+
+    {
+      "model": "oceanic plate",
+      "name": "Overriding Plate",
+      "max depth": 100e3,
+      "min depth": 0.0,
+      "coordinates": [ [1000e3, 1000], [1000e3, -1000], [500e3, -1000], [500e3, 1000] ],
+      "temperature models":
+      [
+        {
+          "model": "half space model", "min depth": 0.0, "max depth": 100e3, "spreading velocity": 0.03, 
+          "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]]
+        }
+      ]
+    },
+
+    {
+      "model": "subducting plate", "name": "initial slab", "coordinates": [[500e3, -1000.0], [500e3, 1000.0]],
+      "dip point": [400e6, 0],
+      "segments": [
+        {"length": 418880.0,
+          "thickness": [100e3],
+          "top truncation": [-100000.0],
+          "angle": [0,60],
+          "composition models": [
+            {
+              "model": "uniform",
+              "compositions": [0],
+              "max distance slab top": 100e3
+            }
+          ]
+        },
+        {
+          "length": 100000.0,
+          "thickness": [100e3],
+          "top truncation": [-100000.0],
+          "angle": [60, 60]
+        }
+      ],
+      "temperature models": [
+        {
+          "model": "mass conserving",
+          "density": 3300,
+          "thermal conductivity": 3.3,
+          "adiabatic heating": true,
+          "plate velocity": 0.05,
+          "subducting velocity": [[0.05, 0.05]],
+          "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]],
+          "coupling depth": 50000.0,
+          "taper distance": 100000.0,
+          "min distance slab top": -100000.0,
+          "max distance slab top": 100e3,
+          "reference model name": "half space model"
+        }
+      ]
+    }
+  ]
+}
diff --git a/tests/gwb-dat/cartesian_no_ridge_drift/screen-output.log b/tests/gwb-dat/cartesian_no_ridge_drift/screen-output.log
@@ -0,0 +1,15 @@
+# x z d T tag 
+500e3 380e3 20e3 1005.87  1
+525e3 380e3 20e3 969.259  1
+550e3 380e3 20e3 882.454  1
+575e3 380e3 20e3 744.612  1
+600e3 380e3 20e3 558.663  1
+625e3 380e3 20e3 635.761  1
+650e3 380e3 20e3 809.857  1
+675e3 380e3 20e3 814.935  1
+700e3 380e3 20e3 804.603  1
+600e3 360e3 40e3 1138.89  1
+650e3 340e3 60e3 1169.63  1
+700e3 320e3 80e3 1035.12  1
+750e3 280e3 120e3 1091.76  1
+800e3 240e3 160e3 926.557  1
diff --git a/tests/gwb-dat/spherical_ridge_distance_check.dat b/tests/gwb-dat/spherical_ridge_distance_check.dat
@@ -0,0 +1,9 @@
+# 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 20 50e3
+6321e3 22.5 22 50e3
diff --git a/tests/gwb-dat/spherical_ridge_distance_check.wb b/tests/gwb-dat/spherical_ridge_distance_check.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" :[[-380, 15],[-380, 30],[380, 30],[380, 15]],
+      "temperature models":[
+        {"model":"half space model", "min depth":0,  "max depth":100e3, "spreading velocity":0.05,
+        "top temperature":273,
+         "ridge coordinates": [[[-370,10], [-371, 16], [-373, 22.5]], [[371, 22.5], [375, 31]]]}]}      
+  ]
+}
diff --git a/tests/gwb-dat/spherical_ridge_distance_check/screen-output.log b/tests/gwb-dat/spherical_ridge_distance_check/screen-output.log
@@ -0,0 +1,3 @@
+# x y z d g T tag 
+6321e3 22.5 20 50e3 1045.9  0
+6321e3 22.5 22 50e3 1044.78  0
diff --git a/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt b/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt
@@ -108,4 +108,33 @@ TITLE
 [104] = 7
 [105] = 8
 [106] = 9
+[107] = 2
+[108] = 10
+[109] = 20
+[110] = 5
+[111] = 1
+[112] = 2
+[113] = 3
+[114] = 4
+[115] = 5
+[116] = 1
+[117] = 0
+[118] = 1
+[119] = 200
+[120] = 1
+[121] = 0
+[122] = 1
+[123] = 101
+[124] = 2
+[125] = 4
+[126] = 10
+[127] = 20
+[128] = 30
+[129] = 40
+[130] = 2
+[131] = 100
+[132] = 200
+[133] = 1
+[134] = 1
+[135] = 200