Merge pull request geodynamics#4470 from jdannberg/geometric_average_…

…only_viscosity Add geometric average only viscosity operation
bangerth · Jan 29, 2022 · 012e56f · 012e56f
2 parents 27e31dc + fc97c2d
commit 012e56f
Show file tree

Hide file tree

Showing 9 changed files with 179 additions and 9 deletions.
diff --git a/doc/modules/changes/20220128b_jdannberg b/doc/modules/changes/20220128b_jdannberg
@@ -0,0 +1,5 @@
+New: There is now an averaging operation "geometric average only viscosity"
+that can be used to average the viscosity values of the different
+points in a cell. This can be used together with the GMG solver.
+<br>
+(Juliane Dannberg, 2022/01/28)
diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h
@@ -689,9 +689,9 @@ namespace aspect
        * quadrature point to \f[ \bar x = {10}^{\frac 1Q \sum_{q=1}^Q \log_{10} x_q} \f]
        * where $x_q$ are the values at the $Q$ quadrature points.
        *
-       * - Harmonic average only viscosity and project to Q1 only viscosity: Like
-       * harmonic averaging and project to Q1, but only
-       * applied to the viscosity.
+       * - Harmonic average only viscosity, Geometric  average only viscosity
+       * and project to Q1 only viscosity: Like harmonic averaging, geometric
+       * averaging and project to Q1, but only applied to the viscosity.
        */
       enum AveragingOperation
       {
@@ -703,6 +703,7 @@ namespace aspect
         project_to_Q1,
         log_average,
         harmonic_average_only_viscosity,
+        geometric_average_only_viscosity,
         project_to_Q1_only_viscosity
       };
 

diff --git a/source/material_model/interface.cc b/source/material_model/interface.cc
@@ -466,7 +466,7 @@ namespace aspect
     {
       std::string get_averaging_operation_names ()
       {
-        return "none|arithmetic average|harmonic average|geometric average|pick largest|project to Q1|log average|harmonic average only viscosity|project to Q1 only viscosity";
+        return "none|arithmetic average|harmonic average|geometric average|pick largest|project to Q1|log average|harmonic average only viscosity|geometric average only viscosity|project to Q1 only viscosity";
       }
 
 
@@ -488,6 +488,8 @@ namespace aspect
           return log_average;
         else if (s == "harmonic average only viscosity")
           return harmonic_average_only_viscosity;
+        else if (s == "geometric average only viscosity")
+          return geometric_average_only_viscosity;
         else if (s == "project to Q1 only viscosity")
           return project_to_Q1_only_viscosity;
         else
@@ -822,6 +824,13 @@ namespace aspect
             return;
           }
 
+        if (operation == geometric_average_only_viscosity)
+          {
+            average_property (geometric_average, projection_matrix, expansion_matrix,
+                              values_out.viscosities);
+            return;
+          }
+
         if (operation == project_to_Q1_only_viscosity)
           {
             average_property (project_to_Q1, projection_matrix, expansion_matrix,

diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc
@@ -203,16 +203,20 @@ namespace aspect
 
         // Functionality to average the additional RHS terms over the cell is not implemented.
         // Consequently, it is only possible to use elasticity with the Material averaging schemes
-        // 'none', 'harmonic average only viscosity', 'project to Q1 only viscosity'.
+        // 'none', 'harmonic average only viscosity', 'geometric average only viscosity', and
+        // 'project to Q1 only viscosity'.
         AssertThrow((this->get_parameters().material_averaging == MaterialModel::MaterialAveraging::none
                      ||
                      this->get_parameters().material_averaging == MaterialModel::MaterialAveraging::harmonic_average_only_viscosity
                      ||
+                     this->get_parameters().material_averaging == MaterialModel::MaterialAveraging::geometric_average_only_viscosity
+                     ||
                      this->get_parameters().material_averaging == MaterialModel::MaterialAveraging::project_to_Q1_only_viscosity),
                     ExcMessage("Material models with elasticity can only be used with the material "
-                               "averaging schemes 'none', 'harmonic average only viscosity', and "
-                               "project to Q1 only viscosity'. This parameter ('Material averaging') "
-                               "is located within the 'Material model' subsection."));
+                               "averaging schemes 'none', 'harmonic average only viscosity', "
+                               "'geometric average only viscosity', and 'project to Q1 only viscosity'. "
+                               "This parameter ('Material averaging') is located within the 'Material "
+                               "model' subsection."));
       }
 
 

diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
@@ -1315,7 +1315,8 @@ namespace aspect
     AssertThrow((sim.parameters.material_averaging &
                  (avg::arithmetic_average | avg::harmonic_average | avg::geometric_average
                   | avg::pick_largest | avg::project_to_Q1 | avg::log_average
-                  | avg::harmonic_average_only_viscosity | avg::project_to_Q1_only_viscosity)) != 0,
+                  | avg::harmonic_average_only_viscosity | avg::geometric_average_only_viscosity
+                  | avg::project_to_Q1_only_viscosity)) != 0,
                 ExcMessage("The matrix-free Stokes solver currently only works if material model averaging "
                            "is enabled. If no averaging is desired, consider using ``project to Q1 only "
                            "viscosity''."));

diff --git a/tests/geometric_average_only_viscosity.prm b/tests/geometric_average_only_viscosity.prm
@@ -0,0 +1,87 @@
+# This is a setup to test geometric averaging of the viscosity
+# and that it works with the GMG.
+# The viscosity in the model decreases exponentially from left to
+# right from 1e14 Pa s to 1e10 Pa s. The model has two cells, so
+# using the geometric average only for the viscosity should lead
+# to a viscosity of 1e11 in the cells on the left and 1e13 in the
+# cells on the right.
+
+set Dimension                              = 2
+set End time                               = 0
+
+set Adiabatic surface temperature          = 100
+
+subsection Solver parameters
+  subsection Stokes solver parameters
+    set Stokes solver type = block GMG
+  end
+end
+
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 100
+    set Y extent = 100
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Variable names = x,z
+    set Function expression = 200 - x
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = top, bottom
+  set List of model names = initial temperature
+end
+
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+
+subsection Material model
+  set Model name = simple
+  set Material averaging = geometric average only viscosity
+
+  subsection Simple model
+    set Viscosity = 1e14
+    set Thermal viscosity exponent = 9.210340372 #ln(10000)
+    set Reference temperature = 100
+    set Minimum thermal prefactor = 1e-6
+  end
+end
+
+
+subsection Mesh refinement
+  set Initial global refinement                = 1
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set Output format = gnuplot
+    set List of output variables      = material properties
+    set Interpolate output            = false
+  end
+end
diff --git a/tests/geometric_average_only_viscosity/screen-output b/tests/geometric_average_only_viscosity/screen-output
@@ -0,0 +1,16 @@
+
+Vectorization over 2 doubles = 128 bits (SSE2), VECTORIZATION_LEVEL=1
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 84 (50+9+25)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving Stokes system... 8+0 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-geometric_average_only_viscosity/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/geometric_average_only_viscosity/solution/solution-00000.0000.gnuplot b/tests/geometric_average_only_viscosity/solution/solution-00000.0000.gnuplot
@@ -0,0 +1,35 @@
+# This file was generated by the deal.II library.
+# Date =  2022/1/28
+# Time =  15:15:37
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <velocity> <velocity> <p> <T> <density> <thermal_expansivity> <specific_heat> <viscosity> 
+0 0 0 0 3.29517e+06 200 3293.4 2e-05 1250 1e+11 
+50 0 -0.0521887 0 3.2953e+06 150 3296.7 2e-05 1250 1e+11 
+
+0 50 0 0.556488 1.64788e+06 200 3293.4 2e-05 1250 1e+11 
+50 50 7.79423e-16 -0.0222499 1.64788e+06 150 3296.7 2e-05 1250 1e+11 
+
+
+50 0 -0.0521887 0 3.2953e+06 150 3296.7 2e-05 1250 1e+13 
+100 0 0 0 3.29726e+06 100 3300 2e-05 1250 1e+13 
+
+50 50 7.79423e-16 -0.0222499 1.64788e+06 150 3296.7 2e-05 1250 1e+13 
+100 50 0 -0.0447992 1.64788e+06 100 3300 2e-05 1250 1e+13 
+
+
+0 50 0 0.556488 1.64788e+06 200 3293.4 2e-05 1250 1e+11 
+50 50 7.79423e-16 -0.0222499 1.64788e+06 150 3296.7 2e-05 1250 1e+11 
+
+0 100 0 0 585.307 200 3293.4 2e-05 1250 1e+11 
+50 100 0.0521887 0 460.499 150 3296.7 2e-05 1250 1e+11 
+
+
+50 50 7.79423e-16 -0.0222499 1.64788e+06 150 3296.7 2e-05 1250 1e+13 
+100 50 0 -0.0447992 1.64788e+06 100 3300 2e-05 1250 1e+13 
+
+50 100 0.0521887 0 460.499 150 3296.7 2e-05 1250 1e+13 
+100 100 0 0 -1506.3 100 3300 2e-05 1250 1e+13 
+
+
diff --git a/tests/geometric_average_only_viscosity/statistics b/tests/geometric_average_only_viscosity/statistics
@@ -0,0 +1,12 @@
+# 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: Iterations for temperature solver
+# 8: Iterations for Stokes solver
+# 9: Velocity iterations in Stokes preconditioner
+# 10: Schur complement iterations in Stokes preconditioner
+# 11: Visualization file name
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 0 7 9 9 output-geometric_average_only_viscosity/solution/solution-00000