visco_plastic.h

Go to the documentation of this file.

/*
  Copyright (C) 2020 - 2023 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_material_model_rheology_visco_plastic_h
#define _aspect_material_model_rheology_visco_plastic_h

#include <aspect/global.h>
#include <aspect/material_model/interface.h>
#include <aspect/material_model/utilities.h>
#include <aspect/material_model/rheology/strain_dependent.h>
#include <aspect/material_model/rheology/friction_models.h>
#include <aspect/material_model/rheology/diffusion_creep.h>
#include <aspect/material_model/rheology/dislocation_creep.h>
#include <aspect/material_model/rheology/frank_kamenetskii.h>
#include <aspect/material_model/rheology/peierls_creep.h>
#include <aspect/material_model/rheology/constant_viscosity_prefactors.h>
#include <aspect/material_model/rheology/drucker_prager.h>
#include <aspect/material_model/rheology/elasticity.h>
#include <aspect/simulator_access.h>

#include<deal.II/fe/component_mask.h>

namespace aspect
{
  namespace MaterialModel
  {
    using namespace dealii;

    template <int dim>
    class PlasticAdditionalOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PlasticAdditionalOutputs(const unsigned int n_points);

        std::vector<double> get_nth_output(const unsigned int idx) const override;

        std::vector<double> cohesions;

        std::vector<double> friction_angles;

        std::vector<double> yield_stresses;

        std::vector<double> yielding;

    };

    struct IsostrainViscosities
    {
      std::vector<double> composition_viscosities;

      std::vector<bool> composition_yielding;

      std::vector<double> current_friction_angles;

      std::vector<double> current_cohesions;
    };

    namespace Rheology
    {

      template <int dim>
      class ViscoPlastic : public ::aspect::SimulatorAccess<dim>
      {
        public:
          ViscoPlastic();

          IsostrainViscosities
          calculate_isostrain_viscosities ( const MaterialModel::MaterialModelInputs<dim> &in,
                                            const unsigned int i,
                                            const std::vector<double> &volume_fractions,
                                            const std::vector<double> &phase_function_values = std::vector<double>(),
                                            const std::vector<unsigned int> &n_phase_transitions_per_composition =
                                              std::vector<unsigned int>()) const;

          void compute_viscosity_derivatives(const unsigned int point_index,
                                             const std::vector<double> &volume_fractions,
                                             const std::vector<double> &composition_viscosities,
                                             const MaterialModel::MaterialModelInputs<dim> &in,
                                             MaterialModel::MaterialModelOutputs<dim> &out,
                                             const std::vector<double> &phase_function_values = std::vector<double>(),
                                             const std::vector<unsigned int> &n_phase_transitions_per_composition =
                                               std::vector<unsigned int>()) const;

          ComponentMask get_volumetric_composition_mask() const;

          static
          void
          declare_parameters (ParameterHandler &prm);

          void
          parse_parameters (ParameterHandler &prm,
                            const std::unique_ptr<std::vector<unsigned int>> &expected_n_phases_per_composition = nullptr);

          void
          create_plastic_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const;

          void fill_plastic_outputs(const unsigned int point_index,
                                    const std::vector<double> &volume_fractions,
                                    const bool plastic_yielding,
                                    const MaterialModel::MaterialModelInputs<dim> &in,
                                    MaterialModel::MaterialModelOutputs<dim> &out,
                                    const IsostrainViscosities &isostrain_viscosities) const;

          double min_strain_rate;

          MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;

          Rheology::StrainDependent<dim> strain_rheology;

          Rheology::FrictionModels<dim> friction_models;

          Rheology::Elasticity<dim> elastic_rheology;


        private:

          double ref_strain_rate;

          std::vector<double> minimum_viscosity;
          std::vector<double> maximum_viscosity;

          enum ViscosityScheme
          {
            diffusion,
            dislocation,
            frank_kamenetskii,
            composite
          } viscous_flow_law;

          enum YieldScheme
          {
            stress_limiter,
            drucker_prager
          } yield_mechanism;

          bool allow_negative_pressures_in_plasticity;

          bool use_adiabatic_pressure_in_creep;

          std::vector<double> exponents_stress_limiter;

          double adiabatic_temperature_gradient_for_viscosity;

          Rheology::DiffusionCreep<dim> diffusion_creep;
          Rheology::DislocationCreep<dim> dislocation_creep;
          std::unique_ptr<Rheology::FrankKamenetskii<dim>> frank_kamenetskii_rheology;

          bool use_peierls_creep;

          std::unique_ptr<Rheology::PeierlsCreep<dim>> peierls_creep;

          Rheology::ConstantViscosityPrefactors<dim> constant_viscosity_prefactors;

          /*
           * Object for computing plastic stresses, viscosities, and additional outputs
           */
          Rheology::DruckerPrager<dim> drucker_prager_plasticity;

          /*
           * Input parameters for the drucker prager plasticity.
           */
          Rheology::DruckerPragerParameters drucker_prager_parameters;

      };
    }
  }
}
#endif