Clean up and rearrange

Author: Wang-yijun

cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md

@@ -1,3 +1,11 @@
+```{tags}
+category:cookbook
+feature:3d
+feature:cartesian
+feature:modular-equations
+feature:compositional-fields
+feature:particles
+
 (sec:cookbooks:CPO_induced_anisotropic_viscosity)=
 # CPO induced anisotropic viscosity
 

cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc

@@ -113,17 +113,14 @@ namespace aspect
       ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(AnisotropicStress,
                                                   "Anisotropic stress",
                                                   "A visualization output object that generates output "
-                                                  "for the 6 (in 3d) components of the shear stress "
-                                                  "tensor, i.e., for the components of the tensor "
-                                                  "$-2\\eta\\varepsilon(\\mathbf u)$ "
-                                                  "in the incompressible case and "
-                                                  "$-2\\eta\\left[\\varepsilon(\\mathbf u)-"
-                                                  "\\tfrac 13(\\textrm{tr}\\;\\varepsilon(\\mathbf u))\\mathbf I\\right]$ "
-                                                  "in the compressible case. If elasticity is used, the "
-                                                  "elastic contribution is being accounted for. The shear "
-                                                  "stress differs from the full stress tensor "
-                                                  "by the absence of the pressure. Note that the convention "
-                                                  "of positive compressive stress is followed. ")
+                                                  "for the 6 (in 3d) components of the anisotropic stress "
+                                                  "tensor. The anisotropic stress is defined as $2 \eta "
+                                                  "(\varepsilon(\mathbf u) - \tfrac 13 \textrm{trace}\ "
+                                                  "\varepsilon(\mathbf u) \mathbf 1) = 2\eta (\varepsilon(\mathbf u) - "
+                                                  "\frac 13 (\nabla \cdot \mathbf u) \mathbf I)$ * stress_strain_directors, and differs from the "
+                                                  "full stress by the absence of the pressure.  The second term in the "
+                                                  "difference is zero if the model is incompressible. "
+                                                  "If elasticity is used, the elastic contribution is being accounted for. ")
     }
   }
 }

cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc

@@ -787,42 +787,6 @@ namespace aspect
 
 
 
-
-
-
-    template <int dim>
-    void
-    CPO_AV_3D<dim>::parse_parameters (ParameterHandler &prm)
-    {
-      prm.enter_subsection("Material model");
-      {
-        prm.enter_subsection("CPO-induced Anisotropic Viscosity");
-        {
-
-          equation_of_state.parse_parameters (prm);
-          eta = prm.get_double("Reference viscosity");
-          min_strain_rate = prm.get_double("Minimum strain rate");
-          grain_size = prm.get_double("Grain size");
-          CnI_F = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for F")));
-          CnI_G = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for G")));
-          CnI_H = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for H")));
-          CnI_L = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for L")));
-          CnI_M = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for M")));
-          CnI_N = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for N")));
-          compositional_delta_rho    = prm.get_double ("Density differential for compositional field 1");
-        }
-        prm.leave_subsection();
-      }
-      prm.leave_subsection();
-
-
-
-      // Declare dependence
-      this->model_dependence.density = NonlinearDependence::compositional_fields;
-    }
-
-
-
     template <int dim>
     void
     CPO_AV_3D<dim>::declare_parameters (ParameterHandler &prm)
@@ -880,6 +844,38 @@ namespace aspect
 
 
 
+    template <int dim>
+    void
+    CPO_AV_3D<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Material model");
+      {
+        prm.enter_subsection("CPO-induced Anisotropic Viscosity");
+        {
+
+          equation_of_state.parse_parameters (prm);
+          eta = prm.get_double("Reference viscosity");
+          min_strain_rate = prm.get_double("Minimum strain rate");
+          grain_size = prm.get_double("Grain size");
+          CnI_F = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for F")));
+          CnI_G = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for G")));
+          CnI_H = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for H")));
+          CnI_L = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for L")));
+          CnI_M = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for M")));
+          CnI_N = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for N")));
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+
+
+
+      // Declare dependence
+      this->model_dependence.density = NonlinearDependence::compositional_fields;
+    }
+
+
+
     template <int dim>
     void
     CPO_AV_3D<dim>::create_additional_named_outputs(MaterialModel::MaterialModelOutputs<dim> &out) const

cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h

@@ -75,26 +75,56 @@ namespace aspect
     {
       public:
         void initialize() override;
+
         void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
                        MaterialModel::MaterialModelOutputs<dim> &out) const override;
+
+        bool is_compressible () const override;
+
         static void declare_parameters (ParameterHandler &prm);
+
         void parse_parameters (ParameterHandler &prm) override;
-        bool is_compressible () const override;
+        
         void create_additional_named_outputs(MaterialModel::MaterialModelOutputs<dim> &out) const override;
+
       private:
-        double eta; //reference viscosity
+
+        /**
+         * Reference viscosity.
+         */
+        double eta;
+
         /**
          * Defining a minimum strain rate stabilizes the viscosity calculation,
          * which involves a division by the strain rate. Units: 1/s.
          */
         double min_strain_rate;
+        
+        /**
+         * These are arrays that store eigenvalues of the olivine textures in a, b, and c axis, and 
+         * the olivine texture represented by Euler angles. For more details, please refer to the
+         * cpo_bingham_average particle property. To use the anisotropic viscosity plugin in this 
+         * cookbook, the CPO Bingham Average particle property must be included and Use rotation matrix
+         * must be set to false. The resulting arrays are:
+         * cpo_bingham_avg_a = [phi, eigenvalue 1 for a-axis, eigenvalue 2 for a-axis, eigenvalue 3 for a-axis]
+         * cpo_bingham_avg_b = [phi, eigenvalue 1 for b-axis, eigenvalue 2 for b-axis, eigenvalue 3 for b-axis]
+         * cpo_bingham_avg_c = [phi, eigenvalue 1 for c-axis, eigenvalue 2 for c-axis, eigenvalue 3 for c-axis]
+         * They are used in computing rotation matrix with regards to the CPO reference frame, and 
+         * the anisotropic Hill coefficients FGHLMN.
+         */
         std::vector<double> cpo_bingham_avg_a, cpo_bingham_avg_b, cpo_bingham_avg_c;
-        double grain_size;
+        
+        /**
+         * These are arrays that store coefficients used to compute the anisotropic Hill coefficients FGHLMN from
+         * a certain olivine texture represented with the eigenvalues of its a-, b-, and c-axis. Each array therefore contatains
+         * 9 coefficients and 1 constant. This empirical relationship is found using linear regression.
+         */
         std::vector<double> CnI_F, CnI_G, CnI_H, CnI_L, CnI_M, CnI_N;
-        double compositional_delta_rho;
 
-        unsigned int n_grains;
+        double grain_size;
+
         EquationOfState::LinearizedIncompressible<dim> equation_of_state;
+
         void set_assemblers(const SimulatorAccess<dim> &,
                             Assemblers::Manager<dim> &assemblers) const;
 

cookbooks/CPO_induced_anisotropic_viscosity/shearbox_cpo_av.prm

@@ -1,15 +1,18 @@
-#This cookbook demonstrate how CPO induced anisotropic viscosity evolves during simple shearing. The model consist of 1 cell and one particle in the middle. The particle includes 500 grains with initially random orientations (isotropic texture), which is then evolves into a CPO aligned with the shear direction
+# This cookbook demonstrates how CPO induced anisotropic viscosity evolves during 
+# simple shearing. The model consists of 1 cell and one particle in the middle. 
+# The particle includes 500 grains with initially random orientations (isotropic 
+# texture), which then evolves into a CPO aligned with the shear direction.
+
 set Additional shared libraries = ./plugin/libCPO_induced_anisotropic_viscosity.debug.so
-set CFL number = 0.1 ##Annotation CFL number
+set CFL number = 0.
 set Timing output frequency  = 10000
 set Dimension                = 3
 set Pressure normalization   = surface
 set Surface pressure         = 0
 set Nonlinear solver scheme  = single Advection, single Stokes
 set End time                 = 0.1
-# set Resume computation       = true
 set Use years instead of seconds = false
-set Output directory  = output_Shearbox_CPO_AV_1ts
+set Output directory  = output_Shearbox_CPO_AV
 
 subsection Initial composition model
   set List of model names = function
@@ -68,7 +71,10 @@ subsection Initial temperature model
   set Model name = function
 
   subsection Function
-  set Function expression = 1600 ## Annotation: Temperature function (note: anisotropic viscosity is not considered at T<1000K)
+
+  # Temperature function (anisotropic viscosity in this material is disabled at T<1000K)
+
+  set Function expression = 1600
   end
 end
 
@@ -87,10 +93,11 @@ subsection Boundary temperature model
 end
 
 subsection Boundary velocity model
-  set Prescribed velocity boundary indicators =  top:function, bottom:function, left:function, right:function ##, , front:function, back:function
+  set Prescribed velocity boundary indicators =  top:function, bottom:function, left:function, right:function
+  
   subsection Function
     set Variable names = x,y,z,t
-    set Function expression = (z-0.5);0;0 ## 13, top
+    set Function expression = (z-0.5);0;0
   end
 end
 
@@ -99,81 +106,90 @@ subsection Material model
   set Model name = CPO-induced Anisotropic Viscosity
 
   subsection CPO-induced Anisotropic Viscosity
-    set Reference viscosity = 1e9 # so that the initial stress is 0.3 MPa
+    # Chosen so that the initial stress is 0.3 MPa
+
+    set Reference viscosity = 1e9
   end
 end
 
 subsection Mesh refinement
   set Initial global refinement                = 0
   set Initial adaptive refinement              = 0
-  set Time steps between mesh refinement       = 10000
+  # For this single cell shear box test, there is no need for AMR
+
+  set Time steps between mesh refinement       = 0
 end
 
 subsection Postprocess
     set List of postprocessors = velocity statistics, composition statistics, memory statistics, visualization, particles, crystal preferred orientation
 
     subsection Visualization
         set Time between graphical output = 0.1
-        set List of output variables = material properties, strain rate,  named additional outputs,  shear stress,  stress
+        set List of output variables = material properties, strain rate, named additional outputs, shear stress, stress
 
       subsection Material properties
         set List of material properties = density, viscosity
       end
     end
+
     subsection Particles
         set Time between data output = 0.1
         set Data output format       = gnuplot, vtu
         set Exclude output properties = rotation_matrix, volume fraction
     end
+
     subsection Crystal Preferred Orientation
         set Time between data output = 0.1
         set Write in background thread = true
         set Compress cpo data files = false
-        set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles
-        set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles
+        set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles
+        set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles
     end
 end
 
 subsection Particles
-    set List of particle properties =  integrated strain, integrated strain invariant, velocity, pT path, strain rate, crystal preferred orientation, cpo bingham average, cpo elastic tensor #integrated strain, integrated strain invariant, velocity, pT path, strain rate, velocity gradient, stress, crystal preferred orientation, cpo bingham average, cpo bingham average euler angle, cpo elastic tensor
+    set List of particle properties =  integrated strain, integrated strain invariant, velocity, pT path, strain rate, crystal preferred orientation, cpo bingham average, cpo elastic tensor
     set Allow cells without particles = false
     set Interpolation scheme = nearest neighbor
     set Particle generator name = ascii file
+
     subsection Generator
         subsection Ascii file
             set Data directory = ./
             set Data file name = particle_one.dat
         end
     end
+
     subsection Crystal Preferred Orientation
+
         subsection Initial grains
             set Minerals = Olivine: A-fabric
             set Volume fractions minerals = 1.0
         end
-        set Number of grains per particle = 1000 ## Annotation grain count
-        set Property advection method = Backward Euler ## Annotation Property advection method
+
+        set Number of grains per particle = 1000
+        set Property advection method = Backward Euler
         set Property advection tolerance = 1e-15
         set CPO derivatives algorithm = D-Rex 2004
+
         subsection D-Rex 2004
-            set Mobility = 10 ## Annotation CPO Mobility
+            set Mobility = 10
             set Stress exponents = 3.5
             set Exponents p = 1.5
-            set Threshold GBS = 0.3 ## Annotation TGBS
+            set Threshold GBS = 0.3
         end
     end
+
     subsection CPO Bingham Average
         set Random number seed = 200
         set Use rotation matrix = false
     end
 end
 
 subsection Solver parameters
-  set Temperature solver tolerance = 1e-10
+  # set Temperature solver tolerance = 1e-10
+
   subsection Stokes solver parameters
     set Stokes solver type = block AMG
   end
 end
-
-# subsection Checkpointing
-#   set Steps between checkpoint = 50
-# end