Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add viscoelastic material model and elastic outputs functionality #1593

Merged
merged 20 commits into from May 25, 2018
Merged

Add viscoelastic material model and elastic outputs functionality #1593

Show file tree
Hide file tree
Changes from 19 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
c72f050
Add viscoelastic material model and elasticity interface
naliboff May 10, 2017
87311fd
Add viscoelastic material model tests
naliboff Feb 7, 2018
a6ca397
Remove grad_phi_u variables
naliboff Feb 25, 2018
2c1c698
Update viscoelastic tests
naliboff Feb 27, 2018
393a8c2
Update viscoelastic documentation and formatting requirements
naliboff Mar 20, 2018
f5cc347
Restructure and cleanup viscoelastic implementation
naliboff Mar 20, 2018
ca59260
Add viscoelastic additional outputs
naliboff Mar 22, 2018
39bb765
Check for access to strain rate
naliboff Mar 22, 2018
37d6e2c
Fix assert error messages
naliboff Mar 22, 2018
a56428f
Add viscoelastic stress build-up benchmark
naliboff Mar 27, 2018
3523e6a
Add viscoelastic bending beam benchmark
naliboff Mar 27, 2018
5401f9b
Move enable_elasticity to Formulation subsection
naliboff Apr 5, 2018
5f00c16
Fix method for obtaining old velocity gradients in viscoelastic mater…
naliboff Apr 5, 2018
ebae58b
Address comments
naliboff May 2, 2018
b5371cc
Add ElasticOutputs for Stokes RHS
naliboff May 3, 2018
dc3669e
Address comments
naliboff May 15, 2018
4cc9f61
Minor cleanup of viscoelastic material model
naliboff May 21, 2018
1a804fd
Update viscoelastic tests
naliboff May 21, 2018
06aa82b
Update viscoelastic bending beam benchmark
naliboff May 21, 2018
1c4781e
Minor updates to viscoelastic material model
naliboff May 24, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
184 changes: 184 additions & 0 deletions benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm
View file
Open in desktop
@@ -0,0 +1,184 @@
# This parameter file reproduces a 2D benchmark for viscoelastic deformation
# published in "Numerical modelling of magma dynamics coupled to tectonic
# deformation of lithosphere and crust" by Keller, May and Kaus, 2013,
# Geophys. J. Int., v. 195, p. 1406-1422. The model setup and results are
# described in detail in section B.2.3 and figure B5.
#
# The benchmark examines bending and unbending (e.g., recovery) of a
# viscoelastic beam surrounded by a less dense and viscous (inelastic) fluid.
# Gravitational forces drive initial bending (elastic strain) of the beam for
# 50 Kyr, which then recovers its shape over ~ 500 Kyr if gravity is turned
# off. The recovery is driven by the accumulated elastic stresses and thus
# provides a basic test for the viscoelasticity implementation.
#
# Compositional fields are used to track the viscoelastic stresses and material
# representing the beam. To improve the accuracy of tracking the beam interface,
# a discontinuous discretization and bound preserving limiter is used for the
# compositional fields. Significantly, the time step is limited to 1 Kyr through
# the "Maximum time step" parameter and a relatively high (0.5) CFL value. Using
# a constant time step ensures the effective (viscoelastic) viscosity of the
# viscoelastic beam and viscous fluid remains constants throughout the model run.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does that mean the solution is time step dependent otherwise? Is there a separate elastic timestep from the convection timestep?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is an option in viscoelastic.cc to use a separate elastic time step from the convection time step (I call it 'numerical' time step). However, in the documentation it is strongly recommended that one uses the averaging scheme option to account for the difference between the elastic and numerical time step.

The ideal case is to set the elastic time step equal to the numerical time step, which I have done in this benchmark. For simplicity and testing purposes, I simply restricted the numerical time step so the effective (viscoelastic) viscosities are constant through time.

One only really runs into trouble if the convection time step becomes much smaller than realistic elastic time steps (~ 10^3-10^5 years). Unfortunately, this can happen when using a really fine mesh.

This situation is effectively the opposite of what the operating splitting was designed for (reaction time scale <<<<< numerical time step). Maybe others will have some insight on how to bridge these time scales.

Most time steps in lithospheric deformation simulations are on the order of 10^3-10^4 years, so this issue does not occur too often.

#
# As currently constructed, the model will run for 50 Kyr with a gravitational
# acceleration of 10 m/s^2. To produce unbending of the beam after the model
# has finished, change the parameter "End time" to 500 Kyr and set the
# gravitational acceleration to 0 m/s^2. A restart file is written every 10 Kyr
# and the parameter "Resume computation = auto" specifies that a model should
# restart from a checkpoint file if one is present.

# Global parameters
set Dimension = 2
set Start time = 0
set End time = 50e3
set Use years in output instead of seconds = true
set Resume computation = auto
set CFL number = 0.5
set Maximum time step = 1e3
set Output directory = output
set Pressure normalization = surface
set Surface pressure = 0.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you remove those parameters from the list that are set to default values and not very interesting for this problem? I am thinking of use direct solver, linear solver tolerance, nonlinear solver scheme, nonlinear solver tolerance, timing output frequency.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍


# Solver settings
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 2000
end
end

# Model geometry (7.5x5 km, 0.025 km spacing)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you update the comment, this should be 0.1 km, correct? (7.5 km, and 75 repetitions?)

subsection Geometry model
set Model name = box

subsection Box
set X repetitions = 75
set Y repetitions = 50
set X extent = 7.5e3
set Y extent = 5e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end

# Element types
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 1
set Use locally conservative discretization = false
set Use discontinuous temperature discretization = false
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
end
end

# Formulation classification
subsection Formulation
set Enable elasticity = true
end

# Velocity boundary conditions
subsection Boundary velocity model
set Zero velocity boundary indicators = left
set Tangential velocity boundary indicators = bottom, top, right
end

# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 4
set Names of fields = stress_xx, stress_yy, stress_xy, beam
end

# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
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)
end
end

# Composition boundary conditions
subsection Boundary composition model
set Fixed composition boundary indicators = bottom, top, right
set List of model names = initial composition
end

# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top, left, right
set List of model names = box
subsection Box
set Bottom temperature = 293
set Left temperature = 293
set Right temperature = 293
set Top temperature = 293
end
end

# Temperature initial conditions
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 293
end
end

# Material model
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 Fixed elastic time step = 1e3
set Use fixed elastic time step = false
set Use stress averaging = false
set Viscosity averaging scheme = maximum composition
end

end

# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.
end
end

# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
subsection Visualization
set List of output variables = material properties, strain rate

subsection Material properties
set List of material properties = density, viscosity
end

set Time between graphical output = 0
set Interpolate output = true
end

end

# Termination criteria
subsection Termination criteria
set Termination criteria = end step
set End step = 500
end

subsection Checkpointing
set Steps between checkpoint = 10
end
192 changes: 192 additions & 0 deletions benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm
View file
Open in desktop
@@ -0,0 +1,192 @@
# This parameter file reproduces a 2D analytical benchmark for the build-up
# of stress in an initially unstressed viscoelastic medium subject to a
# constant strain-rate. This benchmark is described in "Robust
# characteristics method for modelling multiphase visco-elasto-plastic
# thermo-mechanical problems" by Gerya and Yuen, 2007, Phys. Earth. Planet.
# Inter., v. 163, p. 83-105. Full details of the benchmark are located in
# section 3.1 and figure 3 of this manuscript.
#
# The analytical solution for viscoelastic stress build-up in
# an incompressible medium with a constant strain-rate is:
# simga_ij = 2 * edot_ii * eta * (1 - e^(-mu*t/eta)),
# where sigma_ij is the elastic deviatoric stress tensor, edot_ij is
# the deviatoric strain-rate, eta is the background viscosity, mu is the
# shear modulus and t is time.
#
# Following the conditions described in section 3.1 and figure 3 from
# Gerya and Yuen (2007), a 100x100 km body is subject to a constant
# strain-rate of 1.e-14 s^-1 in both the horizontal and vertical directions.
# Constant deformation is driven by inflow and outflow, respectively,
# on the right and bottom walls. The top and left walls are free-slip.
# The material has a viscosity of 1e22 Pa s and a shear modulus of 1e10 Pa.
# With these values, the analytical solution predicts a horizontal
# or vertical viscoelastic stress magnitude of ~ 200 MPa after 250 Kyr
# of deformation. Significantly, the effective "viscoelastic" viscosity
# is enforced to be constant through time by using a constant time step.
# This is achieved by setting a maximum time step (1000 years) much lower
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

much lower than

# than the time step size given by the CFL number of 0.5.
#
# This result can be observed by viewing the compositional field values
# representing the horizontal (stress_xx) or vertical (stress_yy) components
# of the viscoselastic stress tensor. Significantly, the stress near the
# model boundaries is incorrect due to the compositional field boundary
# conditions, which are based on the initial compositional values (e.g., zero).
# This leads to oscillations in the stress field near the boundaries, which
# decay towards a constant value in the model interior. This phenomen highlights
# one advantage to tracking viscoelastic stress with particles, which do not require
# defining compositional boundary conditions.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

well the fields do not need BCs either. Have you tried removing the boundary indicators from Fixed composition boundary indicators? Then there is no flux of stress across the boundary. I am not sure what that will create in the model though.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't tried removing the composition boundary indicators yet, but will give it a shot. However, I assumed initially that having no boundary conditions at all would give odd results.


# Global parameters
set Dimension = 2
set Start time = 0
set End time = 250e3
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, single Stokes
set Nonlinear solver tolerance = 1e-5
set Max nonlinear iterations = 1
set CFL number = 0.5
set Maximum time step = 1000
set Output directory = output
set Timing output frequency = 1
set Pressure normalization = surface
set Surface pressure = 0.

# Solver settings
subsection Solver parameters
subsection Stokes solver parameters
set Use direct solver for Stokes system = false
set Linear solver tolerance = 1e-7
set Number of cheap Stokes solver steps = 2000
end
end

# Model geometry (100x100 km, 2 km spacing)
subsection Geometry model
set Model name = box

subsection Box
set X repetitions = 50
set Y repetitions = 50
set X extent = 100e3
set Y extent = 100e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end

# Element types
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 1
end

# Formulation classification
subsection Formulation
set Enable elasticity = true
end

# Velocity boundary conditions
subsection Boundary velocity model
set Tangential velocity boundary indicators = top, left
set Prescribed velocity boundary indicators = bottom y:function, right x:function
subsection Function
set Variable names = x,y
set Function constants = cm=0.01, year=1, vel=3.154
set Function expression = if (x>50e3 , vel*cm/year, 0.); if (y<50e3 , vel*cm/year, 0.);
end
end

# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 3
set Names of fields = stress_xx, stress_yy, stress_xy
end

# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants =
set Function expression = 0; 0; 0;
end
end

# Composition boundary conditions
subsection Boundary composition model
set Fixed composition boundary indicators = bottom, top, left, right
set List of model names = initial composition
end

# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top, left, right
set List of model names = box
subsection Box
set Bottom temperature = 273
set Left temperature = 273
set Right temperature = 273
set Top temperature = 273
end
end

# Temperature initial conditions
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 273
end
end

# Material model
subsection Material model

set Model name = viscoelastic

subsection Viscoelastic
set Densities = 2800
set Viscosities = 1.e22
set Elastic shear moduli = 1.e10
set Use fixed elastic time step = false
set Fixed elastic time step = 1e3
set Use stress averaging = false
set Viscosity averaging scheme = harmonic
end

end

# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 0.
end
end

# Post processing
subsection Postprocess
set List of postprocessors = basic statistics, composition statistics, temperature statistics, velocity statistics, visualization

subsection Visualization
set List of output variables = material properties, strain rate

subsection Material properties
set List of material properties = density, viscosity
end

set Time between graphical output = 10e3
set Interpolate output = true
end

end

# Termination criteria
subsection Termination criteria
set Termination criteria = end time
end