crystal_preferred_orientation.h

Go to the documentation of this file.

/*
 Copyright (C) 2022 - 2024 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_particle_property_cpo_h
#define _aspect_particle_property_cpo_h

#include <aspect/particle/property/interface.h>
#include <aspect/simulator_access.h>
#include <array>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#include <boost/random.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

namespace aspect
{
  namespace Particle
  {
    namespace Property
    {
      enum class DeformationType
      {
        passive, olivine_a_fabric, olivine_b_fabric, olivine_c_fabric, olivine_d_fabric, olivine_e_fabric, enstatite
      };


      enum class DeformationTypeSelector
      {
        passive, olivine_a_fabric, olivine_b_fabric, olivine_c_fabric, olivine_d_fabric, olivine_e_fabric, enstatite, olivine_karato_2008
      };

      enum class AdvectionMethod
      {
        forward_euler, backward_euler
      };

      enum class CPODerivativeAlgorithm
      {
        spin_tensor, drex_2004
      };

      enum class CPOInitialGrainsModel
      {
        uniform_grains_and_random_uniform_rotations, world_builder
      };

      template <int dim>
      class CrystalPreferredOrientation : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
      {
        public:
          CrystalPreferredOrientation() = default;

          void
          initialize () override;

          void
          initialize_one_particle_property (const Point<dim> &position,
                                            std::vector<double> &particle_properties) const override;

          void
          update_particle_properties (const ParticleUpdateInputs<dim> &inputs,
                                      typename ParticleHandler<dim>::particle_iterator_range &particles) const override;

          UpdateTimeFlags
          need_update () const override;

          InitializationModeForLateParticles
          late_initialization_mode () const override;

          UpdateFlags
          get_update_flags (const unsigned int component) const override;

          std::vector<std::pair<std::string, unsigned int>>
          get_property_information() const override;

          std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
          compute_derivatives(const unsigned int cpo_index,
                              const ArrayView<double> &data,
                              const unsigned int mineral_i,
                              const SymmetricTensor<2,3> &strain_rate_3d,
                              const Tensor<2,3> &velocity_gradient_tensor,
                              const Point<dim> &position,
                              const double temperature,
                              const double pressure,
                              const Tensor<1,dim> &velocity,
                              const std::vector<double> &compositions,
                              const SymmetricTensor<2,dim> &strain_rate,
                              const SymmetricTensor<2,dim> &deviatoric_strain_rate,
                              const double water_content) const;

          std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
          compute_derivatives_drex_2004(const unsigned int cpo_index,
                                        const ArrayView<double> &data,
                                        const unsigned int mineral_i,
                                        const SymmetricTensor<2,3> &strain_rate_3d,
                                        const Tensor<2,3> &velocity_gradient_tensor,
                                        const std::array<double,4> ref_resolved_shear_stress,
                                        const bool prevent_nondimensionalization = false) const;


          static
          void
          declare_parameters (ParameterHandler &prm);

          void
          parse_parameters (ParameterHandler &prm) override;

          unsigned int
          get_number_of_grains() const;

          unsigned int
          get_number_of_minerals() const;

          DeformationType
          determine_deformation_type(const DeformationTypeSelector deformation_type_selector,
                                     const Point<dim> &position,
                                     const double temperature,
                                     const double pressure,
                                     const Tensor<1,dim> &velocity,
                                     const std::vector<double> &compositions,
                                     const SymmetricTensor<2,dim> &strain_rate,
                                     const SymmetricTensor<2,dim> &deviatoric_strain_rate,
                                     const double water_content) const;

          DeformationType
          determine_deformation_type_karato_2008(const double stress,
                                                 const double water_content) const;

          std::array<double,4>
          reference_resolved_shear_stress_from_deformation_type(DeformationType deformation_type,
                                                                double max_value = 1e60) const;

          inline
          DeformationType get_deformation_type(const unsigned int cpo_data_position,
                                               const ArrayView<double> &data,
                                               const unsigned int mineral_i) const
          {
            return static_cast<DeformationType>(data[cpo_data_position + 0 + mineral_i * (n_grains * 10 + 2)]);
          }

          inline
          void set_deformation_type(const unsigned int cpo_data_position,
                                    const ArrayView<double> &data,
                                    const unsigned int mineral_i,
                                    const DeformationType deformation_type) const
          {
            data[cpo_data_position + 0 + mineral_i * (n_grains * 10 + 2)] = static_cast<double>(deformation_type);
          }

          inline
          double get_volume_fraction_mineral(const unsigned int cpo_data_position,
                                             const ArrayView<double> &data,
                                             const unsigned int mineral_i) const
          {
            return data[cpo_data_position + 1 + mineral_i *(n_grains * 10 + 2)];
          }

          inline
          void set_volume_fraction_mineral(const unsigned int cpo_data_position,
                                           const ArrayView<double> &data,
                                           const unsigned int mineral_i,
                                           const double volume_fraction_mineral) const
          {
            data[cpo_data_position + 1 + mineral_i *(n_grains * 10 + 2)] = volume_fraction_mineral;
          }

          inline
          double get_volume_fractions_grains(const unsigned int cpo_data_position,
                                             const ArrayView<const double> &data,
                                             const unsigned int mineral_i,
                                             const unsigned int grain_i) const
          {
            return data[cpo_data_position + 2 + grain_i * 10 + mineral_i * (n_grains * 10 + 2)];
          }

          inline
          void set_volume_fractions_grains(const unsigned int cpo_data_position,
                                           const ArrayView<double> &data,
                                           const unsigned int mineral_i,
                                           const unsigned int grain_i,
                                           const double volume_fractions_grains) const
          {
            data[cpo_data_position + 2 + grain_i * 10 + mineral_i * (n_grains * 10 + 2)] = volume_fractions_grains;
          }

          inline
          Tensor<2,3> get_rotation_matrix_grains(const unsigned int cpo_data_position,
                                                 const ArrayView<const double> &data,
                                                 const unsigned int mineral_i,
                                                 const unsigned int grain_i) const
          {
            Tensor<2,3> rotation_matrix;
            for (unsigned int i = 0; i < Tensor<2,3>::n_independent_components ; ++i)
              {
                const ::TableIndices<2> index = Tensor<2,3>::unrolled_to_component_indices(i);
                rotation_matrix[index] = data[cpo_data_position + 3 + grain_i * 10 + mineral_i * (n_grains * 10 + 2) + i];
              }
            return rotation_matrix;
          }

          inline
          void set_rotation_matrix_grains(const unsigned int cpo_data_position,
                                          const ArrayView<double> &data,
                                          const unsigned int mineral_i,
                                          const unsigned int grain_i,
                                          const Tensor<2,3> &rotation_matrix) const
          {
            for (unsigned int i = 0; i < Tensor<2,3>::n_independent_components ; ++i)
              {
                const ::TableIndices<2> index = Tensor<2,3>::unrolled_to_component_indices(i);
                data[cpo_data_position + 3 + grain_i * 10 + mineral_i * (n_grains * 10 + 2) + i] = rotation_matrix[index];
              }
          }


        private:
          void
          compute_random_rotation_matrix(Tensor<2,3> &rotation_matrix) const;

          double
          advect_forward_euler(const unsigned int cpo_data_position,
                               const ArrayView<double> &data,
                               const unsigned int mineral_i,
                               const double dt,
                               const std::pair<std::vector<double>, std::vector<Tensor<2,3>>> &derivatives) const;

          double
          advect_backward_euler(const unsigned int cpo_data_position,
                                const ArrayView<double> &data,
                                const unsigned int mineral_i,
                                const double dt,
                                const std::pair<std::vector<double>, std::vector<Tensor<2,3>>> &derivatives) const;


          std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
          compute_derivatives_spin_tensor(const Tensor<2,3> &velocity_gradient_tensor) const;

          mutable boost::mt19937 random_number_generator;
          unsigned int random_number_seed;

          unsigned int n_grains;

          unsigned int n_minerals;

          unsigned int water_index;

          std::vector<DeformationTypeSelector> deformation_type_selector;

          std::vector<double> volume_fractions_minerals;

          AdvectionMethod advection_method;

          CPODerivativeAlgorithm cpo_derivative_algorithm;

          double property_advection_tolerance;

          unsigned int property_advection_max_iterations;

          double stress_exponent;

          double nucleation_efficiency;

          double exponent_p;

          double threshold_GBS;

          double mobility;

          CPOInitialGrainsModel initial_grains_model;

      };
    }
  }
}

#endif