Merge branch 'master' into new-melt-solver-rebase

README.md

@@ -34,6 +34,9 @@ output of ASPECT can also be found in the ASPECT
 [manual](http://www.math.clemson.edu/~heister/manual.pdf). This manual also
 discusses the structure of the source code.
 
+For getting started, you can also watch our online 
+[tutorial](https://geodynamics.org/cig/events/calendar/2016-cig-all-hands-meeting/aspect-tutorial/tutorial/).
+
 
 
 Contributing to ASPECT

benchmarks/spiegelman_et_al_2016_gcubed/spiegelman_material.cc

@@ -79,7 +79,7 @@ namespace aspect
      */
 
     template <int dim>
-    class DruckerPragerCompositions : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
+    class SpiegelmanMaterial : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
     {
       public:
         std::vector<double> compute_volume_fractions( const std::vector<double> &compositional_fields) const;
@@ -203,7 +203,7 @@ namespace aspect
   {
     template <int dim>
     std::vector<double>
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     compute_volume_fractions( const std::vector<double> &compositional_fields) const
     {
       std::vector<double> volume_fractions( compositional_fields.size()+1);
@@ -235,7 +235,7 @@ namespace aspect
 
     template <int dim>
     double
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     compute_second_invariant(const SymmetricTensor<2,dim> strain_rate, const double min_strain_rate) const
     {
       const double edot_ii_strict = std::sqrt(strain_rate*strain_rate);
@@ -245,7 +245,7 @@ namespace aspect
 
     template <int dim>
     double
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     compute_viscosity(const double edot_ii,
                       const double pressure,
                       const int comp,
@@ -282,7 +282,7 @@ namespace aspect
 
     template <int dim>
     void
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
              MaterialModel::MaterialModelOutputs<dim> &out) const
     {
@@ -479,31 +479,31 @@ namespace aspect
 
     template <int dim>
     double
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     reference_viscosity () const
     {
       return ref_visc;
     }
 
     template <int dim>
     double
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     reference_density () const
     {
       return densities[0];
     }
 
     template <int dim>
     bool
-    DruckerPragerCompositions<dim>::
+    SpiegelmanMaterial<dim>::
     is_compressible () const
     {
       return (reference_compressibility != 0);
     }
 
     template <int dim>
     void
-    DruckerPragerCompositions<dim>::declare_parameters (ParameterHandler &prm)
+    SpiegelmanMaterial<dim>::declare_parameters (ParameterHandler &prm)
     {
       prm.enter_subsection("Compositional fields");
       {
@@ -601,7 +601,7 @@ namespace aspect
 
     template <int dim>
     void
-    DruckerPragerCompositions<dim>::parse_parameters (ParameterHandler &prm)
+    SpiegelmanMaterial<dim>::parse_parameters (ParameterHandler &prm)
     {
       using namespace Utilities;
       // can't use this->n_compositional_fields(), because some
@@ -697,7 +697,7 @@ namespace aspect
 {
   namespace MaterialModel
   {
-    ASPECT_REGISTER_MATERIAL_MODEL(DruckerPragerCompositions,
+    ASPECT_REGISTER_MATERIAL_MODEL(SpiegelmanMaterial,
                                    "spiegelman 2016",
                                    "An implementation of the spiegelman 2016 benchmark paper in gcubed "
                                    "(doi:10.1002/ 2015GC006228). It implements a regularized Drucker Prager "

doc/modules/changes/20180321_gassmoeller

@@ -0,0 +1,4 @@
+Improved: The 'depth average' postprocessor was significantly optimized, and
+now uses approximately ten times less computing time.
+<br>
+(Rene Gassmoeller, 2018/03/21)

doc/update_prm_files_to_2.0.0.sed

@@ -12,10 +12,10 @@
 # sed -i "" -f update_source_files_to_2.0.0.sed FILENAME
 
 # Rename compositional initial conditions
-s/subsection Compositional initial conditions/Initial composition model/g
+s/subsection Compositional initial conditions/subsection Initial composition model/g
 
 # Rename initial (temperature) conditions
-s/subsection Initial conditions/Initial temperature model/g
+s/subsection Initial conditions/subsection Initial temperature model/g
 
 # Rename adiabatic conditions plugin initial profile
 s/subsection Initial profile/subsection Compute profile/g

include/aspect/lateral_averaging.h

@@ -24,10 +24,68 @@
 
 #include <aspect/simulator_access.h>
 
+#include <deal.II/fe/fe_values.h>
+
 namespace aspect
 {
   using namespace dealii;
 
+  namespace internal
+  {
+    /**
+     * This is the base class for all the functors implemented. Each is
+     * used to compute one of the properties that will be laterally
+     * averaged. The point of the base class is to allow handing over
+     * a variable of type
+     * <code> std::vector<std::unique_ptr<FunctorBase<dim> > > </code> to the
+     * LateralAveraging::get_averages() function.
+     */
+    template <int dim>
+    class FunctorBase
+    {
+      public:
+        /**
+         * operator() will have @p in and @p out filled out if @p true. By default
+         * returns false.
+         */
+        virtual
+        bool
+        need_material_properties() const;
+
+        /**
+         * If this functor needs additional material model outputs
+         * create them in here. By default, this does nothing.
+          */
+        virtual
+        void
+        create_additional_material_model_outputs (const unsigned int n_points,
+                                                  MaterialModel::MaterialModelOutputs<dim> &outputs) const;
+
+        /**
+         * Called once at the beginning of compute_lateral_averages() to setup
+         * internal data structures with the number of quadrature points.
+         */
+        virtual
+        void
+        setup(const unsigned int q_points);
+
+        /**
+        * This takes @p in material model inputs and @p out outputs (which are filled
+        * if need_material_properties() == true), an initialized FEValues
+        * object for a cell, and the current solution vector as inputs.
+        * Functions in derived classes should then evaluate the desired quantity
+        * and return the results in the output vector, which is q_points long.
+        */
+        virtual
+        void
+        operator()(const MaterialModel::MaterialModelInputs<dim> &in,
+                   const MaterialModel::MaterialModelOutputs<dim> &out,
+                   const FEValues<dim> &fe_values,
+                   const LinearAlgebra::BlockVector &solution,
+                   std::vector<double> &output) = 0;
+    };
+  }
+
   /**
    * LateralAveraging is a class that performs various averaging operations
    * on the solution.  The functions of this class take a vector as an argument.
@@ -46,6 +104,23 @@ namespace aspect
   class LateralAveraging : public SimulatorAccess<dim>
   {
     public:
+      /**
+       * Fill the @p values with a set of lateral averages of the selected
+       * @p property_names. See the implementation of this function for
+       * a range of accepted names. This function is more efficient than
+       * calling multiple of the other functions that compute one property
+       * each.
+       *
+       * @param n_slices The number of depth slices to perform the averaging in.
+       * @return The output vector of laterally averaged values. Each vector
+       * has the same size of @p n_slices, and there are
+       * as many vectors returned as names in @p property_names.
+       * @param property_names Names of the available quantities to average.
+       * Check the implementation of this function for available names.
+       */
+      std::vector<std::vector<double> >
+      get_averages(const unsigned int n_slices,
+                   const std::vector<std::string> &property_names) const;
 
       /**
        * Fill the argument with a set of lateral averages of the current
@@ -137,61 +212,24 @@ namespace aspect
       get_vertical_heat_flux_averages(std::vector<double> &values) const;
 
     private:
-
       /**
-       * Internal routine to compute the depth average of a certain quantity.
-       *
-       * The functor @p fctr must be an object of a user defined type that can
-       * be arbitrary but has to satisfy certain requirements. In essence,
-       * this class type needs to implement the following interface of member
-       * functions:
-       * @code
-       * template <int dim>
-       * class Functor
-       * {
-       *   public:
-       *     // operator() will have @p in and @p out filled out if @p true
-       *     bool need_material_properties() const;
-       *
-       *     // called once at the beginning with the number of quadrature points
-       *     void setup(const unsigned int q_points);
-       *
-       *     // Fill @p output for each quadrature point.
-       *     // This takes material model inputs and outputs (which are filled
-       *     // if need_material_properties() == true), an initialized FEValues
-       *     // object for a cell, and the current solution vector as inputs.
-       *     // It then evaluates the desired quantity and puts the results in
-       *     // the output vector, which is q_points long.
-       *     void operator()(const MaterialModel::MaterialModelInputs<dim> &in,
-       *                     const MaterialModel::MaterialModelOutputs<dim> &out,
-       *                     FEValues<dim> &fe_values,
-       *                     const LinearAlgebra::BlockVector &solution,
-       *                     std::vector<double> &output);
-       * };
-       * @endcode
+       * Internal routine to compute the depth averages of several quantities.
+       * All of the public functions that compute a single field also call this
+       * function. The vector of functors @p functors must contain one or more
+       * objects of classes that are derived from FunctorBase and are used to
+       * fill the values vectors.
        *
-       * @param values The output vector of depth averaged values. The
-       * function takes the pre-existing size of this vector as the number of
-       * depth slices.
-       * @param fctr Instance of a class satisfying the signature above.
+       * @param n_slices Number of depth slices to be computed.
+       * @param functors Instances of a class derived from FunctorBase
+       * that are used to compute the averaged properties.
+       * @return The output vectors of depth averaged values. The
+       * function returns one vector of doubles per property and uses
+       * @p n_slices as the number of depth slices.
+       * Each returned vector has the same size.
        */
-      template <class FUNCTOR>
-      void compute_lateral_average(std::vector<double> &values,
-                                   FUNCTOR &fctr) const;
-
-      /**
-       * Compute a depth average of the current temperature/composition. The
-       * function fills a vector that contains average
-       * temperatures/compositions over slices of the domain of same depth.
-       *
-       * @param field  Extractor for temperature or compositional field to average.
-       * @param values The output vector of depth averaged values. The
-       * function takes the pre-existing size of this vector as the number of
-       * depth slices.
-       */
-      void get_field_averages(const FEValuesExtractors::Scalar &field,
-                              std::vector<double> &values) const;
-
+      std::vector<std::vector<double> >
+      compute_lateral_averages(const unsigned int n_slices,
+                               std::vector<std_cxx11::unique_ptr<internal::FunctorBase<dim> > > &functors) const;
   };
 }
 

include/aspect/plugins.h

@@ -27,6 +27,9 @@
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/parameter_handler.h>
 #include <deal.II/base/std_cxx11/tuple.h>
+#include <deal.II/base/exceptions.h>
+
+#include <boost/core/demangle.hpp>
 
 #include <string>
 #include <list>
@@ -40,6 +43,55 @@ namespace aspect
 {
   template <int dim> class SimulatorAccess;
 
+  namespace Plugins
+  {
+    using namespace dealii;
+
+    /**
+     * This function returns if a given plugin (e.g. a material model returned
+     * from SimulatorAccess::get_material_model() ) matches a certain plugin
+     * type (e.g. MaterialModel::Simple). This check is needed, because often
+     * it is only possible to get a reference to an Interface, not the actual
+     * plugin type, but the actual plugin type might be important. For example
+     * a radial gravity model might only be implemented for spherical geometry
+     * models, and would want to check if the current geometry is in fact a
+     * spherical shell.
+     */
+    template <typename TestType, typename PluginType>
+    inline
+    bool
+    plugin_type_matches (const PluginType &object)
+    {
+      return (dynamic_cast<const TestType *> (&object) != NULL);
+    }
+
+    /**
+     * This function converts a reference to a type (in particular a reference
+     * to an interface class) into a reference to a different type (in
+     * particular a plugin class). This allows accessing members of the plugin
+     * that are not specified in the interface class. Note that you should
+     * first check if the plugin type is actually convertible by calling
+     * plugin_matches_type() before calling this function. If the plugin is
+     * not convertible this function throws an exception.
+     */
+    template <typename TestType, typename PluginType>
+    inline
+    TestType &
+    get_plugin_as_type (PluginType &object)
+    {
+      AssertThrow(plugin_type_matches<TestType>(object),
+                  ExcMessage("You have requested to convert a plugin of type <"
+                             + boost::core::demangle(typeid(PluginType).name())
+                             + "> into type <"
+                             + boost::core::demangle(typeid(TestType).name()) +
+                             ">, but this cast cannot be performed."));
+
+      // We can safely dereference the pointer, because we checked above that
+      // the object is actually of type TestType, and so the result
+      // is not a nullptr.
+      return *dynamic_cast<TestType *> (&object);
+    }
+  }
 
   namespace internal
   {

include/aspect/postprocess/depth_average.h

@@ -116,12 +116,7 @@ namespace aspect
         /**
          * List of the quantities to calculate for each depth zone.
          */
-        std::vector<std::string> output_variables;
-
-        /**
-         * Whether to calculate all available quantities when averaging.
-         */
-        bool output_all_variables;
+        std::vector<std::string> variables;
 
         /**
          * Whether to use plain ascii text output