Add documentation on grain size model

[jdannberg]

During the hackathon, we discussed that it might be more accurate to keep the stress constant in the substeps we do for grain size evolution, rather than keeping the strain rate constant. I still think this is true, but only if the grain size reduction small. For large grain size reduction, keeping the stress constant can lead to a runaway effect: We have a large stress, this reduces the grain size by a lot and leads to a lower viscosity. Usually, a lower viscosity would also mean a lower stress, but keeping the stress constant overestimates the grain size reduction. This effect is especially bad at boundaries where we prescribe the velocity (in other words, the strain rate).

I tested this both with the grain_size_strain test and with an application model similar to the grain)size_ridge model, and I am attaching two images here to show what I mean.

For small rates of grain size change, the "new" (constant stress) method is slightly more accurate.

For large rates of grain size change, the "new" (constant stress) method can leads to a completely wrong result.

Therefore, I think it's better to keep the old method, and I am documenting this in the grain_size_change function.

• [x] I have followed the instructions for indenting my code.
• [x] I have tested my new feature locally to ensure it is correct.

For completeness, below is the prm I used for the test case.

# This is a test for the grain size evolutionn in a realistic
# model setup, in this case, the grain_size_ridge cookbook.
# We make the problem more well-posed

include $ASPECT_SOURCE_DIR/cookbooks/grain_size_ridge/grain_size_ridge.prm

set End time                    = 20000
set Nonlinear solver tolerance  = 1e-5
set Maximum time step           = 50
set Output directory            = grain-size-ridge-ode-50

set Maximum first time step                = 1e300
set Maximum relative increase in time step = 10

subsection Solver parameters
  subsection Stokes solver parameters
    set GMRES solver restart length = 200
    set Number of cheap Stokes solver steps = 5000
    set Stokes solver type = block GMG
    set Linear solver tolerance = 1e-8
  end
end

subsection Mesh refinement
  set Initial global refinement = 4
  set Initial adaptive refinement = 2
  set Minimum refinement level = 4
  set Time steps between mesh refinement = 0
  set Strategy = minimum refinement function

  subsection Minimum refinement function
    set Coordinate system  = cartesian
    set Variable names      = x,y
    set Function expression = if(x*x+(y-200000)*(y-200000) < 60000*60000, 6, 4)
    set Function expression = if((y-200000)*(y-200000) < 60000*60000, 6, 4)
  end
end

subsection Geometry model
  set Model name = box
  subsection Box
    set X extent  = 410000
    set Y extent  = 200000
    set X repetitions = 2
    set Y repetitions = 1
  end
end

subsection Boundary velocity model
  set Tangential velocity boundary indicators = left, top
  set Prescribed velocity boundary indicators = bottom x: function, right y:function

  subsection Function
    set Variable names      = x,y
    set Function expression = 0; 0     #0.05*0.5*(1.0+tanh(0.00005*(y-140000))); 0
  end
end

subsection Boundary traction model
  set Prescribed traction boundary indicators = right x:initial lithostatic pressure, bottom y:initial lithostatic pressure

  subsection Initial lithostatic pressure
    set Representative point         = 60000, 0
    set Number of integration points = 10000
  end
end

# Larger conductivity to make up for the smaller model domain.
subsection Material model
  set Material averaging = harmonic average only viscosity
  subsection Grain size model
    set Thermal conductivity  = 20
    set Maximum viscosity     = 1e23
    set Work fraction for boundary area change      = 0.5
  end
end

subsection Initial temperature model
  set List of model names = adiabatic
  subsection Adiabatic
    set Age bottom boundary layer   = 0.0
    set Top boundary layer age model = function

    subsection Age function
      set Function expression = max((x-10000) / 0.05, 0)
    end
  end
end

# Make the grain size constant everywhere.
subsection Initial composition model
  set Model name = function
  subsection Function
    set Function expression = 5e-3
    set Variable names      = x,y
  end
end

subsection Postprocess
  set List of postprocessors = composition statistics, velocity statistics, mass flux statistics, particles

  subsection Particles
    set Number of particles = 500000
    set Minimum particles per cell = 40
    set Maximum particles per cell = 640
    set Time between data output = 1e100
    set Data output format       = none
  end

  subsection Visualization
    set List of output variables      = material properties, nonadiabatic temperature, nonadiabatic pressure, strain rate, named additional outputs, shear stress
  end
end