Add compositional dependence on viscosity in the Steinberger material model

include/aspect/material_model/steinberger.h

@@ -216,7 +216,7 @@ namespace aspect
         double thermal_conductivity (const double temperature,
                                      const double pressure,
                                      const Point<dim> &position) const;
-
+        MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;
         /**
          * Whether the compositional fields representing mass fractions
          * should be normalized to one when computing their fractions
@@ -276,6 +276,14 @@ namespace aspect
         std::vector<double> saturation_scaling;
         double maximum_conductivity;
 
+
+        /**
+         * viscosity_prefactors are dimensionless quantities which
+         * are multiplied by the reference viscosity profile based on
+         * volume fraction
+         */
+        std::vector<double> viscosity_prefactors;
+
         /**
          * Information about lateral temperature averages.
          */

source/material_model/steinberger.cc

@@ -170,13 +170,15 @@ namespace aspect
     Steinberger<dim>::
     viscosity (const double temperature,
                const double /*pressure*/,
-               const std::vector<double> &,
+               const std::vector<double> &volume_fractions,
                const SymmetricTensor<2,dim> &,
                const Point<dim> &position) const
     {
       const double depth = this->get_geometry_model().depth(position);
       const double adiabatic_temperature = this->get_adiabatic_conditions().temperature(position);
 
+
+
       double delta_temperature;
       if (use_lateral_average_temperature)
         {
@@ -187,13 +189,24 @@ namespace aspect
         delta_temperature = temperature-adiabatic_temperature;
 
       // For an explanation on this formula see the Steinberger & Calderwood 2006 paper
-      const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
+      // We here compute the lateral variation of viscosity due to temperature (temperature_prefactor) as
+      // V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))] as in Eq. 6 of the paper.
+      // We get H/nR from the lateral_viscosity_lookup->lateral_viscosity function.
+      const double log_temperature_prefactor = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
+
       // Limit the lateral viscosity variation to a reasonable interval
-      const double vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),max_lateral_eta_variation),1/max_lateral_eta_variation);
+      const double temperature_prefactor = std::max(std::min(std::exp(log_temperature_prefactor),max_lateral_eta_variation),1/max_lateral_eta_variation);
+
+      //Visc_rT = exp[(H/nR)/T_adiabatic], Eq. 7 of the paper
+      const double eta_ref = radial_viscosity_lookup->radial_viscosity(depth);
 
-      const double vis_radial = radial_viscosity_lookup->radial_viscosity(depth);
+// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
+// std::cout<<"size vis_pre"<<viscosity_prefactors.size()<<std::endl;
+      const double compositional_prefactor = MaterialUtilities::average_value (volume_fractions, viscosity_prefactors, viscosity_averaging);
+      // std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
 
-      return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
+      // Radial viscosity profile is multiplied with lateral and compositional viscosity variation
+      return std::max(std::min(temperature_prefactor * eta_ref * compositional_prefactor,max_eta),min_eta);
     }
 
 
@@ -271,9 +284,6 @@ namespace aspect
 
       for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
         {
-          if (in.requests_property(MaterialProperties::viscosity))
-            out.viscosities[i] = viscosity(in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
-
           out.thermal_conductivities[i] = thermal_conductivity(in.temperature[i], in.pressure[i], in.position[i]);
           for (unsigned int c=0; c<in.composition[i].size(); ++c)
             out.reaction_terms[i][c] = 0;
@@ -295,7 +305,7 @@ namespace aspect
                   chemical_compositions.push_back(in.composition[i][c]);
 
               mass_fractions = MaterialUtilities::compute_composition_fractions(chemical_compositions, *composition_mask);
-
+// std::cout<<"mf"<<mass_fractions<<std::endl;
               // The function compute_volumes_from_masses expects as many mass_fractions as densities.
               // But the function compute_composition_fractions always adds another element at the start
               // of the vector that represents the background field. If there is no lookup table for
@@ -308,6 +318,11 @@ namespace aspect
                                                                                eos_outputs[i].densities,
                                                                                true);
 
+
+          if (in.requests_property(MaterialProperties::viscosity))
+            out.viscosities[i] = viscosity(in.temperature[i], in.pressure[i], volume_fractions[i], in.strain_rate[i], in.position[i]);
+
+
           MaterialUtilities::fill_averaged_equation_of_state_outputs(eos_outputs[i], mass_fractions, volume_fractions[i], i, out);
           fill_prescribed_outputs(i, volume_fractions[i], in, out);
         }
@@ -454,6 +469,17 @@ namespace aspect
                              Patterns::Double (0.),
                              "The maximum thermal conductivity that is allowed in the "
                              "model. Larger values will be cut off.");
+          prm.declare_entry ("Viscosity averaging scheme", "harmonic",
+                             Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"),
+                             "When more than one compositional field is present at a point "
+                             "with different viscosities, we need to come up with an average "
+                             "viscosity at that point. Select a weighted harmonic, arithmetic, "
+                             "geometric, or maximum composition.");
+          prm.declare_entry ("Viscosity prefactors", "1",
+                             Patterns::Anything(),
+                             "List of dimensionless quantities for background mantle and compositional fields,"
+                             "for a total of N+1 values, where N is the number of compositional fields."
+                             "If only one value is given, then all use the same value. Units: \\si{\\pascal\\second}.");
 
           // Table lookup parameters
           EquationOfState::ThermodynamicTableLookup<dim>::declare_parameters(prm);
@@ -483,6 +509,9 @@ namespace aspect
           max_eta              = prm.get_double ("Maximum viscosity");
           max_lateral_eta_variation    = prm.get_double ("Maximum lateral viscosity variation");
           thermal_conductivity_value = prm.get_double ("Thermal conductivity");
+          viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme",
+                                prm);
+
 
           // Rheological parameters
           if (prm.get ("Thermal conductivity formulation") == "constant")
@@ -551,7 +580,15 @@ namespace aspect
                                   + " fields of type chemical composition."));
 
           has_background_field = (equation_of_state.number_of_lookups() == n_chemical_fields + 1);
-
+          std::vector<std::string> list_of_composition_names = this->introspection().get_composition_names();
+          if (equation_of_state.number_of_lookups() == 1)
+            list_of_composition_names.resize(1, "background");
+          if ((equation_of_state.number_of_lookups() == 1) && (has_background_field))
+            list_of_composition_names.resize(0, "background");
+          viscosity_prefactors = Utilities::parse_map_to_double_array (prm.get("Viscosity prefactors"),
+                                                                       list_of_composition_names,
+                                                                       has_background_field,
+                                                                       "Viscosity prefactors");
           prm.leave_subsection();
         }
         prm.leave_subsection();

tests/multicomponent_steinberger.prm

@@ -0,0 +1,128 @@
+## Test the multicomponent steinberger material model. Three blobs with different viscosities move around. The reference viscosity of background is from the radial viscosity profile of Steinberger & Calderwood 2006.
+
+set Dimension = 2
+set CFL number                             = 1.0
+set End time                               = 60000
+set Resume computation                     = false
+set Start time                             = 0
+set Adiabatic surface temperature          = 1600
+set Surface pressure                       = 0
+set Pressure normalization = surface
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme                = single Advection, single Stokes
+
+
+subsection Boundary temperature model
+  set List of model names = initial temperature
+  set Fixed temperature boundary indicators = top, bottom, left, right
+end
+
+
+subsection Discretization
+  set Stokes velocity polynomial degree       = 2
+  set Temperature polynomial degree           = 2
+  set Use locally conservative discretization = false
+  subsection Stabilization parameters
+    set alpha = 2
+    set beta  = 0.078
+    set cR    = 0.5
+  end
+end
+
+
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X extent = 28.e5
+    set Y extent = 29.e5
+    set X repetitions = 4
+  end
+end
+
+
+subsection Gravity model
+  set Model name = vertical
+  subsection Vertical
+    set Magnitude = 10.0
+  end
+end
+
+
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function expression =   1600
+  end
+end
+
+
+subsection Material model
+set Model name = Steinberger
+subsection Steinberger model
+    set Use lateral average temperature for viscosity = false
+    set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
+    set Radial viscosity file name = radial-visc.txt
+    set Lateral viscosity file name = temp-viscosity-prefactor.txt
+    set Viscosity prefactors = 1, 1, 10, 0.1
+   set Material file names   =  pyr-ringwood88.txt, pyr-ringwood88.txt , pyr-ringwood88.txt, pyr-ringwood88.txt
+   set Minimum viscosity = 1e19
+   set Maximum viscosity = 1e24
+ end
+end
+
+
+subsection Mesh refinement
+  set Additional refinement times        =
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 4
+  set Refinement fraction                = 0.3
+  set Coarsening fraction                = 0.05
+  set Strategy                           = composition, viscosity
+  set Time steps between mesh refinement = 5
+end
+
+
+# The parameters below this comment were created by the update script
+# as replacement for the old 'Model settings' subsection. They can be
+# safely merged with any existing subsections with the same name.
+
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+end
+
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function expression = if( sqrt( (x-5.e5 )^2 + (y-10.e5)^2) < 5.e4 , 1,0); if( sqrt((x-12.e5)^2 + (y-12.e5)^2 ) < 5.e4, 1,0); if( sqrt((x-17.e5)^2 + (y-20.e5)^2) < 5.e4, 1, 0)
+  end
+end
+
+subsection Termination criteria
+  set Termination criteria = end step
+  set End step                  = 2
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables = viscosity,density
+    set Number of grouped files       = 1
+    set Output format                 = vtu
+    set Time between graphical output = 60000
+  end
+end
+
+subsection Solver parameters
+  subsection Stokes solver parameters
+    set Linear solver tolerance = 1.e-7
+    set Number of cheap Stokes solver steps = 30
+  end
+end

tests/multicomponent_steinberger/screen-output

@@ -0,0 +1,39 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+Number of active cells: 1,024 (on 5 levels)
+Number of degrees of freedom: 26,647 (8,514+1,105+4,257+4,257+4,257+4,257)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 0 iterations.
+   Solving C_2 system ... 0 iterations.
+   Solving C_3 system ... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 30+48 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-multicomponent_steinberger/solution/solution-00000
+
+*** Timestep 1:  t=60000 years, dt=60000 years
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 8 iterations.
+   Solving C_2 system ... 8 iterations.
+   Solving C_3 system ... 9 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 30+38 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-multicomponent_steinberger/solution/solution-00001
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
++----------------------------------+-----------+------------+------------+
++----------------------------------+-----------+------------+------------+
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------

tests/multicomponent_steinberger/statistics

@@ -0,0 +1,17 @@
+# 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: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for composition solver 2
+# 11: Iterations for composition solver 3
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Visualization file name
+0 0.000000000000e+00 0.000000000000e+00 1024 9619 4257 12771 0 0 0 0 78 1801 216 output-multicomponent_steinberger/solution/solution-00000 
+1 6.000000000000e+04 6.000000000000e+04 1024 9619 4257 12771 0 8 8 9 68 1448 209 output-multicomponent_steinberger/solution/solution-00001