Isostatic initial topography

include/aspect/geometry_model/initial_topography_model/interface.h

@@ -51,6 +51,39 @@ namespace aspect
     class Interface : public Plugins::InterfaceBase
     {
       public:
+        /**
+        * Constructor.
+        * By default, the plugin is considered fully initialized and
+        * does not require an additional initialization step.
+        */
+        Interface ()
+          : required_initialized(true)
+        {}
+
+        /**
+        * Perform an additional initialization step after the material model,
+        * initial temperature model, initial composition model, gravity model,
+        * and adiabatic conditions have been initialized.
+        *
+        * Plugins may override this function if they need to query information
+        * from these modules before they are fully initialized.
+        */
+        virtual
+        void
+        required_initialize()
+        {}
+
+
+        /**
+         * Return whether this plugin has been initialized.
+         */
+        bool
+        is_required_initialized () const
+        {
+          return required_initialized;
+        }
+
+
         /**
          * Return the value of the elevation at the given surface point.
          *
@@ -69,6 +102,23 @@ namespace aspect
          */
         virtual
         double max_topography () const = 0;
+
+      protected:
+        /**
+         * Set the initialization state.
+         */
+        void
+        set_required_initialized (const bool state)
+        {
+          required_initialized = state;
+        }
+
+      private:
+        /**
+        * Whether the plugin has completed its required initialization.
+        */
+        bool required_initialized;
+
     };
 
 

include/aspect/geometry_model/initial_topography_model/isostatic_topography.h

@@ -0,0 +1,108 @@
+#ifndef _aspect_geometry_model_initial_topography_model_isostatic_topography_h
+#define _aspect_geometry_model_initial_topography_model_isostatic_topography_h
+
+#include <aspect/geometry_model/initial_topography_model/interface.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace InitialTopographyModel
+  {
+    /**
+    * An initial topography model that computes isostatic topography
+    * from the initial temperature and compositional fields by
+    * balancing the mass of vertical columns.
+    */
+    template <int dim>
+    class IsostaticTopography : public SimulatorAccess<dim>, public Interface<dim>
+    {
+      public:
+        /**
+        * Constructor.
+        */
+        IsostaticTopography();
+
+        /**
+        * Perform basic initialization and verify that the selected
+        * geometry model is supported.
+        */
+        void initialize() override;
+
+        /**
+        * Perform an additional initialization step after the material model,
+        * initial temperature model, initial composition model, gravity model,
+        * and adiabatic conditions have been initialized.
+        */
+        void required_initialize() override;
+
+        /**
+         * Return the initial topography at the given surface point.
+         */
+        double
+        value (const Point<dim-1> &surface_point) const override;
+
+
+        /**
+         * Return the maximum value of the initial topography.
+         */
+        double
+        max_topography () const override;
+
+        /**
+         * Declare the parameters this class takes.
+         */
+        static
+        void
+        declare_parameters (ParameterHandler &prm);
+
+        /**
+         * Read the parameters this class declares.
+         */
+        void
+        parse_parameters (ParameterHandler &prm) override;
+
+
+      private:
+        /**
+         * Number of equally spaced lateral sampling points used to
+         * construct the isostatic topography profile.
+         */
+        double n_lateral_points;
+
+        /**
+         * Number of equally spaced vertical sampling points used to
+         * integrate the mass of each column.
+         */
+        unsigned n_vertical_points;
+
+        /**
+         * The depth above which the model is assumed to be in
+         * isostatic equilibrium, and where the reference density is evaluated
+         */
+        double compensation_depth;
+
+        /**
+         * Maximum absolute value of the computed isostatic topography.
+         */
+        double max_isostatic_topography;
+
+        /**
+         * The characteristic horizontal length scale (in m) over which
+         * topographic loads are assumed to attain isostatic compensation.
+         */
+        double isostatic_length_scale;
+
+        /**
+         * Maximum computed initial topography.
+         */
+        double maximal_topography;
+
+        /**
+         * Isostatic topography evaluated at the lateral sampling points.
+         */
+        std::vector<double> topography;
+    };
+  }
+}
+
+#endif

source/adiabatic_conditions/compute_profile.cc

@@ -293,7 +293,9 @@ namespace aspect
       if (normalized_distance_to_closest_profile_point >=0.0 && normalized_distance_to_closest_profile_point < 1e-6)
         return property[i];
 
-      Assert (i+1 < property.size(), ExcInternalError());
+      // Assert (i+1 < property.size(), ExcInternalError());
+      if (i+1 >= property.size())
+        return property[property.size()-1];
 
       // now do the linear interpolation
       const double d = normalized_distance_from_surface - i;