fastscape.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 2026 by the authors of the ASPECT code.

  This file is part of ASPECT.

  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.

  ASPECT is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
 */

#ifndef _aspect_mesh_deformation_fastscape_h
#define _aspect_mesh_deformation_fastscape_h

#include <aspect/global.h>

#include <aspect/mesh_deformation/interface.h>
#include <deal.II/base/parsed_function.h>

namespace aspect
{
  namespace MeshDeformation
  {
    template <int dim>
    class FastScape : public Interface<dim>, public SimulatorAccess<dim>
    {
      public:
        ~FastScape() override;

        virtual void initialize () override;

        void update() override;

        virtual
        void
        compute_velocity_constraints_on_boundary(const DoFHandler<dim> &mesh_deformation_dof_handler,
                                                 AffineConstraints<double> &mesh_velocity_constraints,
                                                 const std::set<types::boundary_id> &boundary_ids) const override;

        double boundary_composition (const types::boundary_id boundary_indicator,
                                     const Point<dim> &position,
                                     const unsigned int compositional_field) const override;

        bool needs_surface_stabilization () const override;

        static
        void declare_parameters (ParameterHandler &prm);

        void parse_parameters (ParameterHandler &prm) override;

        enum class FastscapeOutputVariable
        {
          //
          kf,
          kd,
          uplift_rate
        };
        FastscapeOutputVariable additional_output_variable;

        template <class Archive>
        void serialize (Archive &ar, const unsigned int version);

        void save (std::map<std::string, std::string> &status_strings) const override;

        void load (const std::map<std::string, std::string> &status_strings) override;

      private:
        void set_ghost_nodes(std::vector<double> &elevation,
                             std::vector<double> &velocity_x,
                             std::vector<double> &velocity_y,
                             std::vector<double> &velocity_z,
                             std::vector<double> &bedrock_transport_coefficient_array,
                             const double &fastscape_timestep_in_years,
                             const bool init) const;

        bool is_ghost_node(const unsigned int &index,
                           const bool &exclude_boundaries) const;

        void fill_fastscape_arrays(std::vector<double> &elevation,
                                   std::vector<double> &bedrock_transport_coefficient_array,
                                   std::vector<double> &bedrock_river_incision_rate_array,
                                   std::vector<double> &velocity_x,
                                   std::vector<double> &velocity_y,
                                   std::vector<double> &velocity_z,
                                   std::vector<std::vector<double>> &temporary_variables) const;

        std::vector<std::vector<double>> get_aspect_values() const;

        void initialize_fastscape(std::vector<double> &elevation,
                                  std::vector<double> &basement,
                                  std::vector<double> &silt_fraction,
                                  bool restart) const;

        void execute_fastscape(std::vector<double> &elevation,
                               std::vector<double> &extra_vtk_field,
                               std::vector<double> &velocity_x,
                               std::vector<double> &velocity_y,
                               std::vector<double> &velocity_z,
                               std::vector<double> &bedrock_transport_coefficient_array,
                               const double &fastscape_timestep_in_years,
                               const unsigned int &fastscape_iterations) const;

        void apply_orographic_controls(const std::vector<double> &elevation,
                                       std::vector<double> &bedrock_river_incision_rate_array,
                                       std::vector<double> &bedrock_transport_coefficient_array) const;

        Table<dim,double> fill_data_table(const std::vector<double> &values,
                                          const TableIndices<dim> &size_idx,
                                          const unsigned int &fastscape_nx,
                                          const unsigned int &fastscape_ny) const;

        unsigned int fastscape_steps_per_aspect_step;

        double maximum_fastscape_timestep;

        double fastscape_dx;

        double fastscape_dy;

        double fastscape_x_extent;

        double fastscape_y_extent;

        double fastscape_y_extent_2d;

        unsigned int fastscape_nx;

        unsigned int fastscape_ny;

        double vexp;

        unsigned int additional_refinement_levels;

        unsigned int maximum_surface_refinement_level;

        unsigned int surface_refinement_difference;

        bool average_out_of_plane_surface_topography;

        int fastscape_seed;

        std::array<std::pair<double,double>,dim> grid_extent;

        std::array<unsigned int, dim> table_intervals;

        bool use_ghost_nodes;

        double noise_elevation;

        std::vector<double> sediment_rain_rates;

        std::vector<double> sediment_rain_times;

        bool fastscape_advection_uplift;

        double node_tolerance;

        double output_interval;

        mutable double last_output_time;

        unsigned int bottom;

        unsigned int top;

        unsigned int right;

        unsigned int left;

        bool topbottom_ghost_nodes_periodic;

        bool leftright_ghost_nodes_periodic;

        unsigned int fastscape_boundary_conditions;

        double bottom_flux;

        double top_flux;

        double right_flux;

        double left_flux;
        double drainage_area_exponent_m;

        double slope_exponent_n;

        double slope_exponent_p;

        double bedrock_deposition_g;

        double sediment_deposition_g;

        Functions::ParsedFunction<2> kf_distribution_function;

        bool use_kf_distribution_function;

        double constant_bedrock_river_incision_rate;

        double sediment_river_incision_rate;

        Functions::ParsedFunction<2> kd_distribution_function;

        bool use_kd_distribution_function;

        double constant_bedrock_transport_coefficient;

        double sediment_transport_coefficient;
        double sea_level_constant_value;

        Functions::ParsedFunction<1> sea_level_function;

        bool use_sea_level_function;

        bool use_fixed_erosional_base;

        double h_erosional_base;

        double sand_surface_porosity;

        double silt_surface_porosity;

        double sand_efold_depth;

        double silt_efold_depth;

        double incoming_silt_fraction;

        double sand_silt_averaging_depth;

        double sand_transport_coefficient;

        double silt_transport_coefficient;

        bool use_marine_component;
        int flat_elevation;

        int wind_barrier_elevation;

        unsigned int wind_direction;

        double flat_erosional_factor;

        double wind_barrier_erosional_factor;

        bool stack_controls;

        bool use_orographic_controls;
    };
  }
}

#endif