ASPECT
crystal_preferred_orientation.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2022 - 2023 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19  */
20 
21 #ifndef _aspect_particle_property_cpo_h
22 #define _aspect_particle_property_cpo_h
23 
26 #include <array>
27 
28 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
29 #include <boost/random.hpp>
30 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
31 
32 namespace aspect
33 {
34  namespace Particle
35  {
36  namespace Property
37  {
46  enum class DeformationType
47  {
49  };
50 
51 
66  {
68  };
69 
73  enum class AdvectionMethod
74  {
76  };
77 
85  {
87  };
88 
96  {
98  };
99 
129  template <int dim>
131  {
132  public:
136  CrystalPreferredOrientation() = default;
137 
142  void
143  initialize () override;
144 
156  void
157  initialize_one_particle_property (const Point<dim> &position,
158  std::vector<double> &particle_properties) const override;
159 
186  void
187  update_particle_property (const unsigned int data_position,
188  const Vector<double> &solution,
189  const std::vector<Tensor<1,dim>> &gradients,
190  typename ParticleHandler<dim>::particle_iterator &particle) const override;
191 
197  need_update () const override;
198 
203  late_initialization_mode () const override;
204 
209  UpdateFlags
210  get_needed_update_flags () const override;
211 
219  std::vector<std::pair<std::string, unsigned int>>
220  get_property_information() const override;
221 
239  std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
240  compute_derivatives(const unsigned int cpo_index,
241  const ArrayView<double> &data,
242  const unsigned int mineral_i,
243  const SymmetricTensor<2,3> &strain_rate_3d,
244  const Tensor<2,3> &velocity_gradient_tensor,
245  const Point<dim> &position,
246  const double temperature,
247  const double pressure,
248  const Tensor<1,dim> &velocity,
249  const std::vector<double> &compositions,
250  const SymmetricTensor<2,dim> &strain_rate,
251  const SymmetricTensor<2,dim> &deviatoric_strain_rate,
252  const double water_content) const;
253 
269  std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
270  compute_derivatives_drex_2004(const unsigned int cpo_index,
271  const ArrayView<double> &data,
272  const unsigned int mineral_i,
273  const SymmetricTensor<2,3> &strain_rate_3d,
274  const Tensor<2,3> &velocity_gradient_tensor,
275  const std::array<double,4> ref_resolved_shear_stress,
276  const bool prevent_nondimensionalization = false) const;
277 
278 
282  static
283  void
284  declare_parameters (ParameterHandler &prm);
285 
289  void
290  parse_parameters (ParameterHandler &prm) override;
291 
295  unsigned int
296  get_number_of_grains() const;
297 
301  unsigned int
302  get_number_of_minerals() const;
303 
314  determine_deformation_type(const DeformationTypeSelector deformation_type_selector,
315  const Point<dim> &position,
316  const double temperature,
317  const double pressure,
318  const Tensor<1,dim> &velocity,
319  const std::vector<double> &compositions,
320  const SymmetricTensor<2,dim> &strain_rate,
321  const SymmetricTensor<2,dim> &deviatoric_strain_rate,
322  const double water_content) const;
323 
329  determine_deformation_type_karato_2008(const double stress,
330  const double water_content) const;
331 
341  std::array<double,4>
342  reference_resolved_shear_stress_from_deformation_type(DeformationType deformation_type,
343  double max_value = 1e60) const;
344 
352  inline
353  DeformationType get_deformation_type(const unsigned int cpo_data_position,
354  const ArrayView<double> &data,
355  const unsigned int mineral_i) const
356  {
357  return static_cast<DeformationType>(data[cpo_data_position + 0 + mineral_i * (n_grains * 10 + 2)]);
358  }
359 
368  inline
369  void set_deformation_type(const unsigned int cpo_data_position,
370  const ArrayView<double> &data,
371  const unsigned int mineral_i,
372  const DeformationType deformation_type) const
373  {
374  data[cpo_data_position + 0 + mineral_i * (n_grains * 10 + 2)] = static_cast<double>(deformation_type);
375  }
376 
384  inline
385  double get_volume_fraction_mineral(const unsigned int cpo_data_position,
386  const ArrayView<double> &data,
387  const unsigned int mineral_i) const
388  {
389  return data[cpo_data_position + 1 + mineral_i *(n_grains * 10 + 2)];
390  }
391 
400  inline
401  void set_volume_fraction_mineral(const unsigned int cpo_data_position,
402  const ArrayView<double> &data,
403  const unsigned int mineral_i,
404  const double volume_fraction_mineral) const
405  {
406  data[cpo_data_position + 1 + mineral_i *(n_grains * 10 + 2)] = volume_fraction_mineral;
407  }
408 
417  inline
418  double get_volume_fractions_grains(const unsigned int cpo_data_position,
419  const ArrayView<const double> &data,
420  const unsigned int mineral_i,
421  const unsigned int grain_i) const
422  {
423  return data[cpo_data_position + 2 + grain_i * 10 + mineral_i * (n_grains * 10 + 2)];
424  }
425 
435  inline
436  void set_volume_fractions_grains(const unsigned int cpo_data_position,
437  const ArrayView<double> &data,
438  const unsigned int mineral_i,
439  const unsigned int grain_i,
440  const double volume_fractions_grains) const
441  {
442  data[cpo_data_position + 2 + grain_i * 10 + mineral_i * (n_grains * 10 + 2)] = volume_fractions_grains;
443  }
444 
454  inline
455  Tensor<2,3> get_rotation_matrix_grains(const unsigned int cpo_data_position,
456  const ArrayView<const double> &data,
457  const unsigned int mineral_i,
458  const unsigned int grain_i) const
459  {
460  Tensor<2,3> rotation_matrix;
461  for (unsigned int i = 0; i < Tensor<2,3>::n_independent_components ; ++i)
462  {
463  const ::TableIndices<2> index = Tensor<2,3>::unrolled_to_component_indices(i);
464  rotation_matrix[index] = data[cpo_data_position + 3 + grain_i * 10 + mineral_i * (n_grains * 10 + 2) + i];
465  }
466  return rotation_matrix;
467  }
468 
478  inline
479  void set_rotation_matrix_grains(const unsigned int cpo_data_position,
480  const ArrayView<double> &data,
481  const unsigned int mineral_i,
482  const unsigned int grain_i,
483  const Tensor<2,3> &rotation_matrix) const
484  {
485  for (unsigned int i = 0; i < Tensor<2,3>::n_independent_components ; ++i)
486  {
487  const ::TableIndices<2> index = Tensor<2,3>::unrolled_to_component_indices(i);
488  data[cpo_data_position + 3 + grain_i * 10 + mineral_i * (n_grains * 10 + 2) + i] = rotation_matrix[index];
489  }
490  }
491 
492 
493  private:
497  void
498  compute_random_rotation_matrix(Tensor<2,3> &rotation_matrix) const;
499 
515  double
516  advect_forward_euler(const unsigned int cpo_data_position,
517  const ArrayView<double> &data,
518  const unsigned int mineral_i,
519  const double dt,
520  const std::pair<std::vector<double>, std::vector<Tensor<2,3>>> &derivatives) const;
521 
537  double
538  advect_backward_euler(const unsigned int cpo_data_position,
539  const ArrayView<double> &data,
540  const unsigned int mineral_i,
541  const double dt,
542  const std::pair<std::vector<double>, std::vector<Tensor<2,3>>> &derivatives) const;
543 
544 
551  std::pair<std::vector<double>, std::vector<Tensor<2,3>>>
552  compute_derivatives_spin_tensor(const Tensor<2,3> &velocity_gradient_tensor) const;
553 
557  mutable boost::mt19937 random_number_generator;
558  unsigned int random_number_seed;
559 
560  unsigned int n_grains;
561 
562  unsigned int n_minerals;
563 
567  unsigned int water_index;
568 
575  std::vector<DeformationTypeSelector> deformation_type_selector;
576 
580  std::vector<double> volume_fractions_minerals;
581 
586 
591 
597 
603 
612 
618 
622  double exponent_p;
623 
630 
635  double mobility;
636 
641 
644  };
645  }
646  }
647 }
648 
649 #endif
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
Sets the value in the data array representing the volume fraction of a mineral.
DeformationTypeSelector
The type of deformation selector used by the CPO code.
CPODerivativeAlgorithm
The algorithm used to compute the derivatives of the grain size and rotation matrix used in the advec...
DeformationType
The type of deformation used by the CPO code.
double get_volume_fraction_mineral(const unsigned int cpo_data_position, const ArrayView< double > &data, const unsigned int mineral_i) const
Returns the value in the data array representing the volume fraction of a mineral.
void declare_parameters(ParameterHandler &prm)
AdvectionMethod
The type of Advection method used to advect the CPO properties.
CPOInitialGrainsModel
An enum used to determine how the initial grain sizes and orientations are set for all particles...
Definition: compat.h:59
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
Gets the rotation matrix for a grain in a mineral.
void set_deformation_type(const unsigned int cpo_data_position, const ArrayView< double > &data, const unsigned int mineral_i, const DeformationType deformation_type) const
Sets the value in the data array representing the deformation type.
DeformationType get_deformation_type(const unsigned int cpo_data_position, const ArrayView< double > &data, const unsigned int mineral_i) const
Returns the value in the data array representing the deformation type.
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
Sets the rotation matrix for a grain in a mineral.
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
Returns the value in the data array representing the volume fraction of a grain.
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
Sets the value in the data array representing the volume fraction of a grain.