Add include

include/aspect/utilities.h

@@ -1060,6 +1060,13 @@ namespace aspect
      * entry in the list must match one of the allowed operations.
      */
     std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
+
+    /**
+     * Create matrix with unit at independent indices
+     */
+    template <int dim>
+    SymmetricTensor<2,dim> symmetric_independent_component_matrix (const unsigned int k);
+
   }
 }
 

source/material_model/visco_plastic.cc

@@ -22,6 +22,7 @@
 #include <aspect/utilities.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/base/signaling_nan.h>
+#include <aspect/newton.h>
 
 namespace aspect
 {
@@ -354,13 +355,24 @@ namespace aspect
     evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
              MaterialModel::MaterialModelOutputs<dim> &out) const
     {
+      // set up additional output for the derivatives
+      MaterialModel::MaterialModelDerivatives<dim> *derivatives;
+      derivatives = out.template get_additional_output<MaterialModel::MaterialModelDerivatives<dim> >();
+
+      // If we use the Newton solver, get the derivative scaling factor,
+      // otherwise set it to zero.
+      double derivative_scaling_factor = 0.0;
+      if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::Newton_Stokes)
+        derivative_scaling_factor = this->get_newton_handler().get_newton_derivative_scaling_factor();
+
       // Loop through points
       for (unsigned int i=0; i < in.temperature.size(); ++i)
         {
           const double temperature = in.temperature[i];
           const double pressure = in.pressure[i];
           const std::vector<double> &composition = in.composition[i];
           const std::vector<double> volume_fractions = compute_volume_fractions(composition);
+          const SymmetricTensor<2,dim> strain_rate = in.strain_rate[i];
 
           // Averaging composition-field dependent properties
 
@@ -397,14 +409,82 @@ namespace aspect
               // TODO: This is only consistent with viscosity averaging if the arithmetic averaging
               // scheme is chosen. It would be useful to have a function to calculate isostress viscosities.
               const std::vector<double> composition_viscosities =
-                calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, in.strain_rate[i],viscous_flow_law,yield_mechanism);
+                calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism);
+
+              std::vector<SymmetricTensor<2,dim> > composition_viscosities_derivatives(volume_fractions.size());
+              std::vector<double> composition_dviscosities_dpressure(volume_fractions.size());
 
               // The isostrain condition implies that the viscosity averaging should be arithmetic (see above).
               // We have given the user freedom to apply alternative bounds, because in diffusion-dominated
               // creep (where n_diff=1) viscosities are stress and strain-rate independent, so the calculation
               // of compositional field viscosities is consistent with any averaging scheme.
               out.viscosities[i] = average_value(composition, composition_viscosities, viscosity_averaging);
 
+              // compute derivatives if nessesary
+              if (derivative_scaling_factor != 0 && derivatives != NULL)
+                {
+                  const double finite_difference_accuracy = 1e-7;
+
+                  // For each independent component, compute the derivative.
+                  for (unsigned int independent_component_k = 0; independent_component_k < SymmetricTensor<2,dim>::n_independent_components; independent_component_k++)
+                    {
+                      const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (independent_component_k);
+                      SymmetricTensor<2,dim> strain_rate_component = strain_rate + std::max(std::fabs(strain_rate[indices_ij]), min_strain_rate) * finite_difference_accuracy
+                                                                     * Utilities::symmetric_independent_component_matrix<dim>(independent_component_k);
+                      std::vector<double> eta_component = calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate_component,viscous_flow_law,yield_mechanism);
+
+                      // For each composition of the independent component, compute the derivative.
+                      for (unsigned int composition_i = 0; composition_i < eta_component.size(); composition_i++)
+                        {
+                          // compute the difference between the viscosity with and without the strain-rate difference.
+                          double derivative_component = eta_component[composition_i] - composition_viscosities[composition_i];
+                          if (derivative_component != 0)
+                            {
+                              // when the difference is non-zero, divide by the difference.
+                              derivative_component /= std::max(std::fabs(strain_rate_component[indices_ij]), min_strain_rate) * finite_difference_accuracy;
+                            }
+                          composition_viscosities_derivatives[composition_i][indices_ij] = derivative_component;
+                        }
+                    }
+
+                  /**
+                   * Now compute the derivative of the viscosity to the pressure
+                   */
+                  double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);
+
+                  std::vector<double> pressure_difference_eta = calculate_isostrain_viscosities(volume_fractions, pressure_difference, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism);
+
+
+                  for (unsigned int composition_i = 0; composition_i < pressure_difference_eta.size(); composition_i++)
+                    {
+                      double deriv_pressure = pressure_difference_eta[composition_i] - composition_viscosities[composition_i];
+                      if (pressure_difference_eta[composition_i] != 0)
+                        {
+                          if (in.pressure[i] != 0)
+                            {
+                              deriv_pressure /= std::fabs(in.pressure[i]) * finite_difference_accuracy;
+                            }
+                          else
+                            {
+                              deriv_pressure = 0;
+                            }
+                        }
+                      composition_dviscosities_dpressure[composition_i] = deriv_pressure;
+                    }
+
+                  double viscosity_averaging_p = 0; // Geometric
+                  if (viscosity_averaging == harmonic)
+                    viscosity_averaging_p = -1;
+                  if (viscosity_averaging == arithmetic)
+                    viscosity_averaging_p = 1;
+                  if (viscosity_averaging == maximum_composition)
+                    viscosity_averaging_p = 1000;
+
+
+                  derivatives->viscosity_derivative_wrt_strain_rate[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i],volume_fractions, composition_viscosities, composition_viscosities_derivatives, viscosity_averaging_p);
+                  derivatives->viscosity_derivative_wrt_pressure[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i],volume_fractions, composition_viscosities, composition_dviscosities_dpressure, viscosity_averaging_p);
+
+                }
             }
 
           out.densities[i] = density;
@@ -432,7 +512,7 @@ namespace aspect
           double e_ii = 0.;
           if  (use_strain_weakening == true && use_finite_strain_tensor == false && this->get_timestep_number() > 0)
             {
-              edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate);
+              edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(strain_rate)))),min_strain_rate);
               e_ii = edot_ii*this->get_timestep();
               // Update reaction term
               out.reaction_terms[i][0] = e_ii;
@@ -739,7 +819,6 @@ namespace aspect
       {
         prm.enter_subsection ("Visco Plastic");
         {
-
           // Reference and minimum/maximum values
           reference_T = prm.get_double("Reference temperature");
           min_strain_rate = prm.get_double("Minimum strain rate");
@@ -770,6 +849,7 @@ namespace aspect
           if (use_strain_weakening)
             AssertThrow(this->n_compositional_fields() >= 1,
                         ExcMessage("There must be at least one compositional field. "));
+
           use_finite_strain_tensor  = prm.get_bool ("Use finite strain tensor");
           if (use_finite_strain_tensor)
             AssertThrow(this->n_compositional_fields() >= s,

source/utilities.cc

@@ -2536,6 +2536,21 @@ namespace aspect
     }
 
 
+    template <int dim>
+    SymmetricTensor<2,dim>
+    symmetric_independent_component_matrix (const unsigned int k)
+    {
+      // somehow I need this, because otherwise it complains that I am passing 3 arguments and it only takes 2...
+      const bool boolian = k < SymmetricTensor<2,dim>::n_independent_components;
+      Assert(boolian, ExcMessage("The component is larger then the amount of independent components in the matrix.") );
+
+      const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (k);
+
+      Tensor<2,dim> result;
+      result[indices_ij] = 1;
+
+      return symmetrize(result);
+    }
 
 
 // Explicit instantiations
@@ -2622,5 +2637,8 @@ namespace aspect
     template std_cxx11::array<double,2> convert_point_to_array<2>(const Point<2> &point);
     template std_cxx11::array<double,3> convert_point_to_array<3>(const Point<3> &point);
 
+    template SymmetricTensor<2,2> symmetric_independent_component_matrix (const unsigned int k);
+    template SymmetricTensor<2,3> symmetric_independent_component_matrix (const unsigned int k);
+
   }
 }