Update VEP implementation

anne-glerum · Jul 13, 2023 · 5e859c0 · 5e859c0
1 parent 3c58091
commit 5e859c0
Show file tree

Hide file tree

Showing 42 changed files with 2,238 additions and 288 deletions.
diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm
@@ -9,7 +9,7 @@ set Dimension                              = 2
 set Start time                             = 0
 set End time                               = 1500
 set Use years in output instead of seconds = true
-set Nonlinear solver scheme                = single Advection, single Stokes
+set Nonlinear solver scheme                = iterated Advection and Stokes
 set CFL number                             = 0.5
 set Maximum time step                      = 10
 set Output directory                       = output_free_surface_VE_cylinder_2D_loading
@@ -38,11 +38,24 @@ subsection Mesh refinement
   set Strategy         = strain rate
 end
 
+# One operator splitting step to update the stresses
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
+
 # Element types
 subsection Discretization
   set Composition polynomial degree     = 2
   set Stokes velocity polynomial degree = 2
   set Temperature polynomial degree     = 1
+  # DG for viscoelastic stresses
+  set Use discontinuous composition discretization = true
 end
 
 # Formulation classification
@@ -83,8 +96,8 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 3
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Number of fields = 6
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
 end
 
 # Spatial domain of different compositional fields
@@ -93,7 +106,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0;
+    set Function expression = 0; 0; 0; 0; 0; 0;
   end
 end
 
@@ -127,7 +140,7 @@ subsection Material model
     set Viscosities                 = 3.e20
     set Elastic shear moduli        = 1.e10
     set Use fixed elastic time step = false
-    set Fixed elastic time step     = 1e2
+    set Fixed elastic time step     = 10
     set Viscosity averaging scheme  = harmonic
   end
 

diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm
@@ -10,7 +10,7 @@ set Dimension                              = 3
 set Start time                             = 0
 set End time                               = 1500
 set Use years in output instead of seconds = true
-set Nonlinear solver scheme                = single Advection, single Stokes
+set Nonlinear solver scheme                = iterated Advection and Stokes
 set CFL number                             = 0.5
 set Maximum time step                      = 10
 set Output directory                       = output_free_surface_VE_cylinder_3D_loading
@@ -39,11 +39,24 @@ subsection Mesh refinement
   set Strategy         = strain rate
 end
 
+# One operator splitting step to update the stresses
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
+
 # Element types
 subsection Discretization
   set Composition polynomial degree     = 2
   set Stokes velocity polynomial degree = 2
   set Temperature polynomial degree     = 1
+  # DG for viscoelastic stresses
+  set Use discontinuous composition discretization = true
 end
 
 # Formulation classification
@@ -84,8 +97,9 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz
+  set Number of fields = 12
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz, \\
+                         ve_stress_xx_old, ve_stress_yy_old, ve_stress_zz_old, ve_stress_xy_old, ve_stress_xz_old, ve_stress_yz_old
 end
 
 # Spatial domain of different compositional fields
@@ -94,7 +108,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y,z
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0;
   end
 end
 
@@ -128,7 +142,7 @@ subsection Material model
     set Viscosities                 = 3.e20
     set Elastic shear moduli        = 1.e10
     set Use fixed elastic time step = false
-    set Fixed elastic time step     = 1e2
+    set Fixed elastic time step     = 10
     set Viscosity averaging scheme  = harmonic
   end
 

diff --git a/benchmarks/infill_density/infill_ascii.prm b/benchmarks/infill_density/infill_ascii.prm
@@ -23,12 +23,23 @@ set Dimension                              = 2
 set Start time                             = 0
 set End time                               = 120e3
 set Use years in output instead of seconds = true
-set Nonlinear solver scheme                = single Advection, single Stokes
+set Nonlinear solver scheme                = iterated Advection and Stokes
 set CFL number                             = 0.1
 set Maximum time step                      = 750
 set Output directory                       = output
 set Pressure normalization                 = no
 
+# One operator splitting step to update the stresses
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
+
 # Model geometry (1200x1200 km, 80 km spacing in x direction, 40km spacing in y direction)
 subsection Geometry model
   set Model name = box
@@ -70,8 +81,8 @@ subsection Discretization
   set Use discontinuous composition discretization = true
   subsection Stabilization parameters
       set Use limiter for discontinuous composition solution = true
-      set Global composition maximum =  1.e13,  1.e13,  1.e13, 1.0
-      set Global composition minimum = -1.e13, -1.e13, -1.e13, 0.0
+      set Global composition maximum =  1.e13,  1.e13,  1.e13,  1.e13,  1.e13,  1.e13, 1.0
+      set Global composition minimum = -1.e13, -1.e13, -1.e13, -1.e13, -1.e13, -1.e13, 0.0
   end
 end
 # Formulation classification
@@ -118,8 +129,8 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 4
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, lithosphere
+  set Number of fields = 7
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old, lithosphere
 end
 
 #Spatial domain of different compositional fields
@@ -128,7 +139,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; if(y>=1125e3, 1, 0)
+    set Function expression = 0; 0; 0; 0; 0; 0; if(y>=1125e3, 1, 0)
   end
 end
 
@@ -161,7 +172,7 @@ subsection Material model
     # fields representing viscoelastic stresses, but significantly these values are not used in the
     # viscosity calculation (i.e., only values representing lithologies are used). Respectively, the
     # fourth and fifth values are assigned to the compositional fields representing the lithosphere.
-    set Viscosities                 = 1.e21, 1.e21, 1.e21, 1.e21, 1.e26
+    set Viscosities                 = 1.e21, 1.e21, 1.e21, 1.e21, 1.e21, 1.e21, 1.e21, 1.e26
     set Elastic shear moduli        = 1.e10
     set Use fixed elastic time step = true
     set Fixed elastic time step     = 7500

diff --git a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm
@@ -37,6 +37,19 @@ set Maximum time step                      = 1e3
 set Output directory                       = output
 set Pressure normalization                 = surface
 set Surface pressure                       = 0.
+set Nonlinear solver scheme                = iterated Advection and Stokes
+set Nonlinear solver tolerance             = 1e-5
+set Max nonlinear iterations               = 100
+
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
 
 # Solver settings
 subsection Solver parameters
@@ -74,8 +87,8 @@ subsection Discretization
   set Use discontinuous composition discretization = true
   subsection Stabilization parameters
       set Use limiter for discontinuous composition solution = true
-      set Global composition maximum =  1.e11,  1.e11,  1.e11, 1.0
-      set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0
+      set Global composition maximum =  1.e11,  1.e11,  1.e11,  1.e11,  1.e11,  1.e11, 1.0
+      set Global composition minimum = -1.e11, -1.e11, -1.e11, -1.e11, -1.e11, -1.e11, 0.0
   end
 end
 
@@ -92,8 +105,8 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 4
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, beam
+  set Number of fields = 7
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old, beam
 end
 
 # Spatial domain of different compositional fields
@@ -102,7 +115,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; if (x<=4.5e3 && y>=2.5e3 && y<=3.0e3, 1, 0)
+    set Function expression = 0; 0; 0; 0; 0; 0; if (x<=4.5e3 && y>=2.5e3 && y<=3.0e3, 1, 0)
   end
 end
 
@@ -138,9 +151,9 @@ subsection Material model
   set Model name = viscoelastic
 
   subsection Viscoelastic
-    set Densities            =  2800,  2800,  2800,  2800,  3300
-    set Viscosities          = 1.e18, 1.e18, 1.e18, 1.e18, 1.e24
-    set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e10
+    set Densities            =  2800,  2800,  2800,  2800,  2800,  2800,  2800,  3300
+    set Viscosities          = 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e24
+    set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e10
     set Fixed elastic time step     = 1e3
     set Use fixed elastic time step = false
     set Viscosity averaging scheme  = maximum composition

diff --git a/benchmarks/viscoelastic_plastic_shear_bands/gerya_2019/gerya_2019_vep.prm b/benchmarks/viscoelastic_plastic_shear_bands/gerya_2019/gerya_2019_vep.prm
@@ -6,15 +6,32 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_plastic_shear_bands/gerya_201
 #  Global parameters
 set End time                               = 500
 set Output directory                       = output_gerya_2019_vep
+set Nonlinear solver scheme                = iterated Advection and Stokes
+
+# One operator splitting step to update the stresses
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
 
 subsection Formulation
   set Enable elasticity = true
 end
 
+subsection Discretization
+  # DG for viscoelastic stresses
+  set Use discontinuous composition discretization = true
+end
+
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, block, air, inclusion
+  set Number of fields = 9
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old, block, air, inclusion
 end
 
 # Spatial domain of different compositional fields
@@ -23,7 +40,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; \
+    set Function expression = 0; 0; 0; 0; 0; 0; \
                               if ( ( x<43.75e3 && y>25.e3 && y<75.e3) || (x>56.25e3 && y>25.e3 && y<75.e3) || (y>56.25e3 && y<75.e3 && x>=43.75e3 && x<=56.25e3) || (y<43.75e3 && y>25.e3 && x>=43.75e3 && x<=56.25e3), 1, 0); \
                               if (y<=25.e3 || y>=75.e3, 1, 0); \
                               if (y<=56.25e3 && y>=43.75e3 && x>=43.75e3 && x<=56.25e3, 1, 0);
@@ -41,7 +58,7 @@ subsection Material model
 
     # For identification purposes, nonsensical Prefactor values (1e-50) are assigned to the compositional
     # fields tracking elastic stresses, which are not included in the volume fractions computation.
-    set Prefactors for dislocation creep          = 5e-24, 1e-50, 1e-50, 1e-50, 5e-24, 5.e-18, 5e-18
+    set Prefactors for dislocation creep          = 5e-24, 1e-50, 1e-50, 1e-50, 1e-50, 1e-50, 1e-50, 5e-24, 5.e-18, 5e-18
     set Stress exponents for dislocation creep    = 1.0
     set Activation energies for dislocation creep = 0.
     set Activation volumes for dislocation creep  = 0.
@@ -52,8 +69,8 @@ subsection Material model
     set Viscosity averaging scheme  = harmonic
 
     # Nonsensical cohesion values (1e50) are assigned to the compositional fields tracking elastic stresses
-    set Angles of internal friction = 37.,       0.,    0.,    0.,    37.,    0.,    0.
-    set Cohesions                   = 100.e6, 1.e50, 1.e50, 1.e50, 100.e6, 10.e6, 10.e6
+    set Angles of internal friction = 37.,       0.,    0.,    0.,    0.,    0.,    0.,    37.,    0.,    0.
+    set Cohesions                   = 100.e6, 1.e50, 1.e50, 1.e50, 1.e50, 1.e50, 1.e50, 100.e6, 10.e6, 10.e6
 
   end
 

diff --git a/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm b/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm
@@ -22,7 +22,7 @@ set Dimension                              = 2
 set Start time                             = 0
 set End time                               = 20e3
 set Use years in output instead of seconds = true
-set Nonlinear solver scheme                = single Advection, iterated Stokes
+set Nonlinear solver scheme                = iterated Advection and Stokes
 set Nonlinear solver tolerance             = 1e-4
 set Max nonlinear iterations               = 100
 set CFL number                             = 0.5
@@ -31,6 +31,17 @@ set Output directory                       = output_kaus_2010_extension
 set Timing output frequency                = 1
 set Pressure normalization                 = no
 
+# One operator splitting step to update the stresses
+set Use operator splitting                 = true
+
+subsection Solver parameters
+  # Make sure to do only 1 splitting step
+  subsection Operator splitting parameters
+    set Reaction time step                     = 5000 # larger than maximum Stokes time step
+    set Reaction time steps per advection step = 1
+  end
+end
+
 # Solver settings
 subsection Solver parameters
   subsection Stokes solver parameters
@@ -64,6 +75,8 @@ subsection Discretization
   set Composition polynomial degree     = 2
   set Stokes velocity polynomial degree = 2
   set Temperature polynomial degree     = 1
+  # DG for viscoelastic stresses
+  set Use discontinuous composition discretization = true
 end
 
 # Formulation classification
@@ -92,8 +105,8 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 4
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, inclusion
+  set Number of fields = 7
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old, inclusion
 end
 
 # Spatial domain of different compositional fields
@@ -102,7 +115,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; \
+    set Function expression = 0; 0; 0; 0; 0; 0; \
                               if (y<=0.4e3 && x>=19.6e3 && x<=20.4e3, 1, 0);
   end
 end
@@ -147,18 +160,18 @@ subsection Material model
     set Maximum viscosity           = 1.e25
     set Minimum viscosity           = 1.e20
 
-    set Prefactors for dislocation creep          = 5e-26, 5e-26, 5e-26, 5e-26, 5e-21
+    set Prefactors for dislocation creep          = 5e-26, 5e-26, 5e-26, 5e-26, 5e-26, 5e-26, 5e-26, 5e-21
     set Stress exponents for dislocation creep    = 1.0
     set Activation energies for dislocation creep = 0.
     set Activation volumes for dislocation creep  = 0.
 
-    set Elastic shear moduli        = 5.e10,5.e10,5.e10,5.e10,1.e50
+    set Elastic shear moduli        = 5.e10,5.e10,5.e10,5.e10,5.e10,5.e10,5.e10,1.e50
     set Use fixed elastic time step = false
     set Fixed elastic time step     = 1e3
     set Viscosity averaging scheme  = maximum composition
 
     set Angles of internal friction = 30.
-    set Cohesions                   = 40.e6,40.e6,40.e6,40.e6,1.e20
+    set Cohesions                   = 40.e6,40.e6,40.e6,40.e6,40.e6,40.e6,40.e6,1.e20
 
   end