address comments and modify test

geodynamics · Oct 28, 2024 · 4470557 · 4470557
1 parent 0ca4143
commit 4470557
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 11 deletions.
diff --git a/doc/modules/changes/20241028_ranpengli b/doc/modules/changes/20241028_ranpengli
@@ -0,0 +1,4 @@
+Changed: The entropy reader throws when the range of the provided look-up table 
+does not fully cover the entropy-pressure range in the model.
+<br>
+(Ranpeng Li, 2024/10/28)
diff --git a/include/aspect/structured_data.h b/include/aspect/structured_data.h
@@ -177,6 +177,9 @@ namespace aspect
          * @param component The index (starting at 0) of the data column to be
          * returned. The index is therefore less than the number of data
          * columns in the data file (or specified in the constructor).
+         * @param crash_if_not_in_range If set to true, the function will throw
+         * when the requested position is outside the range of the coordinates
+         * provided by the data file.
          */
         double
         get_data(const Point<dim> &position,

diff --git a/source/structured_data.cc b/source/structured_data.cc
@@ -856,21 +856,35 @@ namespace aspect
     double
     StructuredDataLookup<dim>::get_data(const Point<dim> &position,
                                         const unsigned int component,
-                                        bool crash_if_not_in_range) const
+                                        const bool crash_if_not_in_range) const
     {
       Assert(component<n_components, ExcMessage("Invalid component index"));
 
       if (crash_if_not_in_range)
         {
-          const std::vector<double> x_coordinates = get_interpolation_point_coordinates(0);
-
-          Assert (position[0] >= x_coordinates[0] && position[0] <= x_coordinates[x_coordinates.size()-1],
-                  ExcMessage("The requested position is outside the range of the data."));
-
-          const std::vector<double> y_coordinates = get_interpolation_point_coordinates(1);
-
-          Assert (position[1] >= y_coordinates[0] && position[1] <= y_coordinates[y_coordinates.size()-1],
-                  ExcMessage("The requested position is outside the range of the data."));
+          const std::vector<double> &x_coordinates = get_interpolation_point_coordinates(0);
+
+          AssertThrow (position[0] >= x_coordinates[0] && position[0] <= x_coordinates[x_coordinates.size()-1],
+                       ExcMessage("The requested position "
+                                  + std::to_string(position[0])
+                                  + " is outside the range of the data (maximum value = "
+                                  + std::to_string(x_coordinates[0])
+                                  + " , minimum value = "
+                                  + std::to_string(x_coordinates[x_coordinates.size()-1])
+                                  + ")."
+                                 ));
+
+          const std::vector<double> &y_coordinates = get_interpolation_point_coordinates(1);
+
+          AssertThrow (position[1] >= y_coordinates[0] && position[1] <= y_coordinates[y_coordinates.size()-1],
+                       ExcMessage("The requested position "
+                                  + std::to_string(position[1])
+                                  + " is outside the range of the data (maximum value = "
+                                  + std::to_string(y_coordinates[0])
+                                  + " , minimum value = "
+                                  + std::to_string(y_coordinates[y_coordinates.size()-1])
+                                  + ")."
+                                 ));
         }
 
       return data[component]->value(position);

diff --git a/tests/entropy_plasticity.prm b/tests/entropy_plasticity.prm
@@ -108,7 +108,7 @@ subsection Initial composition model
     set Coordinate system = spherical
     set Variable names = r,phi
     set Function expression = 2535.08 + 10 * sin(2*phi)*sin(pi*(r-6371000)/2890000) - 1878.618 * erfc((6371000-r)/(2*sqrt(1.152e-6*5e7*31557600))) + \
-    486.368 * max(erfc((r-3481000)/(2*sqrt(1.152e-6*5e8*31557600))), exp(-((phi - 3/4*pi)^2)/(2*(1/35*pi)^2) - ((r-3481000)^2)/(2*400000^2))  );0
+    100 * max(erfc((r-3481000)/(2*sqrt(1.152e-6*5e8*31557600))), exp(-((phi - 3/4*pi)^2)/(2*(1/35*pi)^2) - ((r-3481000)^2)/(2*400000^2))  );0
 
     #set Function expression = 2687.748 + 10 * sin(2*phi)*sin(pi*(r-6371000)/2890000) - 2031.286 * erfc((6371000-r)/(2*sqrt(1.152e-6*5e7*31557600))) + 333.7 * erfc((r-3481000)/(2*sqrt(1.152e-6*5e8*31557600)));0
   end