add documentation and tests

include/world_builder/features/plume_models/temperature/gaussian.h

@@ -80,8 +80,6 @@ namespace WorldBuilder
 
           private:
             // Gaussian temperature submodule parameters
-            double min_depth;
-            double max_depth;
             std::vector<double> depths;
             std::vector<double> center_temperatures;
             std::vector<double> gaussian_sigmas;

source/world_builder/features/plume.cc

@@ -285,10 +285,8 @@ namespace WorldBuilder
                                                               surface_point);
 
       // If we are in the tip, we have to compute the difference diffently:
-      // TODO: move to utilities
       if (depth >= min_depth && depth < depths.front())
         {
-
           const double a = semi_major_axis_lengths.front();
           const double b = a * std::sqrt(1 - std::pow(eccentricity, 2));
           const double c = depths.front() - min_depth;
@@ -298,9 +296,7 @@ namespace WorldBuilder
           const double z = depths.front() - depth;
 
           // use ellipsoid equation:
-          relative_distance_from_center = std::pow(x, 2) / std::pow(a, 2)
-                                          + std::pow(y, 2) / std::pow(b, 2)
-                                          + std::pow(z, 2) / std::pow(c, 2);
+          relative_distance_from_center = (x*x)/(a*a) + (y*y)/(b*b) + (z*z)/(c*c);
         }
 
       if (depth <= max_depth && depth >= min_depth && relative_distance_from_center <= 1.)

source/world_builder/features/plume_models/temperature/gaussian.cc

@@ -39,8 +39,6 @@ namespace WorldBuilder
       {
         Gaussian::Gaussian(WorldBuilder::World *world_)
           :
-          min_depth(NaN::DSNAN),
-          max_depth(NaN::DSNAN),
           operation(Operations::REPLACE)
         {
           this->world = world_;
@@ -56,23 +54,34 @@ namespace WorldBuilder
           // Document plugin and require entries if needed.
           // Add `max distance fault` center to the required parameters.
           prm.declare_entry("", Types::Object({"centerline temperatures"}),
-                            "Gaussian temperature model. Can be set to use an adiabatic temperature at the boundaries.");
-
-          // Declare entries of this plugin
-          prm.declare_entry("min depth", Types::Double(0),
-                            "The depth in meters from which the temperature of the plume is present.");
-          prm.declare_entry("max depth", Types::Double(std::numeric_limits<double>::max()),
-                            "The depth in meters to which the temperature of the plume is present.");
+                            "Gaussian temperature model. The temperature is interpolated between the plume center "
+                            "and margin (as defined by the plume feature) using a Gaussian function: "
+                            "T(r) = T_center(z) exp(-r^2/(2 sigma^2). "
+                            "The temperature at the plume centerline T_center can be changed with depth by defining "
+                            "an array of depths and centerline temperatures, and temperature is interpolated linearly "
+                            "with depth. Similarly, the sigma of the Gaussian function (relative to the width of "
+                            "the plume as given by the plume feature) can be changed with depth. "
+                            "Temperature is always interpolated in a horizonzal/radial plane, except for the plume "
+                            "head: If the first depth of the plume centerline and the minimum depth of the plume "
+                            "feature are different, an ellipsoidal plume head is created in this depth range. "
+                            "Within this plume head, temperature is interpolated radially, i.e., depending on the "
+                            "distance from the center of the ellipsoid.");
 
           prm.declare_entry("depths", Types::Array(Types::Double(0)),
                             "The temperature at the center of this feature in degree Kelvin."
                             "If the value is below zero, the an adiabatic temperature is used.");
           prm.declare_entry("centerline temperatures", Types::Array(Types::Double(0)),
                             "The temperature at the center of this feature in degree Kelvin."
                             "If the value is below zero, the an adiabatic temperature is used.");
-          prm.declare_entry("gaussian sigmas", Types::Array(Types::Double(-1)),
+          prm.declare_entry("gaussian sigmas", Types::Array(Types::Double(0.3)),
                             "The sigma (standard deviation) of the Gaussian function used to compute the "
-                            "temperature distribution within the plume.");
+                            "temperature distribution within the plume. This sigma is non-dimensinal, i.e. "
+                            "it is defined relative to the distance between the plume center and margin as "
+                            "defined by the plume feature. Choosing a sigma of 1 therefore means that the "
+                            "temperature at the plume margin is set to a fraction of 1/sqrt(e) (approx. 0.61) "
+                            "of the centerline temperature. To achieve a smoother transition between the "
+                            "plume temperature and the outside temperature a smaller values has to be chosen "
+                            "for the gaussian sigmas.");
 
           // TODO: assert that the three have the same length
 
@@ -81,15 +90,18 @@ namespace WorldBuilder
         void
         Gaussian::parse_entries(Parameters &prm)
         {
-          min_depth = prm.get<double>("min depth");
-          max_depth = prm.get<double>("max depth");
-          WBAssert(max_depth >= min_depth, "max depth needs to be larger or equal to min depth.");
-
           operation = string_operations_to_enum(prm.get<std::string>("operation"));
 
           depths = prm.get_vector<double>("depths");
           center_temperatures = prm.get_vector<double>("centerline temperatures");
           gaussian_sigmas = prm.get_vector<double>("gaussian sigmas");
+
+          WBAssert(center_temperatures.size() == depths.size() && gaussian_sigmas.size() == depths.size(),
+                   "The depths, center_temperatures and gaussian_sigmas arrays need to have the same number of entries. "
+                   "At the moment there are: "
+                   << depths.size() << " depth entries, "
+                   << center_temperatures.size() << " centerline temperature entries, and "
+                   << gaussian_sigmas.size() << " gaussian sigma entries!");
         }
 
 
@@ -99,11 +111,11 @@ namespace WorldBuilder
                                   const double depth,
                                   const double gravity_norm,
                                   double temperature_,
-                                  const double /*feature_min_depth*/,
-                                  const double /*feature_max_depth*/,
+                                  const double feature_min_depth,
+                                  const double feature_max_depth,
                                   const double relative_distance_from_center) const
         {
-          if (depth <= max_depth && depth >= min_depth && relative_distance_from_center <= 1.)
+          if (depth <= feature_max_depth && depth >= feature_min_depth && relative_distance_from_center <= 1.)
             {
               // Figure out if the point is within the plume
               auto upper = std::upper_bound(depths.begin(), depths.end(), depth);

tests/gwb-dat/plume_gaussian.dat

@@ -0,0 +1,20 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 3
+# compositions = 2
+# x y z d T
+50e3 10e3 65e3 15e3
+50e3 30e3 65e3 15e3
+50e3 35e3 65e3 15e3
+50e3 40e3 65e3 15e3
+50e3 45e3 65e3 15e3
+50e3 50e3 65e3 15e3
+50e3 60e3 65e3 15e3
+50e3 10e3 0 100e3
+50e3 30e3 0 100e3
+50e3 35e3 0 100e3
+50e3 40e3 0 100e3
+50e3 45e3 0 100e3
+50e3 50e3 0 100e3
+50e3 60e3 0 100e3

tests/gwb-dat/plume_gaussian.wb

@@ -0,0 +1,38 @@
+{
+  "version": "0.6",
+  "surface temperature":273.15,
+  "potential mantle temperature":1573.15,
+  "features":
+  [
+    {"model":"mantle layer", "name":"upper mantle", "min depth":0, "max depth":100e3, "coordinates":[[-1e3,-1e3],[251e3,-1e3],[251e3,101e3],[-1e3,101e3]],
+     "temperature models":[{"model":"linear", "min depth":0, "max depth":100e3, "top temperature":1600, "bottom temperature":1601}]
+    },
+
+    {
+      "model":"plume", 
+      "name":"plume A", 
+      "min depth":10e2,
+      "max depth":150e3,
+      "coordinates":[[50e3, 50e3],[50e3, 50e3],[50e3, 50e3],[50e3, 50e3]],
+      "cross section depths":[40e3, 50e3, 75e3, 150e3],
+      "semi-major axis":[50e3, 35e3, 40e3, 45e3],
+      "eccentricity":[0.5, 0.5, 0.5, 0.5],
+      "rotation angles":[345, 355, 5, 15],
+      "temperature models":
+      [
+        {
+          "model":"gaussian",
+          "operation":"add",
+          "centerline temperatures":[200, 300, 400],
+          "gaussian sigmas":[0.3, 0.3, 0.3],
+          "depths":[40e3, 60e3, 150e3]
+        }
+      ], 
+
+      "composition models":
+      [
+        {"model":"uniform","compositions":[1], "min depth":0, "max depth":100e3}
+      ]
+    }
+  ]
+}

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

@@ -0,0 +1,15 @@
+# x y z d g T c0 c1 tag 
+50e3 10e3 65e3 15e3 1600.15 0 0  0
+50e3 30e3 65e3 15e3 1608.37 0 1  1
+50e3 35e3 65e3 15e3 1612.38 0 1  1
+50e3 40e3 65e3 15e3 1616.4 0 1  1
+50e3 45e3 65e3 15e3 1619.42 0 1  1
+50e3 50e3 65e3 15e3 1620.55 0 1  1
+50e3 60e3 65e3 15e3 1616.4 0 1  1
+50e3 10e3 0 100e3 1602.99 0 1  1
+50e3 30e3 0 100e3 1695.91 0 1  1
+50e3 35e3 0 100e3 1767.82 0 1  1
+50e3 40e3 0 100e3 1850.56 0 1  1
+50e3 45e3 0 100e3 1918.78 0 1  1
+50e3 50e3 0 100e3 1945.44 0 1  1
+50e3 60e3 0 100e3 1850.56 0 1  1

tests/gwb-dat/plume_gaussian_spherical.dat

@@ -0,0 +1,9 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 3
+# compositions = 2
+# x y z d T
+5.96291e+06 465152 1.09904e+06 395408.625
+5.9729e+06 518768 837228 416730.65625
+5.56627e+06 567790 850569 800000

tests/gwb-dat/plume_gaussian_spherical.wb

@@ -0,0 +1,38 @@
+{
+  "version": "0.6",
+  "coordinate system":{"model":"spherical", "depth method":"starting point"},
+  "maximum distance between coordinates":0.01,
+  "thermal expansion coefficient":0.0,
+  "features":
+  [
+    {"model":"mantle layer", "name":"upper mantle", "min depth":0, "max depth":1000e3, "coordinates":[[-1,-1],[41,-1],[41,-1],[-1,-1]],
+     "temperature models":[{"model":"linear", "min depth":0, "max depth":1000e3, "top temperature":1600, "bottom temperature":1600}]
+    },
+
+    {
+      "model":"plume", 
+      "name":"plume A", 
+      "min depth":10e2,
+      "coordinates":[[7.5, 7.5],[7.5, 8],[7.5,9]],
+      "cross section depths":[300e3, 500e3, 1000e3],
+      "semi-major axis":[7, 6, 7],
+      "eccentricity":[0.8, 0.8, 0.8],
+      "rotation angles":[355, 40, 90],
+      "temperature models":
+      [
+        {
+          "model":"gaussian",
+          "operation":"add",
+          "centerline temperatures":[200, 300],
+          "gaussian sigmas":[0.3, 0.3],
+          "depths":[200e3, 1000e3]
+        }
+      ], 
+
+      "composition models":
+      [
+        {"model":"uniform","compositions":[1], "min depth":10e2, "max depth":500e3}
+      ]
+    }
+  ]
+}

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

@@ -0,0 +1,4 @@
+# x y z d g T c0 c1 tag 
+5.96291e+06 465152 1.09904e+06 395408.625 1601.15 0 1  1
+5.9729e+06 518768 837228 416730.65625 1623.3 0 1  1
+5.56627e+06 567790 850569 800000 1776.42 0 0  1