Merge pull request #5970 from gassmoeller/adiabatic_pressure_plasticity

Add option to use adiabatic pressure to determine yield stress in visco plastic rheology

include/aspect/material_model/rheology/visco_plastic.h

@@ -282,6 +282,15 @@ namespace aspect
            */
           bool use_adiabatic_pressure_in_creep;
 
+          /**
+           * Whether to use the adiabatic pressure instead of the full pressure
+           * when calculating the plastic yield stress.
+           * This may be helpful in models where the full pressure has
+           * large variations resulting in solver convergence issues.
+           * Be aware that this setting will change the plastic shear band angle.
+           */
+          bool use_adiabatic_pressure_in_plasticity;
+
           /**
            * List of exponents controlling the behavior of the stress limiter
            * yielding mechanism.

source/material_model/rheology/visco_plastic.cc

@@ -316,8 +316,12 @@ namespace aspect
             // than the lithostatic pressure.
 
             double pressure_for_plasticity = in.pressure[i];
+
+            if (use_adiabatic_pressure_in_plasticity)
+              pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
+
             if (allow_negative_pressures_in_plasticity == false)
-              pressure_for_plasticity = std::max(in.pressure[i],0.0);
+              pressure_for_plasticity = std::max(pressure_for_plasticity,0.0);
 
             // Step 5a: calculate the Drucker-Prager yield stress
             const double yield_stress = drucker_prager_plasticity.compute_yield_stress(current_cohesion,
@@ -580,6 +584,14 @@ namespace aspect
                            "full pressure has an unusually large negative value arising from "
                            "large negative dynamic pressure, resulting in solver convergence "
                            "issue and in some cases a viscosity of zero.");
+        prm.declare_entry ("Use adiabatic pressure in plasticity", "false",
+                           Patterns::Bool (),
+                           "Whether to use the adiabatic pressure instead of the full "
+                           "pressure when calculating plastic yield stress. "
+                           "This may be helpful in models where the "
+                           "full pressure has unusually large variations, resulting "
+                           "in solver convergence issues. Be aware that this setting "
+                           "will change the plastic shear band angle.");
 
         // Diffusion creep parameters
         Rheology::DiffusionCreep<dim>::declare_parameters(prm);
@@ -714,6 +726,7 @@ namespace aspect
                                "'drucker prager' plasticity option."));
 
         allow_negative_pressures_in_plasticity = prm.get_bool ("Allow negative pressures in plasticity");
+        use_adiabatic_pressure_in_plasticity = prm.get_bool("Use adiabatic pressure in plasticity");
         use_adiabatic_pressure_in_creep = prm.get_bool("Use adiabatic pressure in creep viscosity");
 
         // Diffusion creep parameters
@@ -813,8 +826,12 @@ namespace aspect
             const double max_yield_stress = drucker_prager_plasticity.compute_drucker_prager_parameters(0).max_yield_stress;
 
             double pressure_for_plasticity = in.pressure[i];
+
+            if (use_adiabatic_pressure_in_plasticity)
+              pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
+
             if (allow_negative_pressures_in_plasticity == false)
-              pressure_for_plasticity = std::max(in.pressure[i], 0.0);
+              pressure_for_plasticity = std::max(pressure_for_plasticity, 0.0);
 
             // average over the volume volume fractions
             for (unsigned int j = 0; j < volume_fractions.size(); ++j)

tests/visco_plastic_adiabatic_pressure_in_plasticity.prm

@@ -0,0 +1,108 @@
+# A test for using the adiabatic pressure in plastic
+# viscosity calculation for the visco_plastic
+# material model. In this example, the adiabatic pressure is 3.3e9 Pa
+# (= rho * g * depth = 3300 * 10 * 100e3) at the bottom. So, the
+# predicted yield stress can be calculated by the
+# following equation:
+#
+#   eta_yield = cos(phi) * coh + sin(phi) * pressure
+#
+#  eta_yield         # yield stress (Pa)
+#  phi  = 10         # angle of internal friction (degrees)
+#  coh  = 1e6        # cohesion (Pa)
+#  p    = 3.3e9      # pressure (Pa)
+
+# Using these parameter we expect a yield stress of 5.7402e8 Pa at y=0 (P=3.3e9)
+# which is confirmed in the output. Note that the pressure in this model is
+# close to 0 (because there is no gravity), and we prescribe the adiabatic
+# pressure to the value above. So the fact that we match the correct yield
+# stress shows we are using the adiabatic pressure instead of actual pressure
+# to compute the yield stress.
+
+set Dimension                              = 2
+set End time                               = 0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme                = single Advection, single Stokes
+set Max nonlinear iterations               = 1
+
+subsection Adiabatic conditions model
+  set Model name =  function
+
+  subsection Function
+    set Variable names = x
+    set Function expression = 1; 3.3e4*x; 3300
+  end
+end
+
+# Model geometry (100x100 km, 10 km spacing)
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X repetitions = 10
+    set Y repetitions = 10
+    set X extent      = 100e3
+    set Y extent      = 100e3
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 0
+  set Time steps between mesh refinement = 0
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators   = bottom, top, left, right
+  set List of model names = box
+
+  subsection Box
+    set Bottom temperature = 1573
+    set Left temperature   = 1573
+    set Right temperature  = 1573
+    set Top temperature    = 1573
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = bottom, top, left, right
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 1573
+  end
+end
+
+subsection Material model
+  set Model name = visco plastic
+
+  subsection Visco Plastic
+    set Reference strain rate = 1.e-16
+    set Viscous flow law = composite
+
+    set Use adiabatic pressure in plasticity = true
+    set Angles of internal friction = 10
+    set Cohesions = 1e6
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0.0
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables = viscosity, named additional outputs, adiabat
+    set Output format            = gnuplot
+  end
+end

tests/visco_plastic_adiabatic_pressure_in_plasticity/screen-output

@@ -0,0 +1,15 @@
+
+Number of active cells: 100 (on 1 levels)
+Number of degrees of freedom: 1,444 (882+121+441)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving Stokes system... 25+0 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-visco_plastic_adiabatic_pressure_in_plasticity/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+