Apply spline to the mass conserving temperature

GeodynamicWorldBuilder · Feb 27, 2024 · e52488f · e52488f
1 parent d9c2f64
commit e52488f
Show file tree

Hide file tree

Showing 8 changed files with 13,640 additions and 5,906 deletions.
diff --git a/doc/world_builder_declarations.schema.json b/doc/world_builder_declarations.schema.json
diff --git a/doc/world_builder_declarations_closed.md b/doc/world_builder_declarations_closed.md
diff --git a/doc/world_builder_declarations_open.md b/doc/world_builder_declarations_open.md
diff --git a/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h b/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h
@@ -131,6 +131,8 @@ namespace WorldBuilder
             };
             ReferenceModelName reference_model_name;
             const int plate_model_summation_number = 100; // for the plate model
+            bool apply_spline;
+            int spline_n_points;
         };
       } // namespace Temperature
     } // namespace SubductingPlateModels

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
@@ -34,6 +34,7 @@
 #include "world_builder/types/value_at_points.h"
 #include "world_builder/utilities.h"
 #include "world_builder/world.h"
+#include "world_builder/types/int.h"
 
 namespace WorldBuilder
 {
@@ -169,6 +170,15 @@ namespace WorldBuilder
           prm.declare_entry("reference model name",  Types::String("half space model"),
                             "The type of thermal model to use in the mass conserving model of slab temperature. "
                             "Options are half space model and plate model");
+
+          prm.declare_entry("apply spline",  Types::Bool(false),
+                            "Choose to apply a spline to smooth temperature where the two analytical equations "
+                            "meet at the minimum temperature.");
+
+          prm.declare_entry("number of points in spline", Types::Int(5),
+                            "The number of points on each side of the spline centered at the point where the two analytic equations meet."
+                            "A default value of 5, will lead to a total of 11 points in the spline.");
+
         }
 
         void
@@ -238,6 +248,9 @@ namespace WorldBuilder
             reference_model_name = plate_model;
           else if (reference_model_name_str=="half space model")
             reference_model_name = half_space_model;
+
+          apply_spline = prm.get<bool>("apply spline");
+          spline_n_points = prm.get<int>("number of points in spline");
         }
 
         double
@@ -473,8 +486,35 @@ namespace WorldBuilder
                       top_heat_content = top_heat_content * std::erfc(taper_con*theta);
                     }
 
-                  // Call the analytic solution to compute the temperature
-                  temperature = get_temperature_analytic(top_heat_content, min_temperature, background_temperature, temperature_, plate_velocity, effective_plate_age, adjusted_distance);
+                  double nondimensional_adjusted_distance = adjusted_distance / max_depth;
+
+                  if (apply_spline)
+                    {
+                      // A total number of (2 * spline_n_points + 1) points are picked,
+                      // spline_n_points points on each side and one at the center.
+                      // These points cover a range of (-1.0, 1.0) in nondimensional_adjusted_distance.
+                      Utilities::interpolation monotone_cubic_spline;
+
+                      const double interval_spline_distance = 1.0 / spline_n_points;
+                      std::vector<double> i_temperatures (2*(spline_n_points + 1), 0.0);
+
+                      for (int i = 0; i < 2 * spline_n_points + 1; ++i)
+                        {
+                          const double i_adjusted_distance = (i * interval_spline_distance - 1.0) * max_depth;
+                          const double i_temperature = get_temperature_analytic(top_heat_content, min_temperature, background_temperature, temperature_, plate_velocity, effective_plate_age, i_adjusted_distance);
+                          i_temperatures[i] = i_temperature;
+                        }
+
+                      monotone_cubic_spline.set_points(i_temperatures);
+
+                      const double index_distance = (nondimensional_adjusted_distance + 1.0) / interval_spline_distance;
+                      temperature =  monotone_cubic_spline(index_distance);
+                    }
+                  else
+                    {
+                      // Call the analytic solution to compute the temperature
+                      temperature = get_temperature_analytic(top_heat_content, min_temperature, background_temperature, temperature_, plate_velocity, effective_plate_age, adjusted_distance);
+                    }
                 }
               else
                 {

diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_spline.dat b/tests/gwb-dat/mass_conserving_slab_temperature_spline.dat
@@ -0,0 +1,18 @@
+# 2 profiles across a slab, one at 500 km depth
+# and one at 200 km depth
+# Now define parameters:
+# dim = 2
+# compositions = 0
+# x z d T 
+8000e3 0 500e3
+8050e3 0 500e3
+8100e3 0 500e3
+8150e3 0 500e3
+8200e3 0 500e3
+8300e3 0 500e3
+7800e3 0 200e3
+7850e3 0 200e3
+7900e3 0 200e3
+7950e3 0 200e3
+8000e3 0 200e3
+8200e3 0 200e3
diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_spline.wb b/tests/gwb-dat/mass_conserving_slab_temperature_spline.wb
@@ -0,0 +1,164 @@
+{
+  "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": "sp plate",
+      "max depth": 150000.0,
+      "min depth": -100000.0,
+      "coordinates": [
+        [
+          0.0,
+          -1000.0
+        ],
+        [
+          0.0,
+          1000.0
+        ],
+        [
+          7500000.0,
+          1000.0
+        ],
+        [
+          7500000.0,
+          -1000.0
+        ]
+      ],
+      "temperature models": [
+        {
+          "model": "plate model",
+          "min depth": -10000.0,
+          "max depth": 150000.0,
+          "spreading velocity": 0.05,
+          "ridge coordinates": [
+            [
+              [
+                0,
+                -1000.0
+              ],
+              [
+                0,
+                1000.0
+              ]
+            ]
+          ]
+        }
+      ],
+      "composition models": [
+        {
+          "model": "uniform",
+          "min depth": -10000.0,
+          "max depth": 7500.0,
+          "compositions": [
+            0
+          ]
+        },
+        {
+          "model": "uniform",
+          "min depth": 7500.0,
+          "max depth": 35200.0,
+          "compositions": [
+            1
+          ]
+        }
+      ]
+    },
+    {
+      "model": "subducting plate",
+      "name": "initial slab",
+      "coordinates": [
+        [
+          7500000.0,
+          -1000.0
+        ],
+        [
+          7500000.0,
+          1000.0
+        ]
+      ],
+      "dip point": [
+        40000000.0,
+        0.0
+      ],
+      "segments": [
+        {
+          "length": 1000000.0,
+          "thickness": [
+            300000.0
+          ],
+          "top truncation": [
+            -100000.0
+          ],
+          "angle": [
+            0,
+            80
+          ],
+          "composition models": [
+            {
+              "model": "uniform",
+              "compositions": [
+                0
+              ],
+              "max distance slab top": 7500.0
+            },
+            {
+              "model": "uniform",
+              "compositions": [
+                1
+              ],
+              "min distance slab top": 7500.0,
+              "max distance slab top": 35200.0
+            }
+          ]
+        },
+        {
+          "length": 300000.0,
+          "thickness": [
+            300000.0
+          ],
+          "top truncation": [
+            -100000.0
+          ],
+          "angle": [
+            80,
+            80
+          ]
+        }
+      ],
+      "temperature models": [
+        {
+          "model": "mass conserving",
+          "density": 3300,
+          "thermal conductivity": 3.3,
+          "adiabatic heating": true,
+          "plate velocity": 0.05,
+          "ridge coordinates": [
+            [
+              [
+                0,
+                -1000.0
+              ],
+              [
+                0,
+                1000.0
+              ]
+            ]
+          ],
+          "coupling depth": 50000.0,
+          "taper distance": 100000.0,
+          "min distance slab top": -100000.0,
+          "max distance slab top": 150000.0,
+          "reference model name": "plate model",
+          "apply spline": true,
+          "number of points in spline": 5
+        }
+      ]
+    }
+  ]
+}
diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_spline/screen-output.log b/tests/gwb-dat/mass_conserving_slab_temperature_spline/screen-output.log
@@ -0,0 +1,13 @@
+# x z d T tag 
+8000e3 0 500e3 1953.49  1
+8050e3 0 500e3 1652.36  1
+8100e3 0 500e3 1342.93  1
+8150e3 0 500e3 983.156  1
+8200e3 0 500e3 1486.98  1
+8300e3 0 500e3 1953.49  -1
+7800e3 0 200e3 1483.09  1
+7850e3 0 200e3 1272.26  1
+7900e3 0 200e3 1017.16  1
+7950e3 0 200e3 684.892  1
+8000e3 0 200e3 1385.27  1
+8200e3 0 200e3 1780.01  -1