add an option to make it crash when look-up table does not cover the …

…range, and manually set it for entropy method

doc/modules/changes/20250110_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, 2025/01/10)

include/aspect/structured_data.h

@@ -177,10 +177,14 @@ 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,
-                 const unsigned int component) const;
+                 const unsigned int component,
+                 const bool crash_if_not_in_range = false) const;
 
         /**
          * Returns the gradient of the function based on the bilinear

source/material_model/utilities.cc

@@ -837,7 +837,7 @@ namespace aspect
         EntropyReader::temperature(const double entropy,
                                    const double pressure) const
         {
-          const double temperature = material_lookup->get_data({entropy,pressure}, 0);
+          const double temperature = material_lookup->get_data({entropy,pressure}, 0, true);
           return temperature;
         }
 

source/structured_data.cc

@@ -862,9 +862,38 @@ namespace aspect
     template <int dim>
     double
     StructuredDataLookup<dim>::get_data(const Point<dim> &position,
-                                        const unsigned int component) const
+                                        const unsigned int component,
+                                        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);
+
+          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 (minimum value = "
+                                  + std::to_string(x_coordinates[0])
+                                  + " , maximum 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 (minimum value = "
+                                  + std::to_string(y_coordinates[0])
+                                  + " , maximum value = "
+                                  + std::to_string(y_coordinates[y_coordinates.size()-1])
+                                  + ")."
+                                 ));
+        }
+
       return data[component]->value(position);
     }
 

tests/entropy_initial_lookup.prm

@@ -5,7 +5,7 @@ set Dimension                              = 2
 set Use years in output instead of seconds = true
 set End time                               = 0
 set Nonlinear solver scheme                = iterated Advection and Stokes
-set Surface pressure                       = 0
+set Surface pressure                       = 25000
 set Adiabatic surface temperature          = 1600.0
 
 subsection Formulation

tests/entropy_initial_lookup/screen-output

@@ -10,56 +10,56 @@ Number of degrees of freedom: 3,217 (1,206+202+603+603+603)
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 1000+3 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0.000495634, 4.77467e-10, 0, 5.9652
-      Relative nonlinear residual (total system) after nonlinear iteration 1: 5.9652
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0.000495626, 4.77402e-10, 0, 5.96521
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 5.96521
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
-   Solving Stokes system (GMG)... 600+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 8.51534e-10, 5.64728e-10, 0, 0.0685014
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.0685014
+   Solving Stokes system (GMG)... 676+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 8.51395e-10, 5.64621e-10, 0, 0.0685012
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.0685012
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
-   Solving Stokes system (GMG)... 500+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.01722e-09, 4.2244e-10, 0, 0.0168383
-      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.0168383
+   Solving Stokes system (GMG)... 550+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.01699e-09, 4.22363e-10, 0, 0.0168382
+      Relative nonlinear residual (total system) after nonlinear iteration 3: 0.0168382
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
-   Solving Stokes system (GMG)... 286+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.51013e-10, 3.95146e-10, 0, 0.00246445
+   Solving Stokes system (GMG)... 245+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.509e-10, 3.95043e-10, 0, 0.00246445
       Relative nonlinear residual (total system) after nonlinear iteration 4: 0.00246445
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 2 iterations.
    Copying properties into prescribed compositional field density_field... done.
-   Solving Stokes system (GMG)... 194+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.21895e-10, 4.19474e-10, 0, 0.000253486
-      Relative nonlinear residual (total system) after nonlinear iteration 5: 0.000253486
+   Solving Stokes system (GMG)... 190+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.21723e-10, 4.19413e-10, 0, 0.00025346
+      Relative nonlinear residual (total system) after nonlinear iteration 5: 0.00025346
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 76+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.79118e-10, 3.18101e-10, 0, 2.03324e-05
-      Relative nonlinear residual (total system) after nonlinear iteration 6: 2.03324e-05
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 7.79017e-10, 3.17962e-10, 0, 2.03204e-05
+      Relative nonlinear residual (total system) after nonlinear iteration 6: 2.03204e-05
 
    Copying properties into prescribed temperature field... done.
    Solving temperature system... 3 iterations.
    Solving entropy system ... 3 iterations.
    Copying properties into prescribed compositional field density_field... done.
-   Solving Stokes system (GMG)... 55+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 5.78613e-10, 3.56974e-10, 0, 1.35817e-06
-      Relative nonlinear residual (total system) after nonlinear iteration 7: 1.35817e-06
+   Solving Stokes system (GMG)... 56+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 5.78371e-10, 3.56725e-10, 0, 1.35606e-06
+      Relative nonlinear residual (total system) after nonlinear iteration 7: 1.35606e-06
 
 
    Postprocessing:

tests/entropy_initial_lookup/statistics

@@ -0,0 +1,32 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Iterations for temperature solver
+# 10: Iterations for composition solver 1
+# 11: Iterations for composition solver 2
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Visualization file name
+# 16: RMS velocity (m/year)
+# 17: Max. velocity (m/year)
+# 18: Minimal temperature (K)
+# 19: Average temperature (K)
+# 20: Maximal temperature (K)
+# 21: Average nondimensional temperature (K)
+# 22: Outward mass flux through boundary with indicator 0 ("left") (kg/yr)
+# 23: Outward mass flux through boundary with indicator 1 ("right") (kg/yr)
+# 24: Outward mass flux through boundary with indicator 2 ("bottom") (kg/yr)
+# 25: Outward mass flux through boundary with indicator 3 ("top") (kg/yr)
+# 26: Minimal value for composition entropy
+# 27: Maximal value for composition entropy
+# 28: Global mass for composition entropy
+# 29: Minimal value for composition density_field
+# 30: Maximal value for composition density_field
+# 31: Global mass for composition density_field
+0 0.000000000000e+00 0.000000000000e+00 100 1408 603 1206 7 25 20 4294967289 2790 3108 2827 output-entropy_initial_lookup/solution/solution-00000 7.88044505e-03 9.98412057e-03 1.59574424e+03 1.60014975e+03 1.60309788e+03 1.16764184e-03 0.00000000e+00 0.00000000e+00 4.08138251e+05 -4.08137893e+05 2.19666723e+03 2.53347698e+03 3.95671887e+13 3.13952225e+03 4.77911946e+03 6.90810628e+13 

tests/entropy_initial_lookup_wb.prm

@@ -8,7 +8,7 @@ set End time = 0
 set Output directory = output_3
 set Timing output frequency = 10
 set Pressure normalization = surface
-set Surface pressure = 0
+set Surface pressure = 25000
 set Adiabatic surface temperature = 1573.0
 set Resume computation = false
 set Nonlinear solver scheme                = iterated Advection and Stokes

tests/entropy_initial_lookup_wb/screen-output

@@ -39,36 +39,36 @@ Number of degrees of freedom: 47,413 (18,018+2,368+9,009+9,009+9,009)
    Solving entropy system ... 0 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 33+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.60794e-16, 1.69272e-16, 0, 0.00433026
-      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.00433026
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.5721e-16, 1.61845e-16, 0, 0.00433027
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.00433027
 
    Solving temperature system... 0 iterations.
    Solving entropy system ... 0 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 28+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.60794e-16, 1.69272e-16, 0, 0.00148712
-      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.00148712
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.5721e-16, 1.61845e-16, 0, 0.00148713
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 0.00148713
 
    Solving temperature system... 0 iterations.
    Solving entropy system ... 0 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 23+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.60794e-16, 1.69272e-16, 0, 6.99073e-05
-      Relative nonlinear residual (total system) after nonlinear iteration 3: 6.99073e-05
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.5721e-16, 1.61845e-16, 0, 6.99077e-05
+      Relative nonlinear residual (total system) after nonlinear iteration 3: 6.99077e-05
 
    Solving temperature system... 0 iterations.
    Solving entropy system ... 0 iterations.
    Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 18+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.60794e-16, 1.69272e-16, 0, 5.0992e-06
-      Relative nonlinear residual (total system) after nonlinear iteration 4: 5.0992e-06
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.5721e-16, 1.61845e-16, 0, 5.09921e-06
+      Relative nonlinear residual (total system) after nonlinear iteration 4: 5.09921e-06
 
 
    Postprocessing:
      Writing graphical output:           output-entropy_initial_lookup_wb/solution/solution-00003
      RMS, max velocity:                  0.0138 m/year, 0.124 m/year
      Temperature min/avg/max:            273 K, 2184 K, 3500 K
-     Mass fluxes through boundary parts: 0.004614 kg/yr, 0.06578 kg/yr, 0 kg/yr, 2.212e-08 kg/yr
+     Mass fluxes through boundary parts: 0.004614 kg/yr, 0.06577 kg/yr, 0 kg/yr, 2.203e-08 kg/yr
      Compositions min/max/mass:          577.9/2758/1.089e+17 // 3228/5592/1.939e+17
 
 Termination requested by criterion: end time

tests/entropy_initial_lookup_wb/statistics

@@ -0,0 +1,35 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Visualization file name
+# 10: RMS velocity (m/year)
+# 11: Max. velocity (m/year)
+# 12: Minimal temperature (K)
+# 13: Average temperature (K)
+# 14: Maximal temperature (K)
+# 15: Average nondimensional temperature (K)
+# 16: Outward mass flux through boundary with indicator 0 ("bottom") (kg/yr)
+# 17: Outward mass flux through boundary with indicator 1 ("top") (kg/yr)
+# 18: Outward mass flux through boundary with indicator 2 ("west") (kg/yr)
+# 19: Outward mass flux through boundary with indicator 3 ("east") (kg/yr)
+# 20: Minimal value for composition entropy
+# 21: Maximal value for composition entropy
+# 22: Global mass for composition entropy
+# 23: Minimal value for composition density_field
+# 24: Maximal value for composition density_field
+# 25: Global mass for composition density_field
+# 26: Iterations for temperature solver
+# 27: Iterations for composition solver 1
+# 28: Iterations for composition solver 2
+# 29: Iterations for Stokes solver
+# 30: Velocity iterations in Stokes preconditioner
+# 31: Schur complement iterations in Stokes preconditioner
+0 0.000000000000e+00 0.000000000000e+00  384  3739 1649  3298 0 output-entropy_initial_lookup_wb/solution/solution-00000 0.00000000e+00 0.00000000e+00 2.73000000e+02 2.17080114e+03 3.50000000e+03 5.88100755e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.77925312e+02 2.75762818e+03 1.07732840e+17 0.00000000e+00 0.00000000e+00 0.00000000e+00 0 0          0  0   0   0 
+0 0.000000000000e+00 0.000000000000e+00  684  7217 3189  6378 0 output-entropy_initial_lookup_wb/solution/solution-00001 0.00000000e+00 0.00000000e+00 2.73000000e+02 2.18276710e+03 3.50000000e+03 5.91808832e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.77925312e+02 2.75762818e+03 1.08621269e+17 0.00000000e+00 0.00000000e+00 0.00000000e+00 0 0          0  0   0   0 
+0 0.000000000000e+00 0.000000000000e+00 1023 10982 4851  9702 0 output-entropy_initial_lookup_wb/solution/solution-00002 0.00000000e+00 0.00000000e+00 2.73000000e+02 2.18453715e+03 3.50000000e+03 5.92357344e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.77925312e+02 2.75762818e+03 1.08867722e+17 0.00000000e+00 0.00000000e+00 0.00000000e+00 0 0          0  0   0   0 
+0 0.000000000000e+00 0.000000000000e+00 1950 20386 9009 18018 4 output-entropy_initial_lookup_wb/solution/solution-00003 1.37687774e-02 1.23769331e-01 2.73000000e+02 2.18434812e+03 3.50000000e+03 5.92298767e-01 4.61391962e-03 6.57718001e-02 0.00000000e+00 2.20276789e-08 5.77925312e+02 2.75762818e+03 1.08896595e+17 3.22822228e+03 5.59177718e+03 1.93948135e+17 0 0 4294967292 98 106 106 

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