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,22 @@ namespace aspect
         delta_temperature = temperature-adiabatic_temperature;
 
       // For an explanation on this formula see the Steinberger & Calderwood 2006 paper
+      // We here compute the lateral variation of viscosity due to temperature (vis_lateral) 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 vis_lateral_exp = -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);
 
+      //Visc_rT = exp[(H/nR)/T_adiabatic], Eq. 7 of the paper
       const double vis_radial = radial_viscosity_lookup->radial_viscosity(depth);
 
-      return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
+
+      const double vis_compositional = MaterialUtilities::average_value (volume_fractions, viscosity_prefactors, viscosity_averaging);
+
+      // Radial viscosity profile is multiplied with lateral and compositional viscosity variation
+      return std::max(std::min(vis_lateral * vis_radial * vis_compositional,max_eta),min_eta);
     }
 
 
@@ -271,9 +282,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;
@@ -308,6 +316,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 +467,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 +507,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 +578,12 @@ namespace aspect
                                   + " fields of type chemical composition."));
 
           has_background_field = (equation_of_state.number_of_lookups() == n_chemical_fields + 1);
+          const std::vector<std::string> list_of_composition_names = this->introspection().get_composition_names();
 
+          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 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        = 1
+  set Initial global refinement          = 5
+  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,82 @@
+-----------------------------------------------------------------------------
+-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
+--     . version 2.6.0-pre (compositional_viscosity_steinberger, 6f88ba556)
+--     . using deal.II 9.4.0
+--     .       with 32 bit indices and vectorization level 0 (64 bits)
+--     . using Trilinos 12.14.1
+--     . using p4est 2.2.0
+--     . running in OPTIMIZED mode
+--     . running with 1 MPI process
+-----------------------------------------------------------------------------
+
+-----------------------------------------------------------------------------
+-- For information on how to cite ASPECT, see:
+--   https://aspect.geodynamics.org/citing.html?ver=2.6.0-pre&sha=6f88ba556&src=code
+-----------------------------------------------------------------------------
+Number of active cells: 4,096 (on 6 levels)
+Number of degrees of freedom: 104,487 (33,410+4,257+16,705+16,705+16,705+16,705)
+
+*** 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+60 iterations.
+
+Number of active cells: 4,453 (on 7 levels)
+Number of degrees of freedom: 121,518 (38,850+4,968+19,425+19,425+19,425+19,425)
+
+*** 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+57 iterations.
+
+   Postprocessing:
+     Writing graphical output: /home/poulami/software/Aspect_globalscale/aspect/tests/multicomponent_steinberger/solution/solution-00000
+
+*** Timestep 1:  t=60000 years, dt=60000 years
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 9 iterations.
+   Solving C_2 system ... 8 iterations.
+   Solving C_3 system ... 9 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 30+39 iterations.
+
+   Postprocessing:
+     Writing graphical output: /home/poulami/software/Aspect_globalscale/aspect/tests/multicomponent_steinberger/solution/solution-00001
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start     |      39.3s |            |
+|                                              |            |            |
+| Section                          | no. calls |  wall time | % of total |
++----------------------------------+-----------+------------+------------+
+| Assemble Stokes system           |         3 |     0.815s |       2.1% |
+| Assemble composition system      |         9 |      2.65s |       6.8% |
+| Assemble temperature system      |         3 |     0.832s |       2.1% |
+| Build Stokes preconditioner      |         3 |     0.737s |       1.9% |
+| Build composition preconditioner |         9 |     0.203s |      0.52% |
+| Build temperature preconditioner |         3 |    0.0636s |      0.16% |
+| Initialization                   |         1 |     0.529s |       1.3% |
+| Postprocessing                   |         2 |     0.315s |       0.8% |
+| Refine mesh structure, part 1    |         1 |     0.526s |       1.3% |
+| Refine mesh structure, part 2    |         1 |    0.0393s |       0.1% |
+| Setup dof systems                |         2 |    0.0935s |      0.24% |
+| Setup initial conditions         |         2 |     0.483s |       1.2% |
+| Setup matrices                   |         2 |      0.51s |       1.3% |
+| Solve Stokes system              |         3 |      31.3s |        80% |
+| Solve composition system         |         9 |    0.0598s |      0.15% |
+| Solve temperature system         |         3 |    0.0111s |         0% |
++----------------------------------+-----------+------------+------------+
+
+-- Total wallclock time elapsed including restarts: 39s
+-----------------------------------------------------------------------------
+-- For information on how to cite ASPECT, see:
+--   https://aspect.geodynamics.org/citing.html?ver=2.6.0-pre&sha=6f88ba556&src=code
+-----------------------------------------------------------------------------

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 4453 43818 19425 58275 0 0 0 0 87 2333 421 /home/poulami/software/Aspect_globalscale/aspect/tests/multicomponent_steinberger/solution/solution-00000 
+1 6.000000000000e+04 6.000000000000e+04 4453 43818 19425 58275 0 9 8 9 69 1600 319 /home/poulami/software/Aspect_globalscale/aspect/tests/multicomponent_steinberger/solution/solution-00001