Update VEP implementation

geodynamics · Mar 12, 2024 · 0d916b3 · 0d916b3
1 parent 845dd6d
commit 0d916b3
Show file tree

Hide file tree

Showing 70 changed files with 2,775 additions and 469 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
@@ -84,8 +97,9 @@ 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
+  set Types of fields  = stress, stress, stress, stress, stress, stress
 end
 
 # Spatial domain of different compositional fields
@@ -95,7 +109,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
 
@@ -129,7 +143,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
 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
@@ -85,8 +98,10 @@ 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
+  set Types of fields  = stress, stress, stress, stress, stress, stress, stress, stress, stress, stress, stress, stress
 end
 
 # Spatial domain of different compositional fields
@@ -96,7 +111,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
 
@@ -130,7 +145,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
 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
@@ -73,8 +84,8 @@ subsection Discretization
 
   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
 
@@ -126,8 +137,9 @@ 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
+  set Types of fields  = stress, stress, stress, stress, stress, stress, chemical composition
 end
 
 #Spatial domain of different compositional fields
@@ -137,7 +149,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
 
@@ -173,7 +185,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.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
@@ -75,8 +88,8 @@ subsection Discretization
 
   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
 
@@ -93,8 +106,9 @@ 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
+  set Types of fields  = stress, stress, stress, stress, stress, stress, chemical composition
 end
 
 # Spatial domain of different compositional fields
@@ -104,7 +118,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
 
@@ -141,9 +155,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, 3300
+    set Viscosities          = 1.e18, 1.e24
+    set Elastic shear moduli = 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_bending_beam/viscoelastic_bending_beam_particles.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam_particles.prm
@@ -6,12 +6,16 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_bending_beam/viscoelastic_ben
 
 set Output directory = output_viscoelastic_bending_beam_particles
 
+# On particles, the operator splitting is handled by the particle property elastic stress.
+set Use operator splitting                 = false
+
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 4
-  set Compositional field methods =           particles,           particles,           particles,         particles
-  set Names of fields             =           ve_stress_xx,        ve_stress_yy,        ve_stress_xy,      beam
-  set Mapped particle properties  = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy, beam:initial beam
+  set Number of fields = 7
+  set Compositional field methods =           particles,           particles,           particles,         particles,           particles,           particles,         particles
+  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
+  set Types of fields             =           stress,              stress,              stress,            stress,              stress,              stress,            chemical composition
+  set Mapped particle properties  = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy, ve_stress_xx_old:ve_stress_xx_old, ve_stress_yy_old:ve_stress_yy_old, ve_stress_xy_old:ve_stress_xy_old, beam:initial beam
 end
 
 # Post processing

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
@@ -7,15 +7,33 @@ 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
+  set Types of fields  = stress, stress, stress, stress, stress, stress, chemical composition, chemical composition, chemical composition
 end
 
 # Spatial domain of different compositional fields
@@ -25,7 +43,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);
@@ -38,9 +56,7 @@ subsection Material model
   set Model name = visco plastic
 
   subsection Visco Plastic
-    # 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, 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.
@@ -49,9 +65,8 @@ subsection Material model
     set Fixed elastic time step     = 20
     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.,    37.,    0.,    0.
+    set Cohesions                   = 100.e6, 100.e6, 10.e6, 10.e6
   end
 end