From 6463b7fa4b440e6c971f0629ea30eef1b12f0eab Mon Sep 17 00:00:00 2001
From: Yimin Jin <jinym@cea-igp.ac.cn>
Date: Sun, 30 Jun 2024 22:37:16 +0800
Subject: [PATCH 1/4] Add dilation relevant terms to the Newton solver

---
 .../material_model/rheology/drucker_prager.h  |  17 +++
 .../material_model/rheology/visco_plastic.h   |   9 +-
 include/aspect/material_model/visco_plastic.h |   4 +-
 include/aspect/newton.h                       |  10 ++
 include/aspect/stokes_matrix_free.h           |  23 ++-
 .../material_model/rheology/drucker_prager.cc |  50 ++++++-
 .../material_model/rheology/visco_plastic.cc  | 136 ++++++++++++------
 source/material_model/visco_plastic.cc        |  15 +-
 source/simulator/assemblers/newton_stokes.cc  |  30 ++++
 source/simulator/core.cc                      |   4 +-
 source/simulator/newton.cc                    |   2 +
 source/simulator/stokes_matrix_free.cc        |  82 ++++++++---
 12 files changed, 306 insertions(+), 76 deletions(-)

diff --git a/include/aspect/material_model/rheology/drucker_prager.h b/include/aspect/material_model/rheology/drucker_prager.h
index d787f070357..904aa86e160 100644
--- a/include/aspect/material_model/rheology/drucker_prager.h
+++ b/include/aspect/material_model/rheology/drucker_prager.h
@@ -45,6 +45,11 @@ namespace aspect
          */
         double angle_internal_friction;
 
+        /**
+         * Dilation angle (psi) for the current composition and phase
+         */
+        double angle_dilation;
+
         /**
          * Cohesion for the current composition and phase
          */
@@ -138,9 +143,21 @@ namespace aspect
           compute_derivative (const double angle_internal_friction,
                               const double effective_strain_rate) const;
 
+          /**
+           * Compute the plastic dilation term, which is defined as the product of
+           * the plastic multiplier and the derivative of plastic potential with
+           * respect to minus pressure.
+           */
+          double
+          compute_plastic_dilation (const double angle_dilation,
+                                    const double non_yielding_viscosity,
+                                    const double effective_viscosity,
+                                    const double effective_strain_rate) const;
+
         private:
 
           std::vector<double> angles_internal_friction;
+          std::vector<double> angles_dilation;
           std::vector<double> cohesions;
           double max_yield_stress;
 
diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h
index e4f6b171102..9879c0e2939 100644
--- a/include/aspect/material_model/rheology/visco_plastic.h
+++ b/include/aspect/material_model/rheology/visco_plastic.h
@@ -108,6 +108,13 @@ namespace aspect
        * The current cohesion.
        */
       std::vector<double> current_cohesions;
+
+      /**
+       * The plastic dilation terms, defined as the product of the
+       * plastic multiplier and the derivative of plastic potential
+       * with respect to negative pressure.
+       */
+      std::vector<double> plastic_dilation;
     };
 
     namespace Rheology
@@ -149,7 +156,7 @@ namespace aspect
            */
           void compute_viscosity_derivatives(const unsigned int point_index,
                                              const std::vector<double> &volume_fractions,
-                                             const std::vector<double> &composition_viscosities,
+                                             const IsostrainViscosities &isostrain_values,
                                              const MaterialModel::MaterialModelInputs<dim> &in,
                                              MaterialModel::MaterialModelOutputs<dim> &out,
                                              const std::vector<double> &phase_function_values = std::vector<double>(),
diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h
index 9ba9bf9fa67..62637a82048 100644
--- a/include/aspect/material_model/visco_plastic.h
+++ b/include/aspect/material_model/visco_plastic.h
@@ -23,7 +23,7 @@
 
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/interface.h>
-#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
+#include <aspect/material_model/equation_of_state/multicomponent_compressible.h>
 #include <aspect/material_model/rheology/visco_plastic.h>
 
 #include<deal.II/fe/component_mask.h>
@@ -268,7 +268,7 @@ namespace aspect
         /**
          * Object for computing the equation of state.
          */
-        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
+        EquationOfState::MulticomponentCompressible<dim> equation_of_state;
 
         /**
          * Object that handles phase transitions.
diff --git a/include/aspect/newton.h b/include/aspect/newton.h
index f7e4fc70c46..591c24a73f8 100644
--- a/include/aspect/newton.h
+++ b/include/aspect/newton.h
@@ -63,6 +63,16 @@ namespace aspect
          * derivatives when material averaging is applied.
          */
         std::vector<double> viscosity_derivative_averaging_weights;
+
+        /**
+         * The derivatives of the prescribed dilation with respect to strain rate.
+         */
+        std::vector<SymmetricTensor<2,dim>> dilation_derivative_wrt_strain_rate;
+
+        /**
+         * The derivatives of the prescribed dilation term with respect to pressure.
+         */
+        std::vector<double> dilation_derivative_wrt_pressure;
     };
   }
 
diff --git a/include/aspect/stokes_matrix_free.h b/include/aspect/stokes_matrix_free.h
index e90f857cc4a..8e2be8c0380 100644
--- a/include/aspect/stokes_matrix_free.h
+++ b/include/aspect/stokes_matrix_free.h
@@ -125,6 +125,12 @@ namespace aspect
        */
       bool apply_stabilization_free_surface_faces;
 
+      /**
+       * If true, derivatives of the prescribed dilation are
+       * part of the operator.
+       */
+      bool enable_prescribed_dilation;
+
       /**
        * Table which stores viscosity values for each cell.
        *
@@ -144,7 +150,7 @@ namespace aspect
        * variables: viscosity derivative with respect to pressure,
        * the Newton derivative scaling factor, and the averaging weight.
        */
-      Table<2, VectorizedArray<number>> newton_factor_wrt_pressure_table;
+      Table<2, VectorizedArray<number>> viscosity_derivative_wrt_pressure_table;
 
       /**
        * Table which stores the product of the following four
@@ -154,7 +160,20 @@ namespace aspect
        * is PD or SPD, otherwise, it is 1.
        */
       Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
-      newton_factor_wrt_strain_rate_table;
+      viscosity_derivative_wrt_strain_rate_table;
+
+      /**
+       * Table which stores the product of the dilation derivative
+       * with respect to pressure and the Newton derivative scaling factor.
+       */
+      Table<2, VectorizedArray<number>> dilation_derivative_wrt_pressure_table;
+
+      /**
+       * Table which stores the product of the dilation derivative
+       * with respect to strain rate and the Newton derivative scaling factor.
+       */
+      Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
+      dilation_derivative_wrt_strain_rate_table;
 
       /**
        * Table which stores the product of the pressure perturbation
diff --git a/source/material_model/rheology/drucker_prager.cc b/source/material_model/rheology/drucker_prager.cc
index aa5856a7a5e..67436cc4a8e 100644
--- a/source/material_model/rheology/drucker_prager.cc
+++ b/source/material_model/rheology/drucker_prager.cc
@@ -34,6 +34,7 @@ namespace aspect
     {
       DruckerPragerParameters::DruckerPragerParameters()
         : angle_internal_friction (numbers::signaling_nan<double>()),
+          angle_dilation (numbers::signaling_nan<double>()),
           cohesion  (numbers::signaling_nan<double>()),
           max_yield_stress (numbers::signaling_nan<double>())
       {}
@@ -60,6 +61,7 @@ namespace aspect
           {
             // no phases
             drucker_prager_parameters.angle_internal_friction = angles_internal_friction[composition];
+            drucker_prager_parameters.angle_dilation = angles_dilation[composition];
             drucker_prager_parameters.cohesion = cohesions[composition];
           }
         else
@@ -67,6 +69,8 @@ namespace aspect
             // Average among phases
             drucker_prager_parameters.angle_internal_friction = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
                                                                 angles_internal_friction, composition);
+            drucker_prager_parameters.angle_dilation = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
+                                                       angles_dilation, composition);
             drucker_prager_parameters.cohesion = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
                                                  cohesions, composition);
           }
@@ -181,6 +185,26 @@ namespace aspect
 
 
 
+      template <int dim>
+      double
+      DruckerPrager<dim>::compute_plastic_dilation(const double angle_dilation,
+                                                   const double non_yielding_viscosity,
+                                                   const double effective_viscosity,
+                                                   const double effective_strain_rate) const
+      {
+        const double sin_psi = std::sin(angle_dilation);
+        const double minus_dQ_dp = (dim == 3 ?
+                                    6.0 * sin_psi / (std::sqrt(3.0) * (3.0 + sin_psi)) :
+                                    sin_psi);
+
+        const double plastic_multiplier =
+          std::max(0.0, (1.0 - effective_viscosity / non_yielding_viscosity) * 2.0 * effective_strain_rate);
+
+        return minus_dQ_dp * plastic_multiplier;
+      }
+
+
+
       template <int dim>
       void
       DruckerPrager<dim>::declare_parameters (ParameterHandler &prm)
@@ -192,6 +216,13 @@ namespace aspect
                            "those corresponding to chemical compositions. "
                            "For a value of zero, in 2d the von Mises criterion is retrieved. "
                            "Angles higher than 30 degrees are harder to solve numerically. Units: degrees.");
+        prm.declare_entry ("Angles of dilation", "0.",
+                           Patterns::Anything(),
+                           "List of angles of plastic dilation, $\\psi$, for background material and compositional fields, "
+                           "for a total of N+1 values, where N is the number of all compositional fields or only "
+                           "those corresponding to chemical compositions. "
+                           "For a value of zero, the von Mises flow rule is retrieved. "
+                           "The dilation angle should never exceed the internal friction angle.");
         prm.declare_entry ("Cohesions", "1e20",
                            Patterns::Anything(),
                            "List of cohesions, $C$, for background material and compositional fields, "
@@ -248,10 +279,25 @@ namespace aspect
 
         angles_internal_friction = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of internal friction"),
                                                                                     options);
+        angles_dilation = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of dilation"),
+                                                                           options);
 
         // Convert angles from degrees to radians
-        for (double &angle : angles_internal_friction)
-          angle *= constants::degree_to_radians;
+        for (unsigned int i = 0; i < angles_internal_friction.size(); ++i)
+          {
+            angles_internal_friction[i] *= constants::degree_to_radians;
+            angles_dilation[i]          *= constants::degree_to_radians;
+
+            AssertThrow(angles_dilation[i] <= angles_internal_friction[i],
+                        ExcMessage("The dilation angle is greater than the internal friction angle for "
+                                   "composition " + Utilities::to_string(i) + "."));
+
+            if (angles_dilation[i] > 0.0)
+              AssertThrow(this->get_parameters().enable_prescribed_dilation == true,
+                          ExcMessage("ASPECT detected a nonzero dilation angle, but dilation is "
+                                     "not enabled. Please set parameter entry 'Enable prescribed "
+                                     "dilation' to 'true'."));
+          }
 
         options.property_name = "Cohesions";
         cohesions = Utilities::MapParsing::parse_map_to_double_array(prm.get("Cohesions"),
diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc
index 2fe6259e7d2..55fe2cab300 100644
--- a/source/material_model/rheology/visco_plastic.cc
+++ b/source/material_model/rheology/visco_plastic.cc
@@ -109,6 +109,7 @@ namespace aspect
         output_parameters.composition_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());
         output_parameters.current_friction_angles.resize(volume_fractions.size(), numbers::signaling_nan<double>());
         output_parameters.current_cohesions.resize(volume_fractions.size(), numbers::signaling_nan<double>());
+        output_parameters.plastic_dilation.resize(volume_fractions.size(), numbers::signaling_nan<double>());
 
         // Assemble stress tensor if elastic behavior is enabled
         SymmetricTensor<2,dim> stress_old = numbers::signaling_nan<SymmetricTensor<2,dim>>();
@@ -383,7 +384,23 @@ namespace aspect
                                                                MaterialModel::MaterialUtilities::PhaseUtilities::logarithmic
                                                              );
             output_parameters.composition_viscosities[j] = std::min(std::max(effective_viscosity, minimum_viscosity_for_composition), maximum_viscosity_for_composition);
+
+            // Compute the plastic dilation if necessary.
+            if (this->get_parameters().enable_prescribed_dilation == true)
+              {
+                if (output_parameters.composition_yielding[j] == true)
+                  {
+                    const double current_dilation = drucker_prager_parameters.angle_dilation * weakening_factors[1];
+                    output_parameters.plastic_dilation[j] = drucker_prager_plasticity.compute_plastic_dilation (current_dilation,
+                                                            non_yielding_viscosity,
+                                                            output_parameters.composition_viscosities[j],
+                                                            effective_edot_ii);
+                  }
+                else
+                  output_parameters.plastic_dilation[j] = 0;
+              }
           }
+
         return output_parameters;
       }
 
@@ -394,7 +411,7 @@ namespace aspect
       ViscoPlastic<dim>::
       compute_viscosity_derivatives(const unsigned int i,
                                     const std::vector<double> &volume_fractions,
-                                    const std::vector<double> &composition_viscosities,
+                                    const IsostrainViscosities &current_isostrain_values,
                                     const MaterialModel::MaterialModelInputs<dim> &in,
                                     MaterialModel::MaterialModelOutputs<dim> &out,
                                     const std::vector<double> &phase_function_values,
@@ -403,16 +420,22 @@ namespace aspect
         MaterialModel::MaterialModelDerivatives<dim> *derivatives =
           out.template get_additional_output<MaterialModel::MaterialModelDerivatives<dim>>();
 
+        const bool enable_dilation = this->get_parameters().enable_prescribed_dilation;
+
         if (derivatives != nullptr)
           {
             // compute derivatives if necessary
-            std::vector<SymmetricTensor<2,dim>> composition_viscosities_derivatives(volume_fractions.size());
-            std::vector<double> composition_dviscosities_dpressure(volume_fractions.size());
+            const unsigned int n_compositions = volume_fractions.size();
+            std::vector<SymmetricTensor<2,dim>> composition_viscosity_derivatives_wrt_strain_rate(n_compositions);
+            std::vector<double> composition_viscosity_derivatives_wrt_pressure(n_compositions);
+
+            std::vector<SymmetricTensor<2,dim>> composition_dilation_derivatives_wrt_strain_rate(enable_dilation ? n_compositions : 0);
+            std::vector<double> composition_dilation_derivatives_wrt_pressure(enable_dilation ? n_compositions : 0);
 
             const double finite_difference_accuracy = 1e-7;
 
             // A new material model inputs variable that uses the strain rate and pressure difference.
-            MaterialModel::MaterialModelInputs<dim> in_derivatives = in;
+            MaterialModel::MaterialModelInputs<dim> in_forward = in;
 
             Assert(std::isfinite(in.strain_rate[i].norm()),
                    ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
@@ -425,59 +448,62 @@ namespace aspect
               {
                 const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component);
 
-                const SymmetricTensor<2,dim> strain_rate_difference = deviatoric_strain_rate
-                                                                      + std::max(std::fabs(deviatoric_strain_rate[strain_rate_indices]), min_strain_rate)
-                                                                      * finite_difference_accuracy
-                                                                      * Utilities::nth_basis_for_symmetric_tensors<dim>(component);
+                const double strain_rate_difference = std::max(std::fabs(deviatoric_strain_rate[strain_rate_indices]), min_strain_rate)
+                                                      * finite_difference_accuracy;
+                const SymmetricTensor<2,dim> forward_strain_rate = deviatoric_strain_rate
+                                                                   + strain_rate_difference * Utilities::nth_basis_for_symmetric_tensors<dim>(component);
 
-                in_derivatives.strain_rate[i] = strain_rate_difference;
+                in_forward.strain_rate[i] = forward_strain_rate;
 
-                std::vector<double> eta_component =
-                  calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
-                                                  phase_function_values, n_phase_transitions_per_composition).composition_viscosities;
+                const IsostrainViscosities forward_isostrain_values =
+                  calculate_isostrain_viscosities(in_forward, i, volume_fractions,
+                                                  phase_function_values, n_phase_transitions_per_composition);
 
                 // For each composition of the independent component, compute the derivative.
-                for (unsigned int composition_index = 0; composition_index < eta_component.size(); ++composition_index)
+                for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index)
                   {
-                    // compute the difference between the viscosity with and without the strain-rate difference.
-                    double viscosity_derivative = eta_component[composition_index] - composition_viscosities[composition_index];
-                    if (viscosity_derivative != 0)
-                      {
-                        // when the difference is non-zero, divide by the difference.
-                        viscosity_derivative /= std::max(std::fabs(strain_rate_difference[strain_rate_indices]), min_strain_rate)
-                                                * finite_difference_accuracy;
-                      }
-                    composition_viscosities_derivatives[composition_index][strain_rate_indices] = viscosity_derivative;
+                    composition_viscosity_derivatives_wrt_strain_rate[composition_index][strain_rate_indices]
+                      = ( forward_isostrain_values.composition_viscosities[composition_index] -
+                          current_isostrain_values.composition_viscosities[composition_index]
+                        ) / strain_rate_difference;
+
+                    if (enable_dilation)
+                      composition_dilation_derivatives_wrt_strain_rate[composition_index][strain_rate_indices]
+                        = ( forward_isostrain_values.plastic_dilation[composition_index] -
+                            current_isostrain_values.plastic_dilation[composition_index]
+                          ) / strain_rate_difference;
                   }
               }
 
             // Now compute the derivative of the viscosity to the pressure
-            const double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);
+            const double pressure_difference = std::fabs(in.pressure[i]) * finite_difference_accuracy;
+            const double forward_pressure = in.pressure[i] + pressure_difference;
 
-            in_derivatives.pressure[i] = pressure_difference;
+            in_forward.pressure[i] = forward_pressure;
 
             // Modify the in_derivatives object again to take the original strain rate.
-            in_derivatives.strain_rate[i] = in.strain_rate[i];
+            in_forward.strain_rate[i] = in.strain_rate[i];
 
-            const std::vector<double> viscosity_difference =
-              calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
-                                              phase_function_values, n_phase_transitions_per_composition).composition_viscosities;
+            const IsostrainViscosities forward_isostrain_values =
+              calculate_isostrain_viscosities(in_forward, i, volume_fractions,
+                                              phase_function_values, n_phase_transitions_per_composition);
 
-            for (unsigned int composition_index = 0; composition_index < viscosity_difference.size(); ++composition_index)
+            for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index)
               {
-                double viscosity_derivative = viscosity_difference[composition_index] - composition_viscosities[composition_index];
-                if (viscosity_difference[composition_index] != 0)
-                  {
-                    if (in.pressure[i] != 0)
-                      {
-                        viscosity_derivative /= std::fabs(in.pressure[i]) * finite_difference_accuracy;
-                      }
-                    else
-                      {
-                        viscosity_derivative = 0;
-                      }
-                  }
-                composition_dviscosities_dpressure[composition_index] = viscosity_derivative;
+                composition_viscosity_derivatives_wrt_pressure[composition_index]
+                  = pressure_difference > 0.0 ?
+                    ( forward_isostrain_values.composition_viscosities[composition_index] -
+                      current_isostrain_values.composition_viscosities[composition_index]
+                    ) / pressure_difference :
+                    0.0;
+
+                if (enable_dilation)
+                  composition_dilation_derivatives_wrt_pressure[composition_index]
+                    = pressure_difference > 0.0 ?
+                      ( forward_isostrain_values.plastic_dilation[composition_index] -
+                        current_isostrain_values.plastic_dilation[composition_index]
+                      ) / pressure_difference :
+                      0.0;
               }
 
             double viscosity_averaging_p = 0; // Geometric
@@ -491,15 +517,33 @@ namespace aspect
             derivatives->viscosity_derivative_wrt_strain_rate[i] =
               Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i],
                                                                volume_fractions,
-                                                               composition_viscosities,
-                                                               composition_viscosities_derivatives,
+                                                               current_isostrain_values.composition_viscosities,
+                                                               composition_viscosity_derivatives_wrt_strain_rate,
                                                                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,
+                                                               current_isostrain_values.composition_viscosities,
+                                                               composition_viscosity_derivatives_wrt_pressure,
                                                                viscosity_averaging_p);
+
+            if (enable_dilation)
+              {
+                double dummy_plastic_dilation = numbers::signaling_nan<double>();
+                // always use arithmetic average for plastic dilation terms
+                derivatives->dilation_derivative_wrt_strain_rate[i] =
+                  Utilities::derivative_of_weighted_p_norm_average(dummy_plastic_dilation,
+                                                                   volume_fractions,
+                                                                   current_isostrain_values.plastic_dilation,
+                                                                   composition_dilation_derivatives_wrt_strain_rate,
+                                                                   1);
+                derivatives->dilation_derivative_wrt_pressure[i] =
+                  Utilities::derivative_of_weighted_p_norm_average(dummy_plastic_dilation,
+                                                                   volume_fractions,
+                                                                   current_isostrain_values.plastic_dilation,
+                                                                   composition_dilation_derivatives_wrt_pressure,
+                                                                   1);
+              }
           }
       }
 
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index cef81e360e1..c0a2b33638d 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -247,7 +247,7 @@ namespace aspect
                     out.template get_additional_output<MaterialModel::MaterialModelDerivatives<dim>>())
 
                 rheology->compute_viscosity_derivatives(i, volume_fractions,
-                                                        isostrain_viscosities.composition_viscosities,
+                                                        isostrain_viscosities,
                                                         in, out, phase_function_values,
                                                         n_phase_transitions_for_each_chemical_composition);
             }
@@ -297,6 +297,14 @@ namespace aspect
                   elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i];
                 }
             }
+
+          // Compute the arithmetic average of plastic dilation term if necessary
+          if (PrescribedPlasticDilation<dim> *plastic_dilation =
+                out.template get_additional_output<PrescribedPlasticDilation<dim>>())
+            plastic_dilation->dilation[i]
+              = MaterialUtilities::average_value (volume_fractions,
+                                                  isostrain_viscosities.plastic_dilation,
+                                                  MaterialUtilities::arithmetic);
         }
 
       // If we use the full strain tensor, compute the change in the individual tensor components.
@@ -340,7 +348,7 @@ namespace aspect
         {
           MaterialUtilities::PhaseFunction<dim>::declare_parameters(prm);
 
-          EquationOfState::MulticomponentIncompressible<dim>::declare_parameters (prm);
+          EquationOfState::MulticomponentCompressible<dim>::declare_parameters (prm);
 
           Rheology::ViscoPlastic<dim>::declare_parameters(prm);
 
@@ -402,8 +410,7 @@ namespace aspect
 
           // Equation of state parameters
           equation_of_state.initialize_simulator (this->get_simulator());
-          equation_of_state.parse_parameters (prm,
-                                              std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+          equation_of_state.parse_parameters(prm);
 
           // Make options file for parsing maps to double arrays
           std::vector<std::string> chemical_field_names = this->introspection().chemical_composition_field_names();
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
index 174bab9cbc2..6a9c5cbab65 100644
--- a/source/simulator/assemblers/newton_stokes.cc
+++ b/source/simulator/assemblers/newton_stokes.cc
@@ -67,6 +67,8 @@ namespace aspect
         AssertThrow(elastic_out != nullptr,
                     ExcMessage("Error: The Newton method requires ElasticOutputs when elasticity is enabled."));
 
+      const bool enable_prescribed_dilation = this->get_parameters().enable_prescribed_dilation;
+
       // 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*/)
@@ -229,6 +231,17 @@ namespace aspect
                          )
                          * JxW;
                     }
+
+              if (enable_prescribed_dilation)
+                {
+                  for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+                    for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
+                      data.local_matrix(i,j) -= pressure_scaling * pressure_scaling
+                                                * scratch.phi_p[i] * scratch.phi_p[j]
+                                                * derivatives->dilation_derivative_wrt_pressure[q]
+                                                * JxW;
+
+                }
             }
         }
 #if DEBUG
@@ -524,6 +537,23 @@ namespace aspect
                                              " = " + Utilities::to_string(eta)));
                         }
                     }
+
+                  if (enable_prescribed_dilation)
+                    {
+                      std::vector<double> dilation_differentiations(stokes_dofs_per_cell);
+                      for (unsigned int k=0; k<stokes_dofs_per_cell; ++k)
+                        dilation_differentiations[k] =
+                          ( derivatives->dilation_derivative_wrt_strain_rate[q] * scratch.grads_phi_u[k]
+                            + derivatives->dilation_derivative_wrt_pressure[q] * pressure_scaling * scratch.phi_p[k]
+                          ) * derivative_scaling_factor;
+
+                      for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+                        for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
+                          data.local_matrix(i,j) += pressure_scaling
+                                                    * scratch.phi_p[i]
+                                                    * dilation_differentiations[j]
+                                                    * JxW;
+                    }
                 }
             }
         }
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index b1a35c03fa0..b9e5b480ac9 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -1009,7 +1009,9 @@ namespace aspect
             // For equal-order interpolation, we need a stabilization term
             // in the bottom right of Stokes matrix. Make sure we have the
             // necessary entries.
-            if (parameters.use_equal_order_interpolation_for_stokes == true)
+            if (parameters.use_equal_order_interpolation_for_stokes == true ||
+                (assemble_newton_stokes_system == true &&
+                 parameters.enable_prescribed_dilation == true))
               coupling[x.pressure][x.pressure] = DoFTools::always;
           }
       }
diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc
index 9297ebd3c15..d0fea7e73e8 100644
--- a/source/simulator/newton.cc
+++ b/source/simulator/newton.cc
@@ -37,6 +37,8 @@ namespace aspect
       : viscosity_derivative_wrt_pressure(n_points, numbers::signaling_nan<double>())
       , viscosity_derivative_wrt_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
       , viscosity_derivative_averaging_weights(n_points, numbers::signaling_nan<double>())
+      , dilation_derivative_wrt_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
+      , dilation_derivative_wrt_pressure(n_points, numbers::signaling_nan<double>())
     {}
   }
 
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 646e3d8437b..6858d601f8d 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -342,9 +342,11 @@ namespace aspect
     OperatorCellData<dim,number>::memory_consumption() const
     {
       return viscosity.memory_consumption()
-             + newton_factor_wrt_pressure_table.memory_consumption()
+             + viscosity_derivative_wrt_pressure_table.memory_consumption()
              + strain_rate_table.memory_consumption()
-             + newton_factor_wrt_strain_rate_table.memory_consumption();
+             + viscosity_derivative_wrt_strain_rate_table.memory_consumption()
+             + dilation_derivative_wrt_pressure_table.memory_consumption()
+             + dilation_derivative_wrt_strain_rate_table.memory_consumption();
     }
 
 
@@ -355,9 +357,11 @@ namespace aspect
     {
       enable_newton_derivatives = false;
       viscosity.clear();
-      newton_factor_wrt_pressure_table.clear();
+      viscosity_derivative_wrt_pressure_table.clear();
       strain_rate_table.clear();
-      newton_factor_wrt_strain_rate_table.clear();
+      viscosity_derivative_wrt_strain_rate_table.clear();
+      dilation_derivative_wrt_pressure_table.clear();
+      dilation_derivative_wrt_strain_rate_table.clear();
     }
   }
 
@@ -455,9 +459,19 @@ namespace aspect
             const VectorizedArray<number> val_p_q = p_eval.get_value(q);
 
             // Terms to be tested by phi_p:
-            const VectorizedArray<number> pressure_terms =
+            VectorizedArray<number> pressure_terms =
               -cell_data->pressure_scaling * div_u_q;
 
+            if (cell_data->enable_newton_derivatives &&
+                cell_data->enable_prescribed_dilation)
+              {
+                pressure_terms += ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q)
+                                    * sym_grad_u_q )
+                                  +
+                                  ( cell_data->dilation_derivative_wrt_pressure_table(cell,q)
+                                    * cell_data->pressure_scaling * val_p_q );
+              }
+
             // Terms to be tested by the symmetric gradients of phi_u:
             SymmetricTensor<2,dim,VectorizedArray<number>>
             velocity_terms = viscosity_x_2 * sym_grad_u_q;
@@ -477,9 +491,9 @@ namespace aspect
                 VectorizedArray<number> deta_dp_times_p(0.);
                 for (const unsigned int r : u_eval.quadrature_point_indices())
                   {
-                    deta_deps_times_sym_grad_u += cell_data->newton_factor_wrt_strain_rate_table(cell,r)
+                    deta_deps_times_sym_grad_u += cell_data->viscosity_derivative_wrt_strain_rate_table(cell,r)
                                                   * sym_grad_u[r];
-                    deta_dp_times_p += cell_data->newton_factor_wrt_pressure_table(cell,r) * val_p[r];
+                    deta_dp_times_p += cell_data->viscosity_derivative_wrt_pressure_table(cell,r) * val_p[r];
                     if (cell_data->symmetrize_newton_system)
                       eps_times_sym_grad_u += cell_data->strain_rate_table(cell,r) * sym_grad_u[r];
                   }
@@ -487,7 +501,7 @@ namespace aspect
                 velocity_terms +=
                   ( cell_data->symmetrize_newton_system ?
                     ( cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u +
-                      cell_data->newton_factor_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) :
+                      cell_data->viscosity_derivative_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) :
                     2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u )
                   +
                   2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p;
@@ -1375,8 +1389,15 @@ namespace aspect
           const unsigned int n_q_points = quadrature_formula.size();
 
           active_cell_data.strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points));
-          active_cell_data.newton_factor_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points));
-          active_cell_data.newton_factor_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points));
+          active_cell_data.viscosity_derivative_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points));
+          active_cell_data.viscosity_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points));
+
+          if (sim.parameters.enable_prescribed_dilation)
+            {
+              active_cell_data.enable_prescribed_dilation = true;
+              active_cell_data.dilation_derivative_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points));
+              active_cell_data.dilation_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points));
+            }
 
           for (unsigned int cell=0; cell<n_cells; ++cell)
             {
@@ -1435,12 +1456,12 @@ namespace aspect
                                             :
                                             1.0;
 
-                      active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]
+                      active_cell_data.viscosity_derivative_wrt_pressure_table(cell,q)[i]
                         = derivatives->viscosity_derivative_wrt_pressure[q] *
                           derivatives->viscosity_derivative_averaging_weights[q] *
                           newton_derivative_scaling_factor;
-                      Assert(std::isfinite(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]),
-                             ExcMessage("active_cell_data.newton_factor_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]) +
+                      Assert(std::isfinite(active_cell_data.viscosity_derivative_wrt_pressure_table(cell,q)[i]),
+                             ExcMessage("active_cell_data.viscosity_derivative_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.viscosity_derivative_wrt_pressure_table(cell,q)[i]) +
                                         ". Relevant variables are derivatives->viscosity_derivative_wrt_pressure[q] = " + std::to_string(derivatives->viscosity_derivative_wrt_pressure[q]) +
                                         ", derivatives->viscosity_derivative_averaging_weights[q] = " + std::to_string(derivatives->viscosity_derivative_averaging_weights[q]) +
                                         ", and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor)));
@@ -1451,16 +1472,39 @@ namespace aspect
                             active_cell_data.strain_rate_table(cell, q)[m][n][i]
                               = effective_strain_rate[m][n];
 
-                            active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]
+                            active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i]
                               = derivatives->viscosity_derivative_wrt_strain_rate[q][m][n] *
                                 derivatives->viscosity_derivative_averaging_weights[q] *
                                 newton_derivative_scaling_factor * alpha;
 
                             Assert(std::isfinite(active_cell_data.strain_rate_table(cell, q)[m][n][i]),
                                    ExcMessage("active_cell_data.strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.strain_rate_table(cell, q)[m][n][i])));
-                            Assert(std::isfinite(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]),
-                                   ExcMessage("active_cell_data.newton_factor_wrt_strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i])));
+                            Assert(std::isfinite(active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i]),
+                                   ExcMessage("active_cell_data.viscosity_derivative_wrt_strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i])));
                           }
+
+                      if (active_cell_data.enable_prescribed_dilation == true)
+                        {
+                          active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]
+                            = derivatives->dilation_derivative_wrt_pressure[q] *
+                              newton_derivative_scaling_factor;
+                          Assert(std::isfinite(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]),
+                                 ExcMessage("active_cell_data.dilation_derivative_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]) +
+                                            ". Relevant variables are derivatives->dilation_derivative_wrt_pressure[q] = " + std::to_string(derivatives->dilation_derivative_wrt_pressure[q]) +
+                                            " and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor)));
+
+                          for (unsigned int m=0; m<dim; ++m)
+                            for (unsigned int n=0; n<dim; ++n)
+                              {
+                                active_cell_data.dilation_derivative_wrt_strain_rate_table(cell,q)[m][n][i]
+                                  = derivatives->dilation_derivative_wrt_strain_rate[q][m][n] *
+                                    newton_derivative_scaling_factor;
+                                Assert(std::isfinite(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]),
+                                       ExcMessage("active_cell_data.dilation_derivative_wrt_strain_rate_table is not finite: " + std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]) +
+                                                  ". Relevant variables are derivatives->dilation_derivative_wrt_strain_rate[q] = " + std::to_string(derivatives->dilation_derivative_wrt_pressure[q]) +
+                                                  " and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor)));
+                              }
+                        }
                     }
                 }
             }
@@ -1476,9 +1520,11 @@ namespace aspect
           // delete data used for Newton derivatives if necessary
           // TODO: use Table::clear() once implemented in 10.0.pre
           active_cell_data.enable_newton_derivatives = false;
-          active_cell_data.newton_factor_wrt_pressure_table.reinit(TableIndices<2>(0,0));
+          active_cell_data.viscosity_derivative_wrt_pressure_table.reinit(TableIndices<2>(0,0));
           active_cell_data.strain_rate_table.reinit(TableIndices<2>(0,0));
-          active_cell_data.newton_factor_wrt_strain_rate_table.reinit(TableIndices<2>(0,0));
+          active_cell_data.viscosity_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(0,0));
+          active_cell_data.dilation_derivative_wrt_pressure_table.reinit(TableIndices<2>(0,0));
+          active_cell_data.dilation_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(0,0));
 
           for (unsigned int level=0; level<n_levels; ++level)
             level_cell_data[level].enable_newton_derivatives = false;

From b0e8fa41efb85b2e7f8cb50c346873948c9b601d Mon Sep 17 00:00:00 2001
From: Yimin Jin <jinym@cea-igp.ac.cn>
Date: Thu, 25 Jul 2024 16:54:36 +0800
Subject: [PATCH 2/4] fix bugs: missing pressure_scaling in local_apply()

---
 source/simulator/stokes_matrix_free.cc | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 6858d601f8d..0d6da231d95 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -465,11 +465,14 @@ namespace aspect
             if (cell_data->enable_newton_derivatives &&
                 cell_data->enable_prescribed_dilation)
               {
-                pressure_terms += ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q)
-                                    * sym_grad_u_q )
-                                  +
-                                  ( cell_data->dilation_derivative_wrt_pressure_table(cell,q)
-                                    * cell_data->pressure_scaling * val_p_q );
+                pressure_terms += (
+                                    ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q)
+                                      * sym_grad_u_q )
+                                    +
+                                    ( cell_data->dilation_derivative_wrt_pressure_table(cell,q)
+                                      * cell_data->pressure_scaling * val_p_q )
+                                  ) 
+                                  * cell_data->pressure_scaling;
               }
 
             // Terms to be tested by the symmetric gradients of phi_u:
@@ -504,7 +507,7 @@ namespace aspect
                       cell_data->viscosity_derivative_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) :
                     2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u )
                   +
-                  2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p;
+                  2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p * cell_data->pressure_scaling;
               }
 
             u_eval.submit_symmetric_gradient(velocity_terms, q);

From 251f1441a079a558e0037f4fd6e87fe9e99c15eb Mon Sep 17 00:00:00 2001
From: Yimin Jin <jinym@cea-igp.ac.cn>
Date: Fri, 26 Jul 2024 00:48:29 +0800
Subject: [PATCH 3/4] change the eos back to incompressible model

---
 include/aspect/material_model/visco_plastic.h | 4 ++--
 source/material_model/visco_plastic.cc        | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h
index 62637a82048..9ba9bf9fa67 100644
--- a/include/aspect/material_model/visco_plastic.h
+++ b/include/aspect/material_model/visco_plastic.h
@@ -23,7 +23,7 @@
 
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/interface.h>
-#include <aspect/material_model/equation_of_state/multicomponent_compressible.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 #include <aspect/material_model/rheology/visco_plastic.h>
 
 #include<deal.II/fe/component_mask.h>
@@ -268,7 +268,7 @@ namespace aspect
         /**
          * Object for computing the equation of state.
          */
-        EquationOfState::MulticomponentCompressible<dim> equation_of_state;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
 
         /**
          * Object that handles phase transitions.
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index c0a2b33638d..bb23225e649 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -348,7 +348,7 @@ namespace aspect
         {
           MaterialUtilities::PhaseFunction<dim>::declare_parameters(prm);
 
-          EquationOfState::MulticomponentCompressible<dim>::declare_parameters (prm);
+          EquationOfState::MulticomponentIncompressible<dim>::declare_parameters (prm);
 
           Rheology::ViscoPlastic<dim>::declare_parameters(prm);
 

From b1c748eb309a5f553833267f80c891e9921fbc4e Mon Sep 17 00:00:00 2001
From: Yimin Jin <jinym@cea-igp.ac.cn>
Date: Fri, 26 Jul 2024 00:49:29 +0800
Subject: [PATCH 4/4] treat plastic dilation in the same way as compressible
 model

---
 source/simulator/assemblers/newton_stokes.cc | 10 ----------
 source/simulator/assemblers/stokes.cc        | 14 --------------
 source/simulator/assembly.cc                 |  2 +-
 source/simulator/newton.cc                   |  3 ++-
 source/simulator/stokes_matrix_free.cc       |  9 ++++++---
 5 files changed, 9 insertions(+), 29 deletions(-)

diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
index 6a9c5cbab65..2f23efb9f86 100644
--- a/source/simulator/assemblers/newton_stokes.cc
+++ b/source/simulator/assemblers/newton_stokes.cc
@@ -451,16 +451,6 @@ namespace aspect
                                        * prescribed_dilation->dilation[q]
                                        * scratch.phi_p[i]
                                      ) * JxW;
-
-              // Only assemble this term if we are running incompressible, otherwise this term
-              // is already included on the LHS of the equation.
-              if (enable_prescribed_dilation && !material_model_is_compressible)
-                data.local_rhs(i) += (
-                                       // RHS of momentum eqn: - \int 2/3 eta R, div v
-                                       - 2.0 / 3.0 * eta
-                                       * prescribed_dilation->dilation[q]
-                                       * scratch.div_phi_u[i]
-                                     ) * JxW;
             }
 
           // and then the matrix, if necessary
diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc
index 164a7a7913e..935889fca30 100644
--- a/source/simulator/assemblers/stokes.cc
+++ b/source/simulator/assemblers/stokes.cc
@@ -374,10 +374,6 @@ namespace aspect
                     {
                       scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q);
                     }
-                  else if (prescribed_dilation && !material_model_is_compressible)
-                    {
-                      scratch.div_phi_u[i_stokes]   = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q);
-                    }
                   ++i_stokes;
                 }
               ++i;
@@ -426,16 +422,6 @@ namespace aspect
                                        * scratch.phi_p[i]
                                      ) * JxW;
 
-              // Only assemble this term if we are running incompressible, otherwise this term
-              // is already included on the LHS of the equation.
-              if (prescribed_dilation != nullptr && !material_model_is_compressible)
-                data.local_rhs(i) += (
-                                       // RHS of momentum eqn: - \int 2/3 eta R, div v
-                                       - 2.0 / 3.0 * eta
-                                       * prescribed_dilation->dilation[q]
-                                       * scratch.div_phi_u[i]
-                                     ) * JxW;
-
               if (scratch.rebuild_stokes_matrix)
                 for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
                   {
diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
index 5142237c083..3764efeddbd 100644
--- a/source/simulator/assembly.cc
+++ b/source/simulator/assembly.cc
@@ -74,7 +74,7 @@ namespace aspect
     assemblers->stokes_preconditioner.push_back(std::make_unique<aspect::Assemblers::StokesPreconditioner<dim>>());
     assemblers->stokes_system.push_back(std::make_unique<aspect::Assemblers::StokesIncompressibleTerms<dim>>());
 
-    if (material_model->is_compressible())
+    if (material_model->is_compressible() || parameters.enable_prescribed_dilation)
       {
         // The compressible part of the preconditioner is only necessary if we use the simplified A block
         if (parameters.use_full_A_block_preconditioner == false)
diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc
index d0fea7e73e8..49d2ea4d76d 100644
--- a/source/simulator/newton.cc
+++ b/source/simulator/newton.cc
@@ -52,7 +52,8 @@ namespace aspect
     assemblers.stokes_preconditioner.push_back(std::make_unique<aspect::Assemblers::NewtonStokesPreconditioner<dim>>());
     assemblers.stokes_system.push_back(std::make_unique<aspect::Assemblers::NewtonStokesIncompressibleTerms<dim>>());
 
-    if (this->get_material_model().is_compressible())
+    if (this->get_material_model().is_compressible() ||
+        this->get_parameters().enable_prescribed_dilation)
       {
         // The compressible part of the preconditioner is only necessary if we use the simplified A block
         if (this->get_parameters().use_full_A_block_preconditioner == false)
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 0d6da231d95..ac693e83920 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -482,7 +482,8 @@ namespace aspect
             for (unsigned int d=0; d<dim; ++d)
               velocity_terms[d][d] -= cell_data->pressure_scaling * val_p_q;
 
-            if (cell_data->is_compressible)
+            if (cell_data->is_compressible || 
+                cell_data->enable_prescribed_dilation)
               for (unsigned int d=0; d<dim; ++d)
                 velocity_terms[d][d] -= viscosity_x_2 / 3. * div_u_q;
 
@@ -851,7 +852,8 @@ namespace aspect
           velocity.get_symmetric_gradient (q);
         sym_grad_u *= viscosity_x_2;
 
-        if (cell_data->is_compressible)
+        if (cell_data->is_compressible ||
+            cell_data->enable_prescribed_dilation)
           {
             const VectorizedArray<number> div = trace(sym_grad_u);
             for (unsigned int d=0; d<dim; ++d)
@@ -1636,7 +1638,8 @@ namespace aspect
     // We never include Newton terms in step 0 and after that we solve with zero boundary conditions.
     // Therefore, we don't need to include Newton terms here.
 
-    const bool is_compressible = sim.material_model->is_compressible();
+    const bool is_compressible = sim.material_model->is_compressible() ||
+                                 sim.parameters.enable_prescribed_dilation;
 
     dealii::LinearAlgebra::distributed::BlockVector<double> rhs_correction(2);
     dealii::LinearAlgebra::distributed::BlockVector<double> u0(2);