address comments

Co-Authored-By: naliboff <john.naliboff@gmail.com>

benchmarks/free_surface_tractions/viscoelastic/README

@@ -11,7 +11,9 @@ The benchmark compares their solution for the surface displacement
 
 A surface pressure of rho_l*g*H0 (where rho_l is the load density,
  H0 is the load height) is applied on the surface for r<r0 (where
- r0 is the load radius), for t>0. The load is fully removed by t=t1,
+ r0 is the load radius), for t>0. The load is either instantaneously
+ loaded/emplaced at times t0<t<t1, or thinned linearly over this
+ time interval ("linear unloading"). The load is fully removed by t=t1,
  and the surface relaxes. This is done both in a 2-D and 3-D geometry
  (by symmetry the load is centered on the left boundary or left/front
  corner). The input files are:
@@ -22,7 +24,7 @@ A surface pressure of rho_l*g*H0 (where rho_l is the load density,
 
 Changing the stress instantaneously as the load is applied/removed
  leads to oscillations in the elastic response. This can be addressed
- by setting a fixed elastic timestep larger than the numerical 
+ by setting a fixed elastic timestep equal to the numerical ("Maximum")
  timestep and turning on stress averaging:
    'free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm'
    'free_surface_VE_cylinder_3D_loading_fixed_elastic_dt.prm'
@@ -36,4 +38,19 @@ The 'topography' output files may be compared against the analytical
  deflection through time, 'compare_viscous_def_profile.gnuplot' for
  deflection  of profile through time). Equivalent scripts for the
  unloading case are also provided.
+
+Note that while the analytical and numerical results for the deflection
+ of the surface agree well near the center of the load (left boundary),
+ the solutions do not match as well on the right (free-slip) boundary.
+ The numerical free surface rides up vertically against the right
+ boundary to conserve mass, while the analytical solution assumes an
+ infinite half-space, predicting near-zero displacement far from the
+ load center. This is also true for the viscous benchmark. This problem
+ is less pronounced in 3-D as the extra mass may be distributed over
+ a larger area. An open far (right/back) boundary resolves this problem
+ in 3-D.
 
+The solutions match well in 3-D. In 2-D, the geometries of the loading
+ function are different (Cartesian in ASPECT vs cylindrical analytically).
+ As such, the agreement in 2-D breaks down for small r0 (load width) or
+ if the right boundary is placed further away.

benchmarks/free_surface_tractions/viscoelastic/compare_VE_def.gnuplot

@@ -6,6 +6,6 @@ set title "Depression of surface in response to loading (inst. from 0 to 1000 ye
 set xlabel "Time [years]"
 set ylabel "Maximum depression of surface [m]"
 
-plot 'output_free_surface_VE_cylinder_2D_loading/statistics' using 2:27 title 'Numerical, no stress averaging' with linesp
-replot 'output_free_surface_VE_cylinder_2D_loading_fixed_elastic_dt/statistics' using 2:27 title 'Numerical, with stress averaging' with linesp
+plot 'output_free_surface_VE_cylinder_2D_loading/statistics' using 2:28 title 'Numerical, no stress averaging' with linesp
+replot 'output_free_surface_VE_cylinder_2D_loading_fixed_elastic_dt/statistics' using 2:28 title 'Numerical, with stress averaging' with linesp
 replot 'soln_z0.txt' using 1:2 title 'Analytical (Nakiboglu and Lambeck, 1982)'

benchmarks/free_surface_tractions/viscoelastic/compare_VE_def_unloading.gnuplot

@@ -6,5 +6,5 @@ set title "Depression of surface in response to linear unloading (from 0 to 1000
 set xlabel "Time [years]"
 set ylabel "Maximum depression of surface [m]"
 
-plot 'output_free_surface_VE_cylinder_2D_loading_unloading/statistics' using 2:27 title 'Numerical, no stress averaging' with linesp
+plot 'output_free_surface_VE_cylinder_2D_loading_unloading/statistics' using 2:28 title 'Numerical, no stress averaging' with linesp
 replot 'soln2_z0.txt' using 1:2 title 'Analytical (Nakiboglu and Lambeck, 1982)'

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm

@@ -35,7 +35,7 @@ subsection Mesh refinement
   set Initial adaptive refinement        = 1
   set Initial global refinement          = 1
   set Time steps between mesh refinement = 1
-  set Strategy				 = strain rate
+  set Strategy         = strain rate
 end
 
 # Element types
@@ -50,10 +50,13 @@ subsection Formulation
   set Enable elasticity = true
 end
 
-subsection Free surface
-  set Free surface boundary indicators = top
-  set Surface velocity projection      = normal
+subsection Mesh deformation
   set Additional tangential mesh velocity boundary indicators = left,right
+  set Mesh deformation boundary indicators = top: free surface
+
+  subsection Free surface
+    set Surface velocity projection      = normal
+  end
 end
 
 # Velocity boundary conditions
@@ -72,8 +75,8 @@ subsection Boundary traction model
     #set Function expression = 0; 1.e9 + -20.e3*x
     set Function constants  = r0=100.e3, H0=1.e3, t1=1.e3, rhoi=900, g=9.8, t0=1.e3
         # r0 is load radius, H0 is load height, t1 is time load is fully removed,
-	# rhoi is density of ice/load
-	# option to linearly thin load beginning at time t0.
+  # rhoi is density of ice/load
+  # option to linearly thin load beginning at time t0.
     set Function expression = 0; if (x<r0 ,if(t<t0,-g*H0*rhoi,if(t<t1,-g*H0*rhoi*(1-(t-t0)/t1),0)), 0)
   end
 end

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_fixed_elastic_dt.prm

@@ -1,6 +1,6 @@
-# Same as free_surface_VE_cylinder_2D_loading.prm, instead
-# setting fixed elastic timestep and allowing stress averaging,
-# to smooth out elastic response to instantaneous stress changes.
+# Same as free_surface_VE_cylinder_2D_loading.prm, but uses
+# a fixed elastic timestep in combination with stress averaging
+# to smooth out the elastic response to instantaneous stress changes.
 
 include free_surface_VE_cylinder_2D_loading.prm
 
@@ -21,4 +21,3 @@ subsection Material model
   end
 
 end
-

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading_unloading.prm

@@ -19,8 +19,8 @@ subsection Boundary traction model
     set Variable names = x,y,t
     set Function constants  = r0=100.e3, H0=1.e3, t1=1.e3, rhoi=900, g=9.8, t0=5.e2
         # r0 is load radius, H0 is load height, t1 is time load is fully removed,
-	# rhoi is density of ice/load
-	# option to linearly thin load beginning at time t0.
+  # rhoi is density of ice/load
+  # option to linearly thin load beginning at time t0.
     set Function expression = 0; if (x<r0 ,if(t<t0,-g*H0*rhoi,if(t<t1,-g*H0*rhoi*(1-(t-t0)/t1),0)), 0)
   end
 end

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm

@@ -36,7 +36,7 @@ subsection Mesh refinement
   set Initial adaptive refinement        = 1
   set Initial global refinement          = 1
   set Time steps between mesh refinement = 1
-  set Strategy				 = strain rate
+  set Strategy         = strain rate
 end
 
 # Element types
@@ -51,10 +51,13 @@ subsection Formulation
   set Enable elasticity = true
 end
 
-subsection Free surface
-  set Free surface boundary indicators = top
-  set Surface velocity projection      = normal
-  set Additional tangential mesh velocity boundary indicators = left,right,front,back
+subsection Mesh deformation
+  set Additional tangential mesh velocity boundary indicators = left, right, front, back
+  set Mesh deformation boundary indicators = top: free surface
+
+  subsection Free surface
+    set Surface velocity projection      = normal
+  end
 end
 
 # Velocity boundary conditions
@@ -73,8 +76,8 @@ subsection Boundary traction model
     #set Function expression = 0; 1.e9 + -20.e3*x
     set Function constants  = r0=100.e3, H0=1.e3, t1=1.e3, rhoi=900, g=9.8, t0=1.e3
         # r0 is load radius, H0 is load height, t1 is time load is fully removed,
-	# rhoi is density of ice/load
-	# option to linearly thin load beginning at time t0.
+  # rhoi is density of ice/load
+  # option to linearly thin load beginning at time t0.
     set Function expression = 0; 0; if ((x^2+y^2)<r0^2 ,if(t<t0,-g*H0*rhoi,if(t<t1,-g*H0*rhoi*(1-(t-t0)/t1),0)), 0)
   end
 end

benchmarks/free_surface_tractions/viscous/README

@@ -3,16 +3,19 @@ Benchmark of infinite viscous half-space loaded/unloaded by an
  The Motion of a Viscous Fluid Under a Surface Load (1935), Physics
  6, 265; doi:10.1063/1.1745329
 
+This benchmark is identical to the viscoelastic benchmarks, in the
+ limit of negligible elasticity.
+
 The benchmark compares his solution for the surface displacement
  following the application of an instantaneous cylindrical surface
  load (Eq. 3.34) with the deformation of the free surface in ASPECT
  after imposing surface boundary tractions.
 
 A surface pressure of rho_l*g*H0 (where rho_l is the load density,
- H0 is the load height) is applied on the surface for r<r0 (where
- r0 is the load radius), for t>0. This is done both in a 2-D and 
- 3-D geometry (by symmetry the load is centered on the left boundary
- or left/front corner). The input files are: 
+ H0 is the load height) is applied instantaneously on the surface
+ for r<r0 (where r0 is the load radius), for t>0. This is done both
+ in a 2-D and 3-D geometry (by symmetry the load is centered on the
+ left boundary or left/front corner). The input files are: 
    'free_surface_viscous_cylinder_2D_loading.prm'
    'free_surface_viscous_cylinder_3D_loading.prm'
 
@@ -32,4 +35,20 @@ The 'topography' output files may be compared against an analytical
  ASPECT 'topography' output files, using the provided gnuplot script
  ('compare_viscous_def.gnuplot' for maximum surface deflection
  through time, 'compare_viscous_def_profile.gnuplot' for deflection
- of profile through time). 
+ of profile through time).
+
+Note that while the analytical and numerical results for the deflection
+ of the surface agree well near the center of the load (left boundary),
+ the solutions do not match as well on the right (free-slip) boundary.
+ The numerical free surface rides up vertically against the right
+ boundary to conserve mass, while the analytical solution assumes an
+ infinite half-space, predicting near-zero displacement far from the
+ load center. This is also true for the visoelastic benchmark. This problem
+ is less pronounced in 3-D as the extra mass may be distributed over
+ a larger area. An open far (right/back) boundary resolves this problem
+ in 3-D.
+
+The solutions match well in 3-D. In 2-D, the geometries of the loading
+ function are different (Cartesian in ASPECT vs cylindrical analytically).
+ As such, the agreement in 2-D breaks down for small r0 (load width) or
+ if the right boundary is placed further away.

benchmarks/free_surface_tractions/viscous/free_surface_viscous_cylinder_2D_loading.prm

@@ -5,7 +5,7 @@
 # axisymmetric cylindrical load over a viscous half-space.
 
 # This is done by using the viscoelastic model and 
-# setting a high shear modulus.
+# setting a high shear modulus, which produces viscous behavior. 
 
 include ../viscoelastic/free_surface_VE_cylinder_2D_loading.prm
 
@@ -22,4 +22,3 @@ subsection Material model
   end
 
 end
-

doc/modules/changes/20200205_clerc

@@ -0,0 +1,4 @@
+New: Benchmark of deformation of free surface by traction boundary 
+conditions, over viscous/viscoelastic halfspace.
+<br>
+(F. Clerc, J. Naliboff, C. Thieulot, D. Douglas, G. Ito, 2020/02/05)