Update description and format

Author: Wang-yijun

cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc

@@ -120,7 +120,8 @@ namespace aspect
                                                   "\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. ")
+                                                  "This particle property plugin should only be activated when using the "
+                                                  "CPO induced anisotropic viscosity material model from the cookbook. ")
     }
   }
 }

cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.h

@@ -44,7 +44,6 @@ namespace aspect
        * \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.
        *
        * The member functions are all implementations of those declared in the
        * base class. See there for their meaning.

cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc

@@ -29,6 +29,7 @@
 #include <aspect/heating_model/interface.h>
 #include <aspect/gravity_model/interface.h>
 #include <aspect/simulator/assemblers/stokes.h>
+#include <aspect/simulator/assemblers/stokes_anisotropic_viscosity.h>
 #include <aspect/simulator_signals.h>
 #include <aspect/postprocess/particles.h>
 #include <aspect/introspection.h>
@@ -79,254 +80,6 @@ namespace aspect
   namespace MaterialModel
   {
     using namespace dealii;
-
-    /**
-     * Additional output fields for anisotropic viscosities to be added to
-     * the MaterialModel::MaterialModelOutputs structure and filled in the
-     * MaterialModel::Interface::evaluate() function.
-     */
-
-
-
-    namespace
-    {
-
-      template <int dim>
-      std::vector<std::string> make_AV_additional_outputs_names()
-      {
-        std::vector<std::string> names;
-
-        for (unsigned int i = 0; i < Tensor<4,dim>::n_independent_components ; ++i)
-          {
-            TableIndices<4> indices(Tensor<4,dim>::unrolled_to_component_indices(i));
-            names.push_back("anisotropic_viscosity"+std::to_string(indices[0]+1)+std::to_string(indices[1]+1)+std::to_string(indices[2]+1)+std::to_string(indices[3]+1));
-          }
-        return names;
-      }
-    }
-
-
-
-    template <int dim>
-    AnisotropicViscosity<dim>::AnisotropicViscosity (const unsigned int n_points)
-      :
-      NamedAdditionalMaterialOutputs<dim>(make_AV_additional_outputs_names<dim>()),
-      stress_strain_directors(n_points, dealii::identity_tensor<dim> ())
-    {}
-
-
-
-    template <int dim>
-    std::vector<double>
-    AnisotropicViscosity<dim>::get_nth_output(const unsigned int idx) const
-    {
-      std::vector<double> output(stress_strain_directors.size());
-      for (unsigned int i = 0; i < stress_strain_directors.size() ; ++i)
-        {
-          output[i]= stress_strain_directors[i][Tensor<4,dim>::unrolled_to_component_indices(idx)];
-        }
-      return output;
-    }
-  }
-}
-
-namespace aspect
-{
-  namespace Assemblers
-  {
-    template <int dim>
-    void
-    StokesPreconditionerAnisotropicViscosity<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesPreconditioner<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesPreconditioner<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesPreconditioner<dim>&> (data_base);
-
-      const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity
-        = scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points           = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-
-      // First loop over all dofs and find those that are in the Stokes system
-      // save the component (pressure and dim velocities) each belongs to.
-      for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
-        {
-          if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-            {
-              scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first;
-              ++i_stokes;
-            }
-          ++i;
-        }
-
-      // Loop over all quadrature points and assemble their contributions to
-      // the preconditioner matrix
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.grads_phi_u[i_stokes] =
-                    scratch.finite_element_values[introspection.extractors
-                                                  .velocities].symmetric_gradient(i, q);
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection
-                                                                          .extractors.pressure].value(i, q);
-                  ++i_stokes;
-                }
-              ++i;
-            }
-
-          const double eta = scratch.material_model_outputs.viscosities[q];
-          const double one_over_eta = 1. / eta;
-          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i)
-            for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j)
-              if (scratch.dof_component_indices[i] ==
-                  scratch.dof_component_indices[j])
-                data.local_matrix(i, j) += (2.0 * eta * (scratch.grads_phi_u[i]
-                                                         * stress_strain_director
-                                                         * scratch.grads_phi_u[j])
-                                            + one_over_eta * pressure_scaling
-                                            * pressure_scaling
-                                            * (scratch.phi_p[i]
-                                               * scratch.phi_p[j]))
-                                           * JxW;
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    StokesPreconditionerAnisotropicViscosity<dim>::
-    create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
-    {
-      const unsigned int n_points = outputs.viscosities.size();
-
-      if (outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
-        {
-          outputs.additional_outputs.push_back(
-            std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    StokesIncompressibleTermsAnisotropicViscosity<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
-
-      const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity
-        = scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-
-      const std::shared_ptr<const MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>> force
-        = scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>>();
-
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].value (i,q);
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
-                  if (scratch.rebuild_stokes_matrix)
-                    {
-                      scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q);
-                      scratch.div_phi_u[i_stokes]   = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q);
-                    }
-                  ++i_stokes;
-                }
-              ++i;
-            }
-          // Viscosity scalar
-          const double eta = (scratch.rebuild_stokes_matrix
-                              ?
-                              scratch.material_model_outputs.viscosities[q]
-                              :
-                              numbers::signaling_nan<double>());
-
-          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
-
-          const Tensor<1,dim>
-          gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
-
-          const double density = scratch.material_model_outputs.densities[q];
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
-            {
-              data.local_rhs(i) += (density * gravity * scratch.phi_u[i])
-                                   * JxW;
-
-              if (force != nullptr)
-                data.local_rhs(i) += (force->rhs_u[q] * scratch.phi_u[i]
-                                      + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i])
-                                     * JxW;
-
-              if (scratch.rebuild_stokes_matrix)
-                for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
-                  {
-                    data.local_matrix(i,j) += ( eta * 2.0 * (scratch.grads_phi_u[i] * stress_strain_director * scratch.grads_phi_u[j])
-                                                // assemble \nabla p as -(p, div v):
-                                                - (pressure_scaling *
-                                                   scratch.div_phi_u[i] * scratch.phi_p[j])
-                                                // assemble the term -div(u) as -(div u, q).
-                                                // Note the negative sign to make this
-                                                // operator adjoint to the grad p term:
-                                                - (pressure_scaling *
-                                                   scratch.phi_p[i] * scratch.div_phi_u[j]))
-                                              * JxW;
-                  }
-            }
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    StokesIncompressibleTermsAnisotropicViscosity<dim>::
-    create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
-    {
-      const unsigned int n_points = outputs.viscosities.size();
-
-      if (outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
-        {
-          outputs.additional_outputs.push_back(
-            std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
-        }
-
-      if (this->get_parameters().enable_additional_stokes_rhs
-          && outputs.template has_additional_output_object<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>>() == false)
-        {
-          outputs.additional_outputs.push_back(
-            std::make_unique<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>> (n_points));
-        }
-      Assert(!this->get_parameters().enable_additional_stokes_rhs
-             ||
-             outputs.template get_additional_output_object<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>>()->rhs_u.size()
-             == n_points, ExcInternalError());
-    }
   }
 
   namespace HeatingModel

cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h

@@ -116,8 +116,8 @@ namespace aspect
 
         /**
          * 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.
+         * a certain olivine texture represented with the eigenvalues of its a-, b-, and c-axis. Each array contains
+         * 9 coefficients and 1 constant.
          */
         std::vector<double> CnI_F, CnI_G, CnI_H, CnI_L, CnI_M, CnI_N;