Add TODOs

Update TODOs based on discussion and reviews

Update TODO

Update TODOs

Add reaction_rates

Update according to App E

Add ReactionRates

Update fill_reaction_terms

Update fill_elastic_force_outputs

Update parse_parameters

Update calculate_isostrain_viscosities

Assert no plastic damper

Update comments

Fix 2 bugs

Fix bug and indent

Use same averaging moduli as viscosity

Update comments

Fix bug reaction rate

Update stress visu pp for new vep formulation

Update benchmark prm

Add pp to output the min and max of the stress components

Get and use elastic viscosity

Address comments

Fix reaction terms and update some comments

Fix old evaluation of composition

Add yielding to plastic output

Add stress residual pp

Use evaluator for old composition solution

Make composition evaluator only for indep tensor comp

Use solution and reaction term for particle stress

Use max composition for plastic output angles

Fix precision statistiscs

Add reaction rates to VE MM

Address comments

Update benchmarks

Update stress postprocessors

Indent

Add DG for composition to all prms

Assert discontinuous composition discretization

Use iterated Advection and (Newton) Stokes in all prms

Use operator splitting in all prms

Also assert for viscoelastic material model

Add VE relaxation benchmark

Add VE relaxation benchmark python plotting script

Also plot error VE relaxation benchmark

Add particle version VE relaxation benchmark

Update counter and add labels

Add more runs, markers and labels

Part 1 stress averaging

Part 2 stress averaging

Use correctly scaled viscosities and avoid timestep_ratio of zero at t0

Only output statistics for the current stress

Add function to get elastic timestep

Fix timestep in pp

Use correct timestep in t0

Clean up and indent

Also update viscoelastic material model for stress averaging

Simplify order ve viscosity calculation

Add includes lost in rebase

Indent

Update stress pp for dtc isnot dte

Clean up after review 1

Address review comments 2

Use stored particle property instead of interpolated solution

Rm viscoelastic from stress residual pp

Also estimate timestep size before initialization

Set mask old stress ve mm

Fix typo

Apply initial residual fix

Apply stress reaction to correct particle property

Subtract old stress instead of current linearization point in reaction term

Add missing include

Clean up

Indent

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

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

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

benchmarks/viscoelastic_plastic_shear_bands/gerya_2019/gerya_2019_vep.prm

@@ -6,11 +6,28 @@ 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

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

benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm

@@ -29,11 +29,22 @@ set Maximum time step                          = 0.01
 set Output directory                           = output-viscoelastic_plastic_simple_shear
 set Pressure normalization                     = surface
 set Surface pressure                           = 0.
-set Nonlinear solver scheme                    = single Advection, iterated Stokes
+set Nonlinear solver scheme                    = iterated Advection and Stokes
 set Max nonlinear iterations                   = 30
 set Nonlinear solver tolerance                 = 1e-5
 set Max nonlinear iterations in pre-refinement = 0
 
+# To update the elastic stresses,
+# we do 1 operator splitting step.
+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
@@ -42,6 +53,11 @@ subsection Solver parameters
   end
 end
 
+subsection Discretization
+  # DG for viscoelastic stresses
+  set Use discontinuous composition discretization = true
+end
+
 subsection Geometry model
   set Model name = box
   subsection Box

benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm

@@ -24,12 +24,24 @@ set Dimension                              = 2
 set Start time                             = 0
 set End time                               = 1e3
 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                      = 5
 set Output directory                       = output-viscoelastic_plate_flexure
 set Pressure normalization                 = no
 
+# To update the elastic stresses,
+# we do 1 operator splitting step.
+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 (500 m initial resolution)
 subsection Geometry model
   set Model name = box

benchmarks/viscoelastic_relaxation/stress_xx_over_time.py

@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc("pdf", fonttype=42)
+rc("lines", linewidth=5, markersize=3)
+
+# Change file name modifiers as needed 
+# This script assumes both the timestep (500, 250, 125 yr)
+# and the uniform resolution (0, 1, 2) have been varied.
+names = [
+         "ve_relaxation_dt500yr_dh10km",
+         "ve_relaxation_dt500yr_dh5km",
+         "ve_relaxation_dt500yr_dh2-5km",
+         "ve_relaxation_dt250yr_dh10km",
+         "ve_relaxation_dt250yr_dh5km",
+         "ve_relaxation_dt250yr_dh2-5km",
+         "ve_relaxation_dt125yr_dh10km",
+         "ve_relaxation_dt125yr_dh5km",
+         "ve_relaxation_dt125yr_dh2-5km"
+        ]
+tail = r"/statistics"
+
+# The labels the graphs will get in the plot
+labels = [
+          'dt = 500 yr, dh = 10 km',
+          'dt = 500 yr, dh = 5 km',
+          'dt = 500 yr, dh = 2.5 km',
+          'dt = 250 yr, dh = 10 km',
+          'dt = 250 yr, dh = 5 km',
+          'dt = 250 yr, dh = 2.5 km',
+          'dt = 125 yr, dh = 10 km',
+          'dt = 125 yr, dh = 5 km',
+          'dt = 125 yr, dh = 2.5 km'
+         ]
+# Set the colors available for plotting
+color1=[0.0051932, 0.098238, 0.34984]
+color2=[0.092304, 0.32922, 0.38504]
+color3=[0.32701, 0.4579, 0.28638]
+color4=[0.67824, 0.55071, 0.1778]
+color5=[0.97584, 0.63801, 0.50183]
+color6=[0.98447, 0.78462, 0.93553]
+colors = [color1, color3, color5, color1, color3, color5, color1, color3, color5]
+# Set the colors available for plotting
+linestyles = ['solid', 'solid', 'solid', 'solid', 'solid', 'solid', 'solid', 'solid', 'solid'] 
+linestyles = ['solid', 'dashed', 'dotted', 'solid', 'dashed', 'dotted', 'solid', 'dashed', 'dotted'] 
+linestyles = ['solid', 'solid', 'solid', 'dashed', 'dashed', 'dashed', 'dotted', 'dotted', 'dotted'] 
+markers = ['o', '|', '', 'o', '|', '', 'o', '|', ''] 
+markers = ['', '', '', '', '', '', '', '', ''] 
+
+# Set up a row of two plots, one with absolute stress values
+# and one with the percentage difference to the analytical solution
+fig = plt.figure(figsize=(10, 6))
+ax = [fig.add_subplot(2, 1, i) for i in range(1, 3)]
+
+yr_in_secs = 3600. * 24. * 365.2425
+counter = 0 
+
+# Create file path
+for name in names: 
+  path = name+tail
+
+  # Read in the time and the minimum xx and yy components of the viscoelastic stress,
+  # which are stored on the fields ve_stress_xx and ve_stress_yy.
+  # The correct columns are selected with usecols.
+  time,stress_xx_min = np.genfromtxt(path, comments='#', usecols=(1,15), unpack=True)
+
+  # Plot the stress elements in MPa against time in ky in
+  # categorical batlow colors.
+  ax[0].plot(time/1e3,stress_xx_min/1e6,label=labels[counter],color=colors[counter],linestyle=linestyles[counter],marker=markers[counter])
+  ax[1].plot(time/1e3,(stress_xx_min-(20e6*np.exp(-1e10*time*yr_in_secs/1e22)))/(20e6*np.exp(-1e10*time*yr_in_secs/1e22))*100.,label=labels[counter],color=colors[counter],linestyle=linestyles[counter],marker=markers[counter])
+
+  counter += 1
+
+# Plot the analytical solution
+# tau_xx(t) = tau_xx_t0 * exp(-mu*t/eta_viscous), 
+# with tau_xx_t0 = 20 MPa, eta_viscous = 1e22 Pas, mu = 1e10 Pa.
+ax[0].plot(time/1e3,20*np.exp(-1e10*time*yr_in_secs/1e22),label='analytical',color='black',linestyle='dashed')
+
+# Labelling of plot
+ax[1].set_xlabel("Time [ky]")
+ax[0].set_ylabel(r"Viscoelastic stress $\tau_{xx}$ [MPa]")
+ax[1].set_ylabel(r"Error [%]")
+# Manually place legend in lower right corner. 
+ax[0].legend(loc='upper right')
+#ax[1].legend(loc='lower right')
+# Grid and tickes
+ax[0].grid(axis='x',color='0.95')
+ax[0].set_yticks([0,5,10,15,20])
+ax[0].grid(axis='y',color='0.95')
+ax[1].grid(axis='x',color='0.95')
+ax[1].grid(axis='y',color='0.95')
+ax[1].set_yticks([0,1,2,3,4,5,6,7])
+
+# Ranges of the axes
+ax[0].set_xlim(0,250) # kyr
+ax[0].set_ylim(0,21) # MPa
+ax[1].set_xlim(0,250) # kyr
+ax[1].set_ylim(0,7) # %
+
+# Add labels a) and b)
+ax[0].text(-13,20,"a)")
+ax[1].text(-13,6.8,"b)")
+
+# Add timestep labels
+ax[1].text(150,4.1,"dt = 500 yr", rotation = 13)
+ax[1].text(150,2.2,"dt = 250 yr", rotation = 6)
+ax[1].text(150,0.4,"dt = 125 yr", rotation = 3)
+
+
+plt.tight_layout()
+
+# Save as pdf
+plt.savefig('1_viscoelastic_relaxation.pdf')