Merge pull request #6144 from danieldouglas92/fix_water_fugacity_with…

…_melt_transport

Use the adiabatic temperature when calculating water fugacity during initialization

source/material_model/rheology/compositional_viscosity_prefactors.cc

@@ -23,6 +23,8 @@
 #include <aspect/utilities.h>
 #include <aspect/global.h>
 #include <aspect/geometry_model/interface.h>
+#include <aspect/adiabatic_conditions/interface.h>
+
 
 #include <deal.II/base/signaling_nan.h>
 #include <deal.II/base/parameter_handler.h>
@@ -60,13 +62,28 @@ namespace aspect
               // We calculate the atomic H/Si ppm (C_OH) at each point to compute the water fugacity of
               // olivine assuming a composition of 90 mol% Forsterite and 10 mol% Fayalite from Hirth
               // and Kohlstaedt 2004 10.1029/138GM06.
+              const double temperature_for_fugacity = (this->simulator_is_past_initialization())
+                                                      ?
+                                                      in.temperature[q]
+                                                      :
+                                                      this->get_adiabatic_conditions().temperature(in.position[q]);
+              AssertThrow(temperature_for_fugacity != 0, ExcMessage(
+                            "The temperature used in the calculation for determining the water fugacity is zero. "
+                            "This is not allowed, because this value is used to divide through. It is probably "
+                            "being caused by the temperature being zero somewhere in the model. The relevant "
+                            "values for debugging are: temperature (" + Utilities::to_string(in.temperature[q]) +
+                            "), and adiabatic temperature ("
+                            + Utilities::to_string(this->get_adiabatic_conditions().temperature(in.position[q])) +
+                            "). If the adiabatic temperature is 0, double check that you are correctly defining an "
+                            " 'Adiabatic conditions model'."));
+
               const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
               const double mass_fraction_H2O = std::max(minimum_mass_fraction_water_for_dry_creep[composition_index], in.composition[q][bound_fluid_idx]); // mass fraction of bound water
               const double mass_fraction_olivine = 1 - mass_fraction_H2O; // mass fraction of olivine
               const double COH = (mass_fraction_H2O/molar_mass_H2O) / (mass_fraction_olivine/molar_mass_olivine) * 1e6; // COH in H / Si ppm
               const double point_water_fugacity = COH / A_H2O *
                                                   std::exp((activation_energy_H2O + in.pressure[q]*activation_volume_H2O)/
-                                                           (constants::gas_constant * in.temperature[q]));
+                                                           (constants::gas_constant * temperature_for_fugacity));
               const double r = modified_flow_laws == diffusion
                                ?
                                -diffusion_water_fugacity_exponents[composition_index]

tests/hk04_olivine_composite_hydration_prefactor.prm

@@ -102,24 +102,31 @@ subsection Material model
   set Model name = visco plastic
 
   subsection Visco Plastic
+    # Composite creep, and choose the Hirth & Kohlstaedt 2004 olivine hydration
+    # viscosity prefactor scheme.
+    set Viscous flow law                               = composite
+    set Viscosity prefactor scheme                     = HK04 olivine hydration
+    set Reference strain rate                          = 1e-17
+    set Minimum strain rate                            = 1e-17
+
     # Dislocation creep parameters
     # The water fugacity exponent is 1.2, but we need to divide by
     # the stress exponent, so we input r/n=0.34285714285714286.
-    set Reference strain rate                          = 1e-17
-    set Minimum strain rate                            = 1e-17
-    set Viscous flow law                               = composite
     set Prefactors for dislocation creep               = 1.0095317511683098e-25
     set Stress exponents for dislocation creep         = 3.5
     set Activation energies for dislocation creep      = 520e3
     set Activation volumes for dislocation creep       = 22e-6
     set Water fugacity exponents for dislocation creep = 0.34285714285714286
+
+    # Diffusion creep parameters
+    # The same approach is required as for the dislocation parameters,
+    # but the stress exponent for diffusion creep is 1.
     set Prefactors for diffusion creep                 = 1.5773933612004842e-21
     set Stress exponents for diffusion creep           = 1.0
     set Activation energies for diffusion creep        = 375e3
     set Activation volumes for diffusion creep         = 10e-6
     set Grain size                                     = 1e-3
     set Water fugacity exponents for diffusion creep   = 0.7
-    set Viscosity prefactor scheme                     = HK04 olivine hydration
   end
 end
 

tests/hk04_olivine_composite_hydration_prefactor_with_melt.prm

@@ -0,0 +1,70 @@
+# This test is almost identical to the test hk04_olivine_composite_hydration_prefactor,
+# the only difference being that this model includes melt transport. This test ensures
+# that the water fugacity calculation uses the temperature defined in the Adiabatic
+# conditions model during initialization.
+include $ASPECT_SOURCE_DIR/tests/hk04_olivine_composite_hydration_prefactor.prm
+
+# The reactive fluid transport model requires that there are compositional fields
+# named bound_fluid and porosity
+subsection Compositional fields
+  set Names of fields  = bound_fluid, porosity
+  set Number of fields = 2
+end
+
+# set a bound_fluid of 1% everywhere, and a porosity of 0% everywhere. The bound_fluid
+# composition is what determines the water fugacity.
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.01; 0
+  end
+end
+
+# Define an adiabatic conditions model, the temperature defined here is what will be used
+# during initialization. We set the adiabatic temperature to 1173 K, the adiabatic pressure
+# to 1 GPa, and the density to 3300 kg/m3.
+subsection Adiabatic conditions model
+  set Model name = function
+  subsection Function
+    set Function expression = 1173; 1e9; 3300
+  end
+end
+
+# Include melt transport.
+subsection Melt settings
+  set Include melt transport = true
+end
+
+subsection Material model
+  # Choose the reactive fluid transport model, which utilizes the melt framework.
+  set Model name =  reactive fluid transport
+
+  subsection Visco Plastic
+    # Composite creep, and choose the Hirth & Kohlstaedt 2004 olivine hydration
+    # viscosity prefactor scheme.
+    set Viscous flow law                               = composite
+    set Viscosity prefactor scheme                     = HK04 olivine hydration
+    set Reference strain rate                          = 1e-17
+    set Minimum strain rate                            = 1e-17
+
+    # Dislocation creep parameters
+    # The water fugacity exponent is 1.2, but we need to divide by
+    # the stress exponent, so we input r/n=0.34285714285714286.
+    set Prefactors for dislocation creep               = 1.0095317511683098e-25
+    set Stress exponents for dislocation creep         = 3.5
+    set Activation energies for dislocation creep      = 520e3
+    set Activation volumes for dislocation creep       = 22e-6
+    set Water fugacity exponents for dislocation creep = 0.34285714285714286
+
+    # Diffusion creep parameters
+    # The same approach is required as for the dislocation parameters,
+    # but the stress exponent for diffusion creep is 1.
+    set Prefactors for diffusion creep                 = 1.5773933612004842e-21
+    set Stress exponents for diffusion creep           = 1.0
+    set Activation energies for diffusion creep        = 375e3
+    set Activation volumes for diffusion creep         = 10e-6
+    set Grain size                                     = 1e-3
+    set Water fugacity exponents for diffusion creep   = 0.7
+  end
+end

tests/hk04_olivine_composite_hydration_prefactor_with_melt/screen-output

@@ -0,0 +1,27 @@
+
+Number of active cells: 100 (on 1 levels)
+Number of degrees of freedom: 3,629 (882+421+882+121+441+441+441)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving bound_fluid system ... 0 iterations.
+   Skipping porosity composition solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 5+0 iterations.
+   Solving fluid velocity system... 1 iterations.
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
+
+
+WARNING: The nonlinear solver in the current timestep failed to converge.
+Acting according to the parameter 'Nonlinear solver failure strategy'...
+Continuing to the next timestep even though solution is not fully converged.
+   Postprocessing:
+     Average density / Average viscosity / Total mass:  3198 kg/m^3, 2.628e+20 Pa s, 3.198e+13 kg
+     Writing graphical output:                          output-hk04_olivine_composite_hydration_prefactor_with_melt/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+
+
+WARNING: During this computation 1 nonlinear solver failures occurred!