Read in nonuniform grid

data/compositional-initial-conditions/ascii-data/test/box_2d_nonuniform.txt

@@ -0,0 +1,13 @@
+# Test data for ascii data compositional initial conditions.
+# Only next line is parsed in format: [nx] [ny] because of keyword "POINTS:"
+# POINTS: 3 3
+# Columns: x y composition_1
+0.      0.      0.
+200000  0.      1.
+660000  0.      0.
+0.      330000  1.
+200000  330000  1.
+660000  330000  1.
+0.      660000  0.
+200000  660000  1.
+660000  660000  0.

doc/modules/changes.h

@@ -6,6 +6,11 @@
  *
  * <ol>
  *
+ * <li> Changed: It is now possible to read in ascii data files of
+ * which the coordinates are not equally spaced.
+ * <br>
+ * (Anne Glerum, 2016/08/05)
+ *
  * <li> New: It is now possible to create compositional fields that are
  * not advected by a field method, but interpolated from particle properties.
  * This simplifies the process of using the particles as 'active' particles

include/aspect/utilities.h

@@ -411,16 +411,23 @@ namespace aspect
 
         /**
          * Interpolation functions to access the data.
+         * Either InterpolatedUniformGridData or InterpolatedTensorProductGridData;
+         * the type is determined from the grid specified in the data file.
          */
-        std::vector<Functions::InterpolatedUniformGridData<dim> *> data;
+        std::vector<Function<dim> *> data;
 
         /**
-         * Min and Max coordinates in data file
+         * The coordinate values in each direction as specified in the data file.
+         */
+        std_cxx11::array<std::vector<double>,dim> coordinate_values;
+
+        /**
+         * The min and max of the coordinates in the data file.
          */
         std_cxx11::array<std::pair<double,double>,dim> grid_extent;
 
         /**
-         * Number of points in the data grid.
+         * Number of points in the data grid as specified in the data file.
          */
         TableIndices<dim> table_points;
 

source/utilities.cc

@@ -928,31 +928,95 @@ namespace aspect
 
           i++;
 
-          // TODO: add checks for coordinate ordering in data files
         }
 
       AssertThrow(i == (components + dim) * data_table.n_elements(),
                   ExcMessage (std::string("Number of read in points does not match number of expected points. File corrupted?")));
 
+      // In case the data is specified on a grid that is equidistant
+      // in each coordinate direction, we only need to store
+      // (besides the data) the number of intervals in each direction and
+      // the begin- and endpoints of the coordinates.
+      // In case the grid is not equidistant, we need to keep
+      // all the coordinates in each direction, which is more costly.
+      // Here we fill the data structures needed for both cases,
+      // and check whether the coordinates are equidistant or not.
+      // We also check the requirement that the coordinates are
+      // strictly ascending.
+
+      // The number of intervals in each direction
       std_cxx11::array<unsigned int,dim> table_intervals;
 
+      // Whether or not the grid is equidistant
+      bool equidistant_grid = true;
+
       for (unsigned int i = 0; i < dim; i++)
         {
           table_intervals[i] = table_points[i] - 1;
 
           TableIndices<dim> idx;
-          grid_extent[i].first = data_tables[i](idx);
-          idx[i] = table_points[i] - 1;
-          grid_extent[i].second = data_tables[i](idx);
+          double temp_coord = data_tables[i](idx);
+          double new_temp_coord = 0;
+
+          // The minimum coordinates
+          grid_extent[i].first = temp_coord;
+
+          // The first coordinate value
+          coordinate_values[i].push_back(temp_coord);
+
+          // The grid spacing
+          double first_grid_spacing;
+          double grid_spacing;
+
+          // Loop over the rest of the coordinate points
+          for (unsigned int n = 1; n < table_points[i]; n++)
+            {
+              idx[i] = n;
+              new_temp_coord = data_tables[i](idx);
+              AssertThrow(new_temp_coord > temp_coord,
+                          ExcMessage ("Coordinates in dimension "
+                                      + int_to_string(i)
+                                      + " are not strictly ascending. "));
+
+              // Test whether grid is equidistant
+              if (n == 1)
+                first_grid_spacing = new_temp_coord - temp_coord;
+              else
+                {
+                  grid_spacing = new_temp_coord - temp_coord;
+                  // Compare current grid spacing with first grid spacing,
+                  // taking into account roundoff of the read-in coordinates
+                  if (std::abs(grid_spacing - first_grid_spacing) > 0.005*(grid_spacing+first_grid_spacing))
+                    equidistant_grid = false;
+                }
+
+              // Set the coordinate value
+              coordinate_values[i].push_back(new_temp_coord);
+              temp_coord = new_temp_coord;
+            }
+
+          // The maximum coordinate
+          grid_extent[i].second = temp_coord;
         }
 
+      // For each data component, set up a GridData,
+      // its type depending on the read-in grid.
       for (unsigned int i = 0; i < components; i++)
         {
           if (data[i])
             delete data[i];
-          data[i] = new Functions::InterpolatedUniformGridData<dim> (grid_extent,
-                                                                     table_intervals,
-                                                                     data_tables[dim+i]);
+
+          if (equidistant_grid)
+            data[i] = new Functions::InterpolatedUniformGridData<dim> (grid_extent,
+                                                                       table_intervals,
+                                                                       data_tables[dim+i]);
+          else
+            {
+              if (Utilities::MPI::this_mpi_process(comm) == 0)
+                std::cout << "   Ascii data file coordinates are not equidistant. " << std::endl << std::endl;
+              data[i] = new Functions::InterpolatedTensorProductGridData<dim> (coordinate_values,
+                                                                               data_tables[dim+i]);
+            }
         }
     }
 

tests/ascii_data_nonuniform_initial_composition_2d_box.prm

@@ -0,0 +1,101 @@
+##### simple test for ascii data initial composition
+
+set Dimension                              = 2
+
+set Use years in output instead of seconds = true
+set End time                               = 1e6
+
+set Adiabatic surface temperature          = 1613.0
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 660000
+    set Y extent = 660000
+  end
+end
+
+subsection Initial conditions
+  set Model name = adiabatic
+  subsection Adiabatic
+    set Amplitude = 300
+    set Radius    = 250000
+  end
+end
+
+
+subsection Compositional fields
+  set Number of fields = 1
+end
+
+subsection Compositional initial conditions
+  set Model name = ascii data
+
+  subsection Ascii data model
+    set Data directory       = $ASPECT_SOURCE_DIR/data/compositional-initial-conditions/ascii-data/test/
+    set Data file name       = box_2d_nonuniform.txt
+  end
+end
+
+subsection Boundary temperature model
+  set Model name = initial temperature
+end
+
+subsection Boundary composition model
+  set Model name = initial composition
+end
+
+subsection Model settings
+  set Fixed temperature boundary indicators   = top,bottom,left
+  set Fixed composition boundary indicators   = left
+
+  set Zero velocity boundary indicators       =
+  set Prescribed velocity boundary indicators = bottom:function,left:function,right:function,top:function
+
+  set Tangential velocity boundary indicators = 
+
+  set Include adiabatic heating               = false
+  set Include shear heating                   = false
+end
+
+
+subsection Boundary velocity model
+  subsection Function
+    set Function expression = 1;0
+  end
+end
+
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+
+subsection Material model
+  set Model name = simple
+  subsection Simple model
+    set Viscosity = 1e21
+  end
+end
+
+
+subsection Mesh refinement
+  set Initial global refinement                = 2
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+  set Strategy                                 = temperature
+end
+
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, composition statistics
+  subsection Visualization
+    set Time between graphical output = 0
+  end
+end
+