add viscosity limit to rft

Author: danieldouglas92

doc/modules/changes/20260305_danieldouglas92

@@ -0,0 +1,4 @@
+Added: A parameter to the reactive fluid transport model that limits
+the amount of fluid weakening on the solid viscosity.
+<br>
+(Daniel Douglas, 2026/03/05)

include/aspect/material_model/reactive_fluid_transport.h

@@ -164,6 +164,7 @@ namespace aspect
         double shear_to_bulk_viscosity_ratio;
         double min_compaction_visc;
         double max_compaction_visc;
+        double min_weakened_viscosity;
         double reference_permeability;
         double alpha_phi;
         double reference_T;

source/material_model/reactive_fluid_transport.cc

@@ -191,7 +191,8 @@ namespace aspect
               for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
                 {
                   const double porosity = std::max(in.composition[q][porosity_idx],0.0);
-                  out.viscosities[q] *= (1.0 - porosity) * std::exp(- alpha_phi * porosity);
+                  const double weakened_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
+                  out.viscosities[q] = std::max(weakened_viscosity, min_weakened_viscosity);
                 }
             }
 
@@ -377,6 +378,9 @@ namespace aspect
                              boost::lexical_cast<std::string>(std::numeric_limits<double>::max()),
                              Patterns::Double (0),
                              "Upper cutoff for the compaction viscosity. Units: $\\text{Pa}\\text{s}$.");
+          prm.declare_entry ("Minimum weakened viscosity", "1e17",
+                             Patterns::Double (0),
+                             "Lower cutoff for the viscosity of the solid phase in the presence of fluids. Units: $\\text{Pa}\\text{s}$.");
           prm.declare_entry ("Reference fluid viscosity", "10",
                              Patterns::Double (0),
                              "The value of the constant melt/fluid viscosity $\\eta_f$. Units: $\\text{Pa}\\text{s}$.");
@@ -448,6 +452,7 @@ namespace aspect
           shear_to_bulk_viscosity_ratio     = prm.get_double ("Shear to bulk viscosity ratio");
           max_compaction_visc               = prm.get_double ("Maximum compaction viscosity");
           min_compaction_visc               = prm.get_double ("Minimum compaction viscosity");
+          min_weakened_viscosity            = prm.get_double ("Minimum weakened viscosity");
           eta_f                             = prm.get_double ("Reference fluid viscosity");
           reference_permeability            = prm.get_double ("Reference permeability");
           alpha_phi                         = prm.get_double ("Exponential fluid weakening factor");

tests/reactive_fluid_transport_min_viscosity.prm

@@ -0,0 +1,32 @@
+# This model tests the minimum viscosity within the reactive fluid transport model
+# correctly stops the viscosity from decreasing below a threshold in the presence
+# of fluids. The model prescribes a porosity of 10 % and an exponential fluid 
+# weakening factor of 30. Without a minimum value for the viscosity, this would 
+# result in the material viscosity of 4.481e17 Pa s, which is below the minimum
+# value that is specified in the Visco Plastic base material model. By setting
+# the parameter Minimum weakened viscosity to 1e10 Pa s, we ensure that the viscosity 
+# of the material does not drop below this value, even in the presence of fluids.
+include $ASPECT_SOURCE_DIR/tests/reactive_fluid_transport_zero_solubility.prm
+
+set Output directory = reactive_fluid_transport_min_viscosity
+set End time         = 0
+
+# 10 % porosity initially.
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.1; 0.0
+  end
+end
+
+subsection Material model
+  subsection Reactive Fluid Transport Model
+    set Exponential fluid weakening factor = 30
+    set Minimum weakened viscosity         = 1e19
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = material statistics
+end

tests/reactive_fluid_transport_min_viscosity/screen-output

@@ -0,0 +1,22 @@
+
+Number of active cells: 1 (on 1 levels)
+Number of degrees of freedom: 74 (18+7+18+4+9+9+9)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 0+0 iterations.
+   Solving fluid velocity system... 1 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.3813e-16, 4.89778e-17, 0, 1.55195e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 1.55195e-09
+
+
+   Postprocessing:
+     Average density / Average viscosity / Total mass:  3000 kg/m^3, 1e+19 Pa s, 1.2e+10 kg
+
+Termination requested by criterion: end time
+
+
+

tests/reactive_fluid_transport_min_viscosity/statistics

@@ -0,0 +1,18 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Iterations for temperature solver
+# 10: Iterations for composition solver 1
+# 11: Iterations for composition solver 2
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Average density (kg/m^3)
+# 16: Average viscosity (Pa s)
+# 17: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 1 22 9 18 1 0 0 0 0 0 0 3.00000000e+03 1.00000000e+19 1.20000000e+10