Add fluid pressure to yield strenght

include/aspect/material_model/rheology/drucker_prager.h

@@ -87,7 +87,8 @@ namespace aspect
           compute_yield_stress (const double cohesion,
                                 const double angle_internal_friction,
                                 const double pressure,
-                                const double max_yield_stress) const;
+                                const double max_yield_stress,
+                                const double fluid_pressure = 0.) const;
 
           /**
            * Compute the plastic viscosity with the yield stress and effective strain rate.
@@ -98,7 +99,8 @@ namespace aspect
                              const double pressure,
                              const double effective_strain_rate,
                              const double max_yield_stress,
-                             const double pre_yield_viscosity = std::numeric_limits<double>::infinity()) const;
+                             const double pre_yield_viscosity = std::numeric_limits<double>::infinity(),
+                             const double fluid_pressure = 0.) const;
 
           /**
            * Compute the strain rate and first stress derivative

include/aspect/material_model/rheology/strain_dependent.h

@@ -138,8 +138,7 @@ namespace aspect
                                       const int i,
                                       const double min_strain_rate,
                                       const bool plastic_yielding,
-                                      MaterialModel::MaterialModelOutputs<dim> &out,
-                                      const double old_strainrate = 0.) const;
+                                      MaterialModel::MaterialModelOutputs<dim> &out) const;
 
           /**
            * A function that returns a ComponentMask, which indicates that components

include/aspect/material_model/rheology/visco_plastic.h

@@ -114,7 +114,8 @@ namespace aspect
                                             const std::vector<double> &volume_fractions,
                                             const std::vector<double> &phase_function_values = std::vector<double>(),
                                             const std::vector<unsigned int> &n_phases_per_composition =
-                                              std::vector<unsigned int>()) const;
+                                              std::vector<unsigned int>(),
+                                            const double fluid_pressure = 0.) const;
 
           /**
            * A function that fills the viscosity derivatives in the

include/aspect/material_model/visco_plastic.h

@@ -23,7 +23,6 @@
 
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/interface.h>
-#include <aspect/material_model/melt_boukare.h>
 #include <aspect/material_model/melt_phipps_morgan.h>
 #include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 #include <aspect/material_model/rheology/visco_plastic.h>

source/material_model/rheology/drucker_prager.cc

@@ -34,18 +34,21 @@ namespace aspect
       DruckerPrager<dim>::compute_yield_stress (const double cohesion,
                                                 const double angle_internal_friction,
                                                 const double pressure,
-                                                const double max_yield_stress) const
+                                                const double max_yield_stress,
+                                                const double fluid_pressure) const
       {
         const double sin_phi = std::sin(angle_internal_friction);
         const double cos_phi = std::cos(angle_internal_friction);
         const double stress_inv_part = 1. / (std::sqrt(3.0) * (3.0 + sin_phi));
+        const double fluid_pressure_sin_phi = fluid_pressure > 0. ? (1. - fluid_pressure/pressure)*sin_phi : sin_phi;
 
         // Initial yield stress (no stabilization terms)
         const double yield_stress = ( (dim==3)
                                       ?
-                                      ( 6.0 * cohesion * cos_phi + 6.0 * pressure * sin_phi) * stress_inv_part
+                                      ( 6.0 * cohesion * cos_phi + 6.0 * pressure * fluid_pressure_sin_phi) * stress_inv_part
                                       :
-                                      cohesion * cos_phi + pressure * sin_phi);
+                                      cohesion * cos_phi + pressure * fluid_pressure_sin_phi);
+//        std::cout << "yield_stress " << yield_stress << " pressure " << pressure << " fluid_pressure " << fluid_pressure << std::endl;
 
         return std::min(yield_stress, max_yield_stress);
       }
@@ -59,9 +62,10 @@ namespace aspect
                                              const double pressure,
                                              const double effective_strain_rate,
                                              const double max_yield_stress,
-                                             const double pre_yield_viscosity) const
+                                             const double pre_yield_viscosity,
+                                             const double fluid_pressure) const
       {
-        const double yield_stress = compute_yield_stress(cohesion, angle_internal_friction, pressure, max_yield_stress);
+        const double yield_stress = compute_yield_stress(cohesion, angle_internal_friction, pressure, max_yield_stress, fluid_pressure);
 
         const double strain_rate_effective_inv = 1./(2.*effective_strain_rate);
 

source/material_model/rheology/strain_dependent.cc

@@ -485,9 +485,9 @@ namespace aspect
                              const int i,
                              const double min_strain_rate,
                              const bool plastic_yielding,
-                             MaterialModel::MaterialModelOutputs<dim> &out,
-                             const double old_strainrate) const
+                             MaterialModel::MaterialModelOutputs<dim> &out) const
       {
+
         // If strain weakening is used, overwrite the first reaction term,
         // which represents the second invariant of the (plastic) strain tensor.
         // If plastic strain is tracked (so not the total strain), only overwrite
@@ -496,18 +496,7 @@ namespace aspect
         // Calculate changes in strain and update the reaction terms
         if  (this->simulator_is_past_initialization() && this->get_timestep_number() > 0 && in.requests_property(MaterialProperties::reaction_terms))
           {
-            // If an iterated Advection scheme is used,
-            // use the strain rate of the last timestep to compute
-            // a constant reaction term within one timestep.
-            // TODO plastic_yielding can have changed per nonlinear iteration
-//           std::cout << "current and old strain rate " << std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate) << ", " << old_strainrate << std::endl;
-            double edot_ii = 0;
-            if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::iterated_Advection_and_Stokes
-                && this->get_nonlinear_iteration() > 0)
-              edot_ii = old_strainrate;
-            else
-              edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate);
-
+            const double edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate);
             double delta_e_ii = edot_ii*this->get_timestep();
 
             // Adjusting strain values to account for strain healing without exceeding an unreasonable range

source/material_model/rheology/visco_plastic.cc

@@ -98,7 +98,8 @@ namespace aspect
                                        const unsigned int i,
                                        const std::vector<double> &volume_fractions,
                                        const std::vector<double> &phase_function_values,
-                                       const std::vector<unsigned int> &n_phases_per_composition) const
+                                       const std::vector<unsigned int> &n_phases_per_composition,
+                                       const double fluid_pressure) const
       {
         IsostrainViscosities output_parameters;
 
@@ -254,7 +255,8 @@ namespace aspect
             const double yield_stress = drucker_prager_plasticity.compute_yield_stress(current_cohesion,
                                                                                        current_friction,
                                                                                        pressure_for_plasticity,
-                                                                                       drucker_prager_parameters.max_yield_stress);
+                                                                                       drucker_prager_parameters.max_yield_stress,
+                                                                                       fluid_pressure);
 
             // Step 4b: select if yield viscosity is based on Drucker Prager or stress limiter rheology
             double viscosity_yield = viscosity_pre_yield;
@@ -280,7 +282,8 @@ namespace aspect
                                                                                     pressure_for_plasticity,
                                                                                     current_edot_ii,
                                                                                     drucker_prager_parameters.max_yield_stress,
-                                                                                    viscosity_pre_yield);
+                                                                                    viscosity_pre_yield,
+                                                                                    fluid_pressure);
                       output_parameters.composition_yielding[j] = true;
                     }
                   break;

source/material_model/visco_plastic.cc

@@ -23,6 +23,8 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/base/signaling_nan.h>
 #include <aspect/newton.h>
+#include <aspect/melt.h>
+#include <aspect/simulator.h>
 #include <aspect/adiabatic_conditions/interface.h>
 #include <aspect/gravity_model/interface.h>
 #include <aspect/introspection.h>
@@ -181,29 +183,29 @@ namespace aspect
       // While the number of phases is fixed, the value of the phase function is updated for every point
       std::vector<double> phase_function_values(phase_function.n_phase_transitions(), 0.0);
 
-      std::vector<SymmetricTensor<2,dim> > old_strainrate_values;
-      std::vector<SymmetricTensor<2,dim> > old_old_strainrate_values;
-      if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::iterated_Advection_and_Stokes &&
-          this->get_timestep_number() > 0 &&
-          in.requests_property(MaterialProperties::reaction_terms) &&
+      // Get the fluid pressures to modify the yield strength
+      std::vector<double> fluid_pressures(in.position.size());
+      if (this->include_melt_transport() &&
+          in.requests_property(MaterialProperties::viscosity) &&
+          this->get_melt_handler().is_melt_cell(in.current_cell) &&
           in.current_cell.state() == IteratorState::valid)
         {
-          const QGauss<dim> quadrature_formula (this->get_fe()
-                                                .base_element(this->introspection().base_elements.velocities).degree+1);
-          const unsigned int n_q_points = quadrature_formula.size();
-          old_strainrate_values.resize(n_q_points);
-          old_old_strainrate_values.resize(n_q_points);
+          std::vector<Point<dim> > quadrature_positions(in.position.size());
+
+          for (unsigned int i=0; i < in.position.size(); ++i)
+            quadrature_positions[i] = this->get_mapping().transform_real_to_unit_cell(in.current_cell, in.position[i]);
+
+          // FEValues requires a quadrature and we provide the default quadrature
+          // as we only need to evaluate the solution and gradients.
           FEValues<dim> fe_values (this->get_mapping(),
                                    this->get_fe(),
-                                   quadrature_formula,
-                                   update_gradients   |
-                                   update_quadrature_points);
+                                   Quadrature<dim>(quadrature_positions),
+                                   update_values | update_gradients);
           fe_values.reinit (in.current_cell);
-          fe_values[this->introspection().extractors.velocities].get_function_symmetric_gradients (this->get_old_solution(),old_strainrate_values);
-          if (this->get_timestep_number() > 1)
-            {
-              fe_values[this->introspection().extractors.velocities].get_function_symmetric_gradients (this->get_old_old_solution(),old_old_strainrate_values);
-            }
+          // get fluid pressure from the current solution
+          const FEValuesExtractors::Scalar extractor_pressure = this->introspection().variable("fluid pressure").extractor_scalar();
+          fe_values[extractor_pressure].get_function_values (this->get_solution(),
+                                                             fluid_pressures);
         }
 
       // Loop through all requested points
@@ -285,7 +287,7 @@ 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 IsostrainViscosities isostrain_viscosities =
-                rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, phase_function.n_phase_transitions_for_each_composition());
+                rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, phase_function.n_phase_transitions_for_each_composition(),fluid_pressures[i]);
 
               // 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
@@ -309,22 +311,7 @@ namespace aspect
             }
 
           // Calculate changes in strain invariants and update the reaction terms
-          double old_strainrate  = 0.;
-          if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::iterated_Advection_and_Stokes &&
-              this->get_timestep_number() > 0 &&
-              in.requests_property(MaterialProperties::reaction_terms) &&
-              in.current_cell.state() == IteratorState::valid)
-            {
-              if (this->get_timestep_number() > 1)
-                {
-                  SymmetricTensor<2,dim> linearization_strainrate = (1 + this->get_timestep()/this->get_old_timestep()) * old_strainrate_values[i]
-                                                                    - this->get_timestep()/ this->get_old_timestep() * old_old_strainrate_values[i];
-                  old_strainrate  = std::max(sqrt(std::fabs(second_invariant(deviator(linearization_strainrate)))),rheology->min_strain_rate);
-                }
-              else
-                old_strainrate  = std::max(sqrt(std::fabs(second_invariant(deviator(old_strainrate_values[i])))),rheology->min_strain_rate);
-            }
-          rheology->strain_rheology.fill_reaction_outputs(in, i, rheology->min_strain_rate, plastic_yielding, out, old_strainrate);
+          rheology->strain_rheology.fill_reaction_outputs(in, i, rheology->min_strain_rate, plastic_yielding, out);
 
           // Fill plastic outputs if they exist.
           rheology->fill_plastic_outputs(i,volume_fractions,plastic_yielding,in,out);