composite_visco_plastic.h

Go to the documentation of this file.

// ------------------------------------------------------
//
// SPDX-License-Identifier: GPL-2.0-or-later
//
// SPDX-FileCopyrightText: Copyright (C) 2020-2026 by the ASPECT authors.
//
// This file is part of ASPECT.
//
// Detailed license information governing the source code and
// contributions can be found in the folder LICENSES and in
// LICENSE and CONTRIBUTING.md at the top level directory.
//
// ------------------------------------------------------

#ifndef _aspect_material_model_rheology_composite_visco_plastic_h
#define _aspect_material_model_rheology_composite_visco_plastic_h

#include <aspect/global.h>
#include <aspect/material_model/interface.h>
#include <aspect/material_model/utilities.h>
#include <aspect/material_model/rheology/diffusion_creep.h>
#include <aspect/material_model/rheology/dislocation_creep.h>
#include <aspect/material_model/rheology/peierls_creep.h>
#include <aspect/material_model/rheology/drucker_prager_power.h>
#include <aspect/simulator_access.h>

namespace aspect
{
  namespace MaterialModel
  {
    namespace Rheology
    {
      namespace ViscosityAveraging
      {
        enum Kind
        {
          isostress,
          isostrain
        };

        ViscosityAveraging::Kind
        parse (const std::string &parameter_name,
               const ParameterHandler &prm);
      }

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

          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);

          double
          compute_viscosity (const double pressure,
                             const double temperature,
                             const double grain_size,
                             const std::vector<double> &volume_fractions,
                             const SymmetricTensor<2,dim> &strain_rate,
                             std::vector<double> &partial_strain_rates,
                             const std::vector<double> &phase_function_values = {},
                             const std::vector<unsigned int> &n_phase_transitions_per_composition = {}) const;

        private:
          double
          compute_isostress_viscosity (const double pressure,
                                       const double temperature,
                                       const double grain_size,
                                       const std::vector<double> &volume_fractions,
                                       const SymmetricTensor<2,dim> &strain_rate,
                                       std::vector<double> &partial_strain_rates,
                                       const std::vector<double> &phase_function_values = {},
                                       const std::vector<unsigned int> &n_phase_transitions_per_composition = {}) const;

          std::pair<double, double>
          calculate_isostress_log_strain_rate_and_derivative(const std::vector<std::array<std::pair<double, double>, 4>> &logarithmic_strain_rates_and_stress_derivatives,
                                                             const double viscoplastic_stress,
                                                             std::vector<double> &partial_strain_rates) const;

          double
          compute_isostrain_viscosity (const double pressure,
                                       const double temperature,
                                       const double grain_size,
                                       const std::vector<double> &volume_fractions,
                                       const SymmetricTensor<2,dim> &strain_rate,
                                       std::vector<double> &partial_strain_rates,
                                       const std::vector<double> &phase_function_values = {},
                                       const std::vector<unsigned int> &n_phase_transitions_per_composition = {}) const;


          double
          compute_composition_viscosity (const double pressure,
                                         const double temperature,
                                         const double grain_size,
                                         const unsigned int composition,
                                         const SymmetricTensor<2,dim> &strain_rate,
                                         std::vector<double> &partial_strain_rates,
                                         const std::vector<double> &phase_function_values = {},
                                         const std::vector<unsigned int> &n_phase_transitions_per_composition = {}) const;

          std::pair<double, double>
          calculate_composition_log_strain_rate_and_derivative(const std::array<std::pair<double, double>, 4> &logarithmic_strain_rates_and_stress_derivatives,
                                                               const double viscoplastic_stress,
                                                               std::vector<double> &partial_strain_rates) const;

          ViscosityAveraging::Kind viscosity_averaging_scheme;

          bool use_diffusion_creep;
          bool use_dislocation_creep;
          bool use_peierls_creep;
          bool use_drucker_prager;

          std::vector<unsigned int> active_flow_mechanisms;

          static constexpr unsigned int n_decomposed_strain_rates = 5;

          std::unique_ptr<Rheology::DiffusionCreep<dim>> diffusion_creep;
          std::unique_ptr<Rheology::DislocationCreep<dim>> dislocation_creep;
          std::unique_ptr<Rheology::PeierlsCreep<dim>> peierls_creep;
          std::unique_ptr<Rheology::DruckerPragerPower<dim>> drucker_prager;

          unsigned int number_of_chemical_compositions;

          double maximum_viscosity;


          double damper_viscosity;

          double strain_rate_scaling_factor;

          double minimum_strain_rate;

          double log_strain_rate_residual_threshold;

          unsigned int stress_max_iteration_number;

          const double logmin = std::log(std::numeric_limits<double>::min());
      };
    }
  }
}
#endif