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 index 6bba7d2ad73..88bcedf4fce 100644 --- 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 index fca02c2d27c..7b6e00b5a16 100644 --- 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 index dbed3116c5c..8c316188acb 100644 --- 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 index 3d14a096e72..ca95861c7af 100644 --- 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 index 8f2fb4c2c68..a5793d32452 100644 --- 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 index ac107463d2e..bfe0923d730 100644 --- 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 diff --git a/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm b/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm index 24bf989f6fb..ca837410698 100644 --- a/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm +++ b/benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm @@ -25,7 +25,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 Newton Stokes +set Nonlinear solver scheme = iterated Advection and Newton Stokes set Nonlinear solver tolerance = 1e-4 set Max nonlinear iterations = 100 set CFL number = 0.5 @@ -34,6 +34,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 @@ -82,6 +93,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 @@ -113,13 +126,16 @@ end # Number and name of compositional fields subsection Compositional fields - set Number of fields = 5 - set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, plastic_strain, inclusion - set Types of fields = stress, stress, stress, strain, chemical composition + set Number of fields = 8 + 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, plastic_strain, inclusion + set Types of fields = stress, stress, stress, stress, stress, stress, strain, chemical composition set Compositional field methods = particles 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, \ plastic_strain: plastic_strain, \ inclusion: initial inclusion end @@ -131,7 +147,7 @@ subsection Initial composition model subsection Function set Variable names = x,y set Function constants = - set Function expression = 0; 0; 0; 0; \ + set Function expression = 0; 0; 0; 0; 0; 0; 0; \ if (y<=0.4e3 && x>=19.6e3 && x<=20.4e3, 1, 0); end end @@ -176,16 +192,19 @@ subsection Material model set Reference strain rate = 1.e-15 set Maximum viscosity = 1.e25 set Minimum viscosity = 1.e20 - set Prefactors for dislocation creep = 5e-26, 5e-26, 5e-26, 5e-26, 5e-26, 5e-21 + + set Prefactors for dislocation creep = 5e-26, 5e-21 set Stress exponents for dislocation creep = 1.0 set Activation energies for dislocation creep = 0. set Activation volumes for dislocation creep = 0. - set Elastic shear moduli = 5.e10, 5.e10, 5.e10, 5.e10, 5.e10, 1.e50 + + set Elastic shear moduli = 5.e10, 1.e50 set Use fixed elastic time step = false set Fixed elastic time step = 1e3 set Viscosity averaging scheme = harmonic set Angles of internal friction = 30. - set Cohesions = 40.e6, 40.e6, 40.e6, 40.e6, 40.e6, 1.e20 + set Cohesions = 40.e6, 1.e20 + set Strain weakening mechanism = plastic weakening with plastic strain only set Start plasticity strain weakening intervals = 0.0 set End plasticity strain weakening intervals = 0.1 diff --git a/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm b/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm index 33285ea8f23..47e6a330a42 100644 --- a/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm +++ b/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm @@ -29,11 +29,23 @@ 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 subsection Stokes solver parameters @@ -41,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 @@ -80,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 # Composition boundary conditions @@ -95,7 +113,7 @@ subsection Initial composition model set Model name = function subsection Function - 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 diff --git a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm b/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm index 8b4f09583d9..f4a545b52c4 100644 --- a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm +++ b/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 @@ -62,8 +74,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, 1.0 - set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0, 0.0 + set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.0, 1.0 + set Global composition minimum = -1.e11, -1.e11, -1.e11, -1.e11, -1.e11, -1.e11, 0.0, 0.0 end end @@ -90,8 +102,9 @@ end # Number and name of compositional fields subsection Compositional fields - set Number of fields = 5 - set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, plate, inclusion + set Number of fields = 8 + 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, plate, inclusion + set Types of fields = stress, stress, stress, stress, stress, stress, chemical composition, chemical composition end # Spatial domain of different compositional fields @@ -101,7 +114,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 (y>12.5e3, 1, 0); \ if (x>=45e3 && y>=7.5e3 && y<=12.5e3, 1, 0); end @@ -133,12 +146,12 @@ subsection Material model set Model name = viscoelastic subsection Viscoelastic - set Densities = 2700, 2700, 2700, 2700, 2700, 1890 - set Viscosities = 1.e17, 1.e17, 1.e17, 1.e17, 1.e35, 1.e35 - set Elastic shear moduli = 1.e50, 1.e50, 1.e50, 1.e50, 3.e10, 3.e10 - set Fixed elastic time step = 5 - set Use fixed elastic time step = false - set Viscosity averaging scheme = maximum composition + set Densities = 2700, 2700, 1890 + set Viscosities = 1.e17, 1.e35, 1.e35 + set Elastic shear moduli = 1.e50, 3.e10, 3.e10 + set Fixed elastic time step = 5 + set Use fixed elastic time step = false + set Viscosity averaging scheme = maximum composition end end diff --git a/benchmarks/viscoelastic_relaxation/stress_xx_over_time.py b/benchmarks/viscoelastic_relaxation/stress_xx_over_time.py new file mode 100644 index 00000000000..7f64dcae3a4 --- /dev/null +++ b/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') diff --git a/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm b/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm new file mode 100644 index 00000000000..93a26fc3bb6 --- /dev/null +++ b/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm @@ -0,0 +1,168 @@ +# This setup models the relaxation of elastic stresses in the absence +# of drivers. +# The analytical solution is: +# \tau_{xx} = \tau_{xx_{t0}} \exp(-t \frac{\mu}{\eta_{viscous}}), +# where \tau_{xx} is the first component of the stored stress tensor, +# \mu is the shear modulus and \eta_{viscous} the viscous viscosity. +# 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 = iterated Advection and Stokes +set Nonlinear solver tolerance = 1e-6 +set Max nonlinear iterations = 20 +set CFL number = 0.5 +set Maximum time step = 125 +set Pressure normalization = surface +set Surface pressure = 0. +set Output directory = ve_relaxation_dt125yr_dh10km + +# New elasticity settings +set Use operator splitting = true + +# Solver settings +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 Geometry model + set Model name = box + + subsection Box + set X repetitions = 10 + set Y repetitions = 10 + 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 Use discontinuous composition discretization = true +end + +# Formulation classification +subsection Formulation + set Enable elasticity = true +end + +# Velocity boundary conditions +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, left, right, bottom +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, 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 +subsection Initial composition model + set Model name = function + subsection Function + set Variable names = x,y + set Function constants = + set Function expression = 20e6; -20e6; 0; 20e6; -20e6; 0 + end +end + +# Composition boundary conditions +# We specify that no boundaries have a fixed composition +# in order to prevent boundary effects from developing +# in the compositional fields tracking viscoelastic stresses. +subsection Boundary composition model + set Fixed composition boundary indicators = +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 + + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 10. + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, velocity statistics, visualization + + subsection Visualization + set List of output variables = material properties, strain rate, named additional outputs, stress + + subsection Material properties + set List of material properties = viscosity + end + + set Time between graphical output = 500 + set Interpolate output = true + set Write higher order output = true + end + +end + + +# Material model +subsection Material model + + set Model name = visco plastic + + subsection Visco Plastic + + # Uniform viscous viscosity of 1e22 Pas + set Viscous flow law = dislocation + set Prefactors for dislocation creep = 5e-23 + set Minimum strain rate = 1e-30 + set Stress exponents for dislocation creep = 1.0 + set Activation energies for dislocation creep = 0. + set Activation volumes for dislocation creep = 0. + + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = false + # Elastic time step for t0 should match the computational timestep + set Fixed elastic time step = 125 + + # Very high cohesions so that plasticity does not play a role + set Angles of internal friction = 0. + set Cohesions = 1e20 + set Maximum yield stress = 1e20 + + set Densities = 2800 + + end + +end diff --git a/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation_particles.prm b/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation_particles.prm new file mode 100644 index 00000000000..3b4023c16f6 --- /dev/null +++ b/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation_particles.prm @@ -0,0 +1,196 @@ +# This setup models the relaxation of elastic stresses in the absence +# of drivers. +# The analytical solution is: +# \tau_{xx} = \tau_{xx_{t0}} \exp(-t \frac{\mu}{\eta_{viscous}}), +# where \tau_{xx} is the first component of the stored stress tensor, +# \mu is the shear modulus and \eta_{viscous} the viscous viscosity. +# 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 = iterated Advection and Stokes +set Nonlinear solver tolerance = 1e-6 +set Max nonlinear iterations = 20 +set CFL number = 0.5 +set Maximum time step = 125 +set Pressure normalization = surface +set Surface pressure = 0. +set Output directory = ve_relaxation_dt125yr_dh10km_particles + +# 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 = 200 + end +end + +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 10 + set Y repetitions = 10 + 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 = 2 + set Use discontinuous composition discretization = true +end + +# Formulation classification +subsection Formulation + set Enable elasticity = true +end + +# Velocity boundary conditions +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, left, right, bottom +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, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old + set Types of fields = stress, stress, stress, stress, stress, stress + set Compositional field methods = particles, particles, particles, particles, particles, particles + 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 +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 = 20e6; -20e6; 0; 20e6; -20e6; 0 + end +end + +# Composition boundary conditions +# We specify that no boundaries have a fixed composition +# in order to prevent boundary effects from developing +# in the compositional fields tracking viscoelastic stresses. +subsection Boundary composition model + set Fixed composition boundary indicators = +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 + + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 10. + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, velocity statistics, visualization, particles + + subsection Visualization + set List of output variables = material properties, strain rate, named additional outputs, stress + + subsection Material properties + set List of material properties = viscosity + end + + set Time between graphical output = 500 + set Interpolate output = true + set Write higher order output = true + end + + subsection Particles + set Time between data output = 500 + set Data output format = vtu + end +end + +subsection Particles + set Minimum particles per cell = 25 + set Maximum particles per cell = 35 + set Particle generator name = random uniform + set Interpolation scheme = cell average + set List of particle properties = elastic stress + set Load balancing strategy = remove and add particles + set Update ghost particles = true + + subsection Generator + subsection Probability density function + set Number of particles = 3000 + end + end +end + + +# Termination criteria +subsection Termination criteria + set Termination criteria = end time +end + + +# Material model +subsection Material model + + set Model name = visco plastic + + subsection Visco Plastic + + # Uniform viscous viscosity of 1e22 Pas + set Viscous flow law = dislocation + set Prefactors for dislocation creep = 5e-23 + set Minimum strain rate = 1e-30 + set Stress exponents for dislocation creep = 1.0 + set Activation energies for dislocation creep = 0. + set Activation volumes for dislocation creep = 0. + + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = false + # Elastic time step for t0 should match the computational timestep + set Fixed elastic time step = 125 + + # Very high cohesions so that plasticity does not play a role + set Angles of internal friction = 0. + set Cohesions = 1e20 + set Maximum yield stress = 1e20 + + set Densities = 2800 + + end + +end diff --git a/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm b/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm index f5c6bf5106f..c2c84662add 100644 --- a/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm +++ b/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm @@ -40,7 +40,7 @@ set Maximum time step = 0.01 set Output directory = output-viscoelastic_sheared_torsion 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 @@ -52,6 +52,22 @@ subsection Solver parameters end end +# 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 + +subsection Discretization + set Use discontinuous composition discretization = true +end + subsection Geometry model set Model name = box @@ -102,8 +118,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 # Composition boundary conditions @@ -117,7 +135,7 @@ subsection Initial composition model set Model name = function subsection Function - set Function expression = 0.0;0.0;0.0;0.0;0.0;0.0 + set Function expression = 0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0 end end diff --git a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm index aa4b2f544da..4cacc3ba429 100644 --- a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm +++ b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm @@ -37,16 +37,27 @@ 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 scheme = iterated Advection and Stokes set Nonlinear solver tolerance = 1e-5 -set Max nonlinear iterations = 1 +set Max nonlinear iterations = 100 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. +# 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 subsection Stokes solver parameters @@ -77,6 +88,7 @@ end # Element types subsection Discretization + set Use discontinuous composition discretization = true set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 1 @@ -101,8 +113,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 @@ -112,7 +125,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 diff --git a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_particles.prm b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_particles.prm index 427f2439829..1709eb69d58 100644 --- a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_particles.prm +++ b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_particles.prm @@ -6,10 +6,14 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_ set Output directory = output_viscoelastic_stress_build-up_particles +# For particles, the operator splitting is done by the particle property elastic stress. +set Use operator splitting = false + # Number and name of compositional fields subsection Compositional fields - set Compositional field methods = particles, particles, particles - set Mapped particle properties = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy + set Compositional field methods = particles, particles, particles, particles, particles, particles + 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 end # Post processing diff --git a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield.prm b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield.prm index a0d946e36fc..f327baf2c0f 100644 --- a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield.prm +++ b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield.prm @@ -17,11 +17,9 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm -set Nonlinear solver scheme = single Advection, iterated Stokes -set Nonlinear solver tolerance = 1e-5 -set Max nonlinear iterations = 100 set End time = 100e3 + # Material model subsection Material model set Model name = visco plastic @@ -41,3 +39,10 @@ subsection Material model set Densities = 2800 end end + +# Post processing +subsection Postprocess + subsection Visualization + set List of output variables = material properties, strain rate, named additional outputs, stress + end +end diff --git a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield_particles.prm b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield_particles.prm index 1d915b02615..36ad3e57389 100644 --- a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield_particles.prm +++ b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up_yield_particles.prm @@ -15,16 +15,17 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm -set Nonlinear solver scheme = single Advection, iterated Stokes -set Nonlinear solver tolerance = 1e-5 -set Max nonlinear iterations = 100 set End time = 100e3 set Output directory = output_viscoelastic_stress_build-up_yield_particles +# For particles, the operator splitting is done by the particle property elastic stress. +set Use operator splitting = false + # Number and name of compositional fields subsection Compositional fields - set Compositional field methods = particles, particles, particles - set Mapped particle properties = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy + set Compositional field methods = particles, particles, particles, particles, particles, particles + 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 end # Material model diff --git a/doc/modules/changes/20231204_glerum b/doc/modules/changes/20231204_glerum new file mode 100644 index 00000000000..ddbc2b6dafc --- /dev/null +++ b/doc/modules/changes/20231204_glerum @@ -0,0 +1,10 @@ + Changed: The implementation of visco-elasticity and + visco-elasto-plasticity has been updated to properly + track stresses over time. This means that iterative + advection schemes need to be used, as well as the DG + method for compositions. For fields, operator splitting + has to be switched on, while in the case of particles, + the particle property 'elastic stress' takes care of + the stress update. +
+ (Anne Glerum, Robert Myhill, Rene Gassmoeller, Juliane Dannberg, John Naliboff, Gerry Puckett, Esther Heckenbach 2023/12/04) diff --git a/include/aspect/heating_model/shear_heating.h b/include/aspect/heating_model/shear_heating.h index ee9cd7bf53c..5ce951688cc 100644 --- a/include/aspect/heating_model/shear_heating.h +++ b/include/aspect/heating_model/shear_heating.h @@ -95,6 +95,9 @@ namespace aspect * Additional output fields for the shear heating computation * to be added to the MaterialModel::MaterialModelOutputs structure * and filled in the MaterialModel::evaluate() function. + * + * This structure allows to modify the shear heating term by + * multiplying it with a factor computed by the material model. */ template class ShearHeatingOutputs : public MaterialModel::NamedAdditionalMaterialOutputs @@ -112,6 +115,31 @@ namespace aspect */ std::vector shear_heating_work_fractions; }; + + /** + * Additional output fields for the shear heating computation + * to be added to the MaterialModel::MaterialModelOutputs structure + * and filled in the MaterialModel::evaluate() function. + * + * This structure allows to prescribe the full shear heating term from + * the material model. + */ + template + class PrescribedShearHeatingOutputs : public MaterialModel::NamedAdditionalMaterialOutputs + { + public: + PrescribedShearHeatingOutputs(const unsigned int n_points); + + std::vector get_nth_output(const unsigned int idx) const override; + + /** + * The viscous dissipation rate contributing to shear heating. + * If this object is created and filled by the material model + * it will replace the default viscous dissipation rate + * computed by the shear heating model. + */ + std::vector prescribed_shear_heating_rates; + }; } } diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h index 0f66ed8ee09..9d400088ddd 100644 --- a/include/aspect/material_model/interface.h +++ b/include/aspect/material_model/interface.h @@ -1178,8 +1178,8 @@ namespace aspect { public: ElasticOutputs(const unsigned int n_points) - : elastic_force(n_points, numbers::signaling_nan>()) - , viscoelastic_strain_rate(n_points, numbers::signaling_nan>()) + : elastic_force(n_points, numbers::signaling_nan>()), + viscoelastic_strain_rate(n_points, numbers::signaling_nan>()) {} ~ElasticOutputs() override diff --git a/include/aspect/material_model/rheology/elasticity.h b/include/aspect/material_model/rheology/elasticity.h index fdb28edb24f..494df2ed955 100644 --- a/include/aspect/material_model/rheology/elasticity.h +++ b/include/aspect/material_model/rheology/elasticity.h @@ -34,9 +34,9 @@ namespace aspect using namespace dealii; /** - * Additional output fields for the elastic shear modulus to be added to - * the MaterialModel::MaterialModelOutputs structure and filled in the - * MaterialModel::Interface::evaluate() function. + * Additional output fields for the elastic shear modulus and other + * elastic outputs to be added to the MaterialModel::MaterialModelOutputs + * structure and filled in the MaterialModel::Interface::evaluate() function. */ template class ElasticAdditionalOutputs : public NamedAdditionalMaterialOutputs @@ -52,6 +52,21 @@ namespace aspect * the current object. */ std::vector elastic_shear_moduli; + + /** + * Elastic viscosity at the evaluation points passed to + * the instance of MaterialModel::Interface::evaluate() that fills + * the current object. + */ + std::vector elastic_viscosity; + + /** + * The deviatoric stress of the current timestep, so including + * the rotation, advection and stress update, at the evaluation points + * passed to the instance of MaterialModel::Interface::evaluate() + * that fills the current object. + */ + std::vector> deviatoric_stress; }; @@ -76,26 +91,40 @@ namespace aspect parse_parameters (ParameterHandler &prm); /** - * Create the additional material model outputs object that contains the - * elastic shear moduli. + * Create the two additional material model output objects that contain the + * elastic shear moduli, elastic viscosity, ratio of computational to elastic timestep, + * and deviatoric stress of the current timestep and the reaction rates. */ void - create_elastic_outputs (MaterialModel::MaterialModelOutputs &out) const; + create_elastic_additional_outputs (MaterialModel::MaterialModelOutputs &out) const; /** * Given the stress of the previous time step in the material model inputs @p in, - * the elastic shear moduli @p average_elastic_shear_moduli a each point, + * the elastic shear moduli @p average_elastic_shear_moduli at each point, * and the (viscous) viscosities given in the material model outputs object @p out, - * fill an additional material model outputs objects with the elastic force terms. + * fill a material model outputs objects with the elastic force terms, viscoelastic + * strain rate and viscous dissipation. */ void fill_elastic_outputs (const MaterialModel::MaterialModelInputs &in, const std::vector &average_elastic_shear_moduli, MaterialModel::MaterialModelOutputs &out) const; + /** + * Given the stress of the previous time step in the material model inputs @p in, + * the elastic shear moduli @p average_elastic_shear_moduli at each point, + * and the (viscous) viscosities given in the material model outputs object @p out, + * fill a material model outputs (ElasticAdditionalOutputs) object with the + * average shear modulus, elastic viscosity, and the deviatoric stress of the current timestep. + */ + void + fill_elastic_additional_outputs (const MaterialModel::MaterialModelInputs &in, + const std::vector &average_elastic_shear_moduli, + MaterialModel::MaterialModelOutputs &out) const; + /** * Given the stress of the previous time step in the material model inputs @p in, - * the elastic shear moduli @p average_elastic_shear_moduli a each point, + * the elastic shear moduli @p average_elastic_shear_moduli at each point, * and the (viscous) viscosities given in the material model outputs object @p out, * compute an update to the elastic stresses and use it to fill the reaction terms * material model output property. @@ -105,6 +134,18 @@ namespace aspect const std::vector &average_elastic_shear_moduli, MaterialModel::MaterialModelOutputs &out) const; + /** + * Given the stress of the previous time step in the material model inputs @p in, + * the elastic shear moduli @p average_elastic_shear_moduli at each point, + * and the (viscous) viscosities given in the material model outputs object @p out, + * compute the update to the elastic stresses of the previous timestep and use it + * to fill the reaction rates material model output property. + */ + void + fill_reaction_rates (const MaterialModel::MaterialModelInputs &in, + const std::vector &average_elastic_shear_moduli, + MaterialModel::MaterialModelOutputs &out) const; + /** * Return the values of the elastic shear moduli for each composition used in the * rheology model. @@ -113,7 +154,7 @@ namespace aspect get_elastic_shear_moduli () const; /** - * Calculates the effective elastic viscosity (this is the equivalent viscosity of + * Calculate the effective elastic viscosity (this is the equivalent viscosity of * a material which was unstressed at the end of the previous timestep). */ double @@ -139,7 +180,9 @@ namespace aspect */ SymmetricTensor<2,dim> calculate_viscoelastic_strain_rate (const SymmetricTensor<2,dim> &strain_rate, - const SymmetricTensor<2,dim> &stored_stress, + const SymmetricTensor<2, dim> &stress_0_advected, + const SymmetricTensor<2, dim> &stress_old, + const double viscosity_pre_yield, const double shear_modulus) const; /** @@ -148,7 +191,23 @@ namespace aspect double elastic_timestep () const; + /** + * Calculate the ratio between the computational timestep and + * the elastic timestep. + */ + double + calculate_timestep_ratio() const; + private: + /** + * Get the stored stress of the previous timestep. For fields, use a + * composition evaluator of the old solution. For particles, get the + * stress directly from the particles, which is available from in.composition. + */ + std::vector> + retrieve_stress_previous_timestep (const MaterialModel::MaterialModelInputs &in, + const std::vector> &quadrature_positions) const; + /** * Viscosity of a damper used to stabilize elasticity. * A value of 0 Pas is equivalent to not using a damper. @@ -169,7 +228,7 @@ namespace aspect bool use_fixed_elastic_time_step; /** - * Double for fixed elastic time step value, read from parameter file + * Double for fixed elastic time step value, read from parameter file. */ double fixed_elastic_time_step; @@ -183,13 +242,15 @@ namespace aspect double stabilization_time_scale_factor; /** - * We cache the evaluator that is necessary to evaluate the old velocity - * gradients. They are required to compute the elastic stresses, but - * are not provided by the material model. - * By caching the evaluator, we can avoid recreating it every time we - * need it. + * We cache the evaluators that are necessary to evaluate the velocity + * gradients and the old compositions. They are required to compute the elastic stresses, + * but are not provided by the material model. + * By caching the evaluators, we can avoid recreating them every time we need them. */ mutable std::unique_ptr> evaluator; + static constexpr unsigned int n_independent_components = SymmetricTensor<2, dim>::n_independent_components; + mutable std::unique_ptr> evaluator_composition; + }; } } diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h index 7b75573ebd2..bb33e1ebd76 100644 --- a/include/aspect/material_model/rheology/visco_plastic.h +++ b/include/aspect/material_model/rheology/visco_plastic.h @@ -72,7 +72,8 @@ namespace aspect std::vector friction_angles; /** - * The plastic yield stress. + * The current plastic yield stress, depending on composition, + * pressure, and strain. */ std::vector yield_stresses; @@ -245,14 +246,16 @@ namespace aspect /** * Enumeration for selecting which type of viscous flow law to use. - * Select between diffusion, dislocation, frank_kamenetskii or composite. + * Select between diffusion, dislocation, frank_kamenetskii, composite, + * or the minimum of the diffusion and dislocation viscosities. */ enum ViscosityScheme { diffusion, dislocation, frank_kamenetskii, - composite + composite, + minimum_diffusion_dislocation } viscous_flow_law; /** diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index cd9666d5fc2..bf1014c9704 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -77,9 +77,12 @@ namespace aspect * representing these components must be named and listed in a very specific * format, which is designed to minimize mislabeling stress tensor components * as distinct 'compositional rock types' (or vice versa). For 2D models, - * the first three compositional fields must be labeled ve_stress_xx, ve_stress_yy - * and ve_stress_xy. In 3D, the first six compositional fields must be labeled - * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz. + * 3+3 consecutive compositional fields must be labeled ve_stress_xx, ve_stress_yy + * ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old. + * In 3D, six consecutive compositional fields must be labeled + * 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. * * Combining this viscoelasticity implementation with non-linear viscous flow * and plasticity produces a constitutive relationship commonly referred to @@ -95,10 +98,11 @@ namespace aspect * * The overview below directly follows Moresi et al. (2003) eqns. 23-32. * However, an important distinction between this material model and - * the studies above is the use of compositional fields, rather than + * the studies above is the option to use compositional fields, rather than * particles, to track individual components of the viscoelastic stress - * tensor. The material model will be updated when an option to track - * and calculate viscoelastic stresses with particles is implemented. + * tensor. Calculating viscoelastic stresses with particles is also implemented, + * and can be switched on by using particles with the particle property + * 'elastic stress'. * Moresi et al. (2003) begins (eqn. 23) by writing the deviatoric * rate of deformation ($\hat{D}$) as the sum of elastic * ($\hat{D_{e}}$) and viscous ($\hat{D_{v}}$) components: @@ -130,12 +134,11 @@ namespace aspect * W^{t}\tau^{t} + \tau^{t}W^{t}$. * In this material model, the size of the time step above ($\Delta t^{e}$) * can be specified as the numerical time step size or an independent fixed time - * step. If the latter case is selected, the user has an option to apply a - * stress averaging scheme to account for the differences between the numerical - * and fixed elastic time step (eqn. 32). If one selects to use a fixed elastic time - * step throughout the model run, an equal numerical and elastic time step can be - * achieved by using CFL and maximum time step values that restrict the numerical - * time step to the fixed elastic time step. + * step. If the latter case is selected, a stress averaging scheme is applied + * to account for the differences between the numerical + * and fixed elastic time step (eqn. 32). A fixed computational time + * step throughout the model run can be achieved by using a large CFL number and + * smaller maximum time step values that restrict the numerical time step to a specific value. * * The formulation above allows rewriting the total rate of deformation (eqn. 29) as * $\tau^{t + \Delta t^{e}} = \eta_{eff} \left ( @@ -150,8 +153,10 @@ namespace aspect * viscosity is reduced relative to the initial viscosity. * * Elastic effects are introduced into the governing Stokes equations through - * an elastic force term (eqn. 30) using stresses from the previous time step: - * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{t}$. + * an elastic force term (eqn. 30 updated to the term in eqn. 5 in Farrington et al. 2014) + * using stresses from the previous time step rotated and advected into the current + * time step: + * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{0adv}$. * This force term is added onto the right-hand side force vector in the * system of equations. * diff --git a/include/aspect/material_model/viscoelastic.h b/include/aspect/material_model/viscoelastic.h index 7003201e0e7..c67cb18862a 100644 --- a/include/aspect/material_model/viscoelastic.h +++ b/include/aspect/material_model/viscoelastic.h @@ -45,9 +45,12 @@ namespace aspect * representing these components must be named and listed in a very specific * format, which is designed to minimize mislabeling stress tensor components * as distinct 'compositional rock types' (or vice versa). For 2D models, - * the first three compositional fields must be labeled ve_stress_xx, ve_stress_yy - * and ve_stress_xy. In 3D, the first six compositional fields must be labeled - * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz. + * 3+3 consecutive compositional fields must be labeled ve_stress_xx, ve_stress_yy, + * ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old. + * In 3D, 6+6 consecutive compositional fields must be labeled + * 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. * * Expanding the model to include non-linear viscous flow (e.g., * diffusion/dislocation creep) and plasticity would produce a constitutive @@ -63,10 +66,11 @@ namespace aspect * * The overview below directly follows Moresi et al. (2003) eqns. 23-32. * However, an important distinction between this material model and - * the studies above is the use of compositional fields, rather than + * the studies above is the option to use compositional fields, rather than * particles, to track individual components of the viscoelastic stress - * tensor. The material model will be updated when an option to track - * and calculate viscoelastic stresses with particles is implemented. + * tensor. Calculating viscoelastic stresses with particles is also implemented, + * and can be switched on by using particles with the particle property + * 'elastic stress'. * * Moresi et al. (2003) begins (eqn. 23) by writing the deviatoric * rate of deformation ($\hat{D}$) as the sum of elastic @@ -92,11 +96,11 @@ namespace aspect * W^{t}\tau^{t} + \tau^{t}W^{t}$. * In this material model, the size of the time step above ($\Delta t^{e}$) * can be specified as the numerical time step size or an independent fixed time - * step. If the latter case is a selected, the user has an option to apply a - * stress averaging scheme to account for the differences between the numerical - * and fixed elastic time step (eqn. 32). If one selects to use a fixed elastic time - * step throughout the model run, this can still be achieved by using CFL and - * maximum time step values that restrict the numerical time step to a specific time. + * step. If the latter case is selected, a stress averaging scheme is applied + * to account for the differences between the numerical + * and fixed elastic time step (eqn. 32). A fixed computational time + * step throughout the model run can be achieved by using a large CFL number and + * smaller maximum time step values that restrict the numerical time step to a specific time. * * The formulation above allows rewriting the total rate of deformation (eqn. 29) as * $\tau^{t + \Delta t^{e}} = \eta_{eff} \left ( @@ -110,9 +114,11 @@ namespace aspect * The magnitude of the shear modulus thus controls how much the effective * viscosity is reduced relative to the initial viscosity. * - * Elastic effects are introduced into the governing stokes equations through - * an elastic force term (eqn. 30) using stresses from the previous time step: - * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{t}$. + * Elastic effects are introduced into the governing Stokes equations through + * an elastic force term (eqn. 30 updated to the term in eqn. 5 in Farrington et al. 2014) + * using stresses from the previous time step rotated and advected into the current + * time step: + * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{0adv}$. * This force term is added onto the right-hand side force vector in the * system of equations. * @@ -128,9 +134,11 @@ namespace aspect * user supplies a comma delimited list of length N+1, where N is the * number of compositional fields. The additional field corresponds to * the value for background material. They should be ordered - * ``background, composition1, composition2...''. However, the first 3 (2D) - * or 6 (3D) composition fields correspond to components of the elastic - * stress tensor and their material values will not contribute to the volume + * ``background, composition1, composition2...''. 3+3 (2D) + * or 6+6 (3D) consecutive composition fields should correspond to components + * of the elastic stress tensor of the previous timestep rotated and advected + * into the current timestep and the tensor of the previous timestep advected + * into the current timestep only. Their material values will not contribute to the volume * fractions. If a single value is given, then all the compositional fields * are given that value. Other lengths of lists are not allowed. For a given * compositional field the material parameters are treated as constant, @@ -176,22 +184,20 @@ namespace aspect /** * Declare the parameters this class takes through input files. */ - static - void - declare_parameters (ParameterHandler &prm); + static void + declare_parameters(ParameterHandler &prm); /** * Read the parameters this class declares from the parameter file. */ void - parse_parameters (ParameterHandler &prm) override; + parse_parameters(ParameterHandler &prm) override; /** * @} */ void - create_additional_named_outputs (MaterialModel::MaterialModelOutputs &out) const override; - + create_additional_named_outputs(MaterialModel::MaterialModelOutputs &out) const override; private: /** diff --git a/include/aspect/particle/property/elastic_stress.h b/include/aspect/particle/property/elastic_stress.h index 9df4f0093df..c9c4388e443 100644 --- a/include/aspect/particle/property/elastic_stress.h +++ b/include/aspect/particle/property/elastic_stress.h @@ -50,6 +50,15 @@ namespace aspect void initialize () override; + /** + * Function to update particles after they have been restored + * to their position and values from the beginning of the timestep. + * This restoring happens at the beginning of nonlinear iterations + * of iterative Advection solver schemes. + */ + void + update_particles (typename Particle::Manager &particle_manager) const; + /** * @copydoc aspect::Particle::Property::Interface::initialize_one_particle_property() */ @@ -94,6 +103,18 @@ namespace aspect */ mutable MaterialModel::MaterialModelInputs material_inputs; mutable MaterialModel::MaterialModelOutputs material_outputs; + + /** + * The indices of the compositional fields that represent components of the + * viscoelastic stress tensors. + */ + std::vector stress_field_indices; + + /** + * The indices of the compositional fields that do not represent components of the + * viscoelastic stress tensors. + */ + std::vector non_stress_field_indices; }; } } diff --git a/include/aspect/postprocess/visualization/stress_residual.h b/include/aspect/postprocess/visualization/stress_residual.h new file mode 100644 index 00000000000..a4abee69c72 --- /dev/null +++ b/include/aspect/postprocess/visualization/stress_residual.h @@ -0,0 +1,65 @@ +/* + Copyright (C) 2011 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#ifndef _aspect_postprocess_visualization_stress_residual_h +#define _aspect_postprocess_visualization_stress_residual_h + +#include +#include + +#include + + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + /** + * A class derived from DataPostprocessor that takes an output vector + * and computes a variable that represents the difference between the + * second moment invariant of the deviatoric stress and the yield stress. + * + * The member functions are all implementations of those declared in the + * base class. See there for their meaning. + */ + template + class StressResidual + : public DataPostprocessorScalar, + public SimulatorAccess, + public Interface + { + public: + /** + * Constructor. + */ + StressResidual (); + + void + evaluate_vector_field(const DataPostprocessorInputs::Vector &input_data, + std::vector> &computed_quantities) const override; + }; + } + } +} + +#endif diff --git a/include/aspect/simulator_signals.h b/include/aspect/simulator_signals.h index 66115523644..12587416234 100644 --- a/include/aspect/simulator_signals.h +++ b/include/aspect/simulator_signals.h @@ -320,6 +320,16 @@ namespace aspect * calling the DataOut member function set_cell_selection(). */ boost::signals2::signal &)> pre_data_out_build_patches; + + /** + * A signal that is triggered after particles have been restored to their position + * and property values from the beginning of the current timestep. This happens + * at the beginning of each nonlinear iteration (except for the first iteration of the timestep) + * of iterative advection schemes. + * Parameters are a reference to a ParticleManager. + */ + boost::signals2::signal &)> post_restore_particles; + }; diff --git a/source/heating_model/shear_heating.cc b/source/heating_model/shear_heating.cc index da8a6debcf0..4ca6bf5b22a 100644 --- a/source/heating_model/shear_heating.cc +++ b/source/heating_model/shear_heating.cc @@ -40,8 +40,20 @@ namespace aspect const ShearHeatingOutputs *shear_heating_out = material_model_outputs.template get_additional_output>(); + const PrescribedShearHeatingOutputs *prescribed_shear_heating_out = + material_model_outputs.template get_additional_output>(); + for (unsigned int q=0; qprescribed_shear_heating_rates[q]; + heating_model_outputs.lhs_latent_heat_terms[q] = 0.0; + continue; + } + const SymmetricTensor<2,dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() ? @@ -72,8 +84,10 @@ namespace aspect heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate; - // If shear heating work fractions are provided, reduce the - // overall heating by this amount (which is assumed to be converted into other forms of energy) + // If shear heating work fractions are provided, multiply the + // overall shear heating by this factor. Work fractions are usually smaller than + // one, and the reduction in shear heating is assumed to be work converted into + // other forms of energy). if (shear_heating_out != nullptr) { heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q]; @@ -176,6 +190,27 @@ namespace aspect return shear_heating_work_fractions; } + + + + template + PrescribedShearHeatingOutputs::PrescribedShearHeatingOutputs (const unsigned int n_points) + : + MaterialModel::NamedAdditionalMaterialOutputs({"prescribed_shear_heating_rates"}), + prescribed_shear_heating_rates(n_points, numbers::signaling_nan()) + {} + + + + template + std::vector + PrescribedShearHeatingOutputs::get_nth_output(const unsigned int idx) const + { + (void) idx; + AssertIndexRange (idx, 1); + + return prescribed_shear_heating_rates; + } } } @@ -194,7 +229,8 @@ namespace aspect "right-hand side of the temperature equation.") #define INSTANTIATE(dim) \ - template class ShearHeatingOutputs; + template class ShearHeatingOutputs; \ + template class PrescribedShearHeatingOutputs; ASPECT_INSTANTIATE(INSTANTIATE) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index b2829f205f0..dc6ef8a4d86 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -21,10 +21,13 @@ #include -#include -#include #include +#include +#include +#include +#include +#include #include @@ -37,7 +40,10 @@ namespace aspect std::vector make_elastic_additional_outputs_names() { std::vector names; + // We do not output all the additional output to + // visualization output, only the ones listed here. names.emplace_back("elastic_shear_modulus"); + names.emplace_back("elastic_viscosity"); return names; } } @@ -46,7 +52,9 @@ namespace aspect ElasticAdditionalOutputs::ElasticAdditionalOutputs (const unsigned int n_points) : NamedAdditionalMaterialOutputs(make_elastic_additional_outputs_names()), - elastic_shear_moduli(n_points, numbers::signaling_nan()) + elastic_shear_moduli(n_points, numbers::signaling_nan()), + elastic_viscosity(n_points, numbers::signaling_nan()), + deviatoric_stress(n_points, SymmetricTensor<2,dim>()) {} @@ -56,7 +64,20 @@ namespace aspect ElasticAdditionalOutputs::get_nth_output(const unsigned int idx) const { (void)idx; // suppress warning in release mode - AssertIndexRange (idx, 1); + AssertIndexRange (idx, 2); + + switch (idx) + { + case 0: + return elastic_shear_moduli; + + case 1: + return elastic_viscosity; + + default: + AssertThrow(false, ExcInternalError()); + } + return elastic_shear_moduli; } @@ -89,8 +110,9 @@ namespace aspect "the model to run the user must select 'true' or 'false'."); prm.declare_entry ("Fixed elastic time step", "1.e3", Patterns::Double (0.), - "The fixed elastic time step $dte$. Units: years if the " - "'Use years in output instead of seconds' parameter is set; " + "The fixed elastic time step $dte$. It is always used during the first " + "timestep; afterwards on if 'Used fixed elastic time step' is true. " + "Units: years if the 'Use years in output instead of seconds' parameter is set; " "seconds otherwise."); prm.declare_entry ("Stabilization time scale factor", "1.", Patterns::Double (1.), @@ -118,6 +140,9 @@ namespace aspect void Elasticity::parse_parameters (ParameterHandler &prm) { + AssertThrow(this->get_parameters().enable_elasticity == true, + ExcMessage("Rheology model elasticity only works if 'Enable elasticity' is set to true")); + // Retrieve the list of composition names std::vector compositional_field_names = this->introspection().get_composition_names(); @@ -154,61 +179,118 @@ namespace aspect if (this->convert_output_to_years()) fixed_elastic_time_step *= year_in_seconds; - AssertThrow(this->get_parameters().enable_elasticity == true, - ExcMessage ("Material model Viscoelastic only works if 'Enable elasticity' is set to true")); - - // Check whether the compositional fields representing the viscoelastic - // stress tensor are both named correctly and listed in the right order. - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_xx") == 0, + // When using the visco_plastic or viscoelastic material model, + // make sure that no damping is applied. Damping could potentially + // improve stability under rapidly changing dynamics, but + // so far it has not been necessary. + AssertThrow(elastic_damper_viscosity == 0., + ExcMessage("The viscoelastic material model and the visco-plastic material model with elasticity enabled require " + "that no elastic damping is applied.")); + // An update of the stored stresses is done in an operator splitting step for fields or by the particle property 'elastic stress'. + AssertThrow(this->get_parameters().use_operator_splitting || (this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx")), + ExcMessage("The viscoelastic material model and the visco-plastic material model with elasticity enabled require " + "operator splitting for stresses tracked on compositional fields or the particle property 'elastic stress' " + "for stresses tracked on particles.")); + if ((this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx"))) + AssertThrow(!this->get_parameters().use_operator_splitting, + ExcMessage("If stresses are tracked on particles, the stress update is applied by the particle property 'elastic stress' " + "and operator splitting should not be turned on. ")); + + // Check that 3+3 in 2D or 6+6 in 3D stress fields exist. + AssertThrow((this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::stress) == 2*SymmetricTensor<2,dim>::n_independent_components), + ExcMessage("Rheology model Elasticity requires 3+3 in 2D or 6+6 in 3D fields of type stress.")); + + // Check that the compositional fields representing the viscoelastic + // stress tensor components are both named correctly and listed in the right order + // as well as use a discontinuous discretization. + std::vector stress_field_names = this->introspection().get_names_for_fields_of_type(CompositionalFieldDescription::stress); + std::vector stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress); + + // The discontinuous element is required to accommodate discontinuous + // strain rates that feed into the stored stresses. + const std::vector use_discontinuous_composition_discretization = this->get_parameters().use_discontinuous_composition_discretization; + for (auto stress_field : stress_field_indices) + AssertThrow(use_discontinuous_composition_discretization[stress_field], + ExcMessage("The viscoelastic material model and the visco-plastic material model with elasticity enabled require " + "the use of discontinuous elements for compositions that represent stress tensor components.")); + + // We require a consecutive range of indices (for example for FEPointEvaluation) + // to extract the fields representing the viscoelastic stress tensor components, + // so check that they are listed without interruption by other fields. + // They do not, however, have to be the first fields listed. + AssertThrow(((stress_field_indices[2*n_independent_components-1] - stress_field_indices[0]) == (2*n_independent_components-1)), + ExcMessage("Rheology model Elasticity requires that the compositional fields representing stress tensor components are listed in consecutive order.")); + + AssertThrow(stress_field_names[0] == "ve_stress_xx", ExcMessage("Rheology model Elasticity only works if the first " - "compositional field is called ve_stress_xx.")); - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_yy") == 1, + "compositional field representing stress tensor components is called ve_stress_xx.")); + AssertThrow(stress_field_names[1] == "ve_stress_yy", ExcMessage("Rheology model Elasticity only works if the second " - "compositional field is called ve_stress_yy.")); + "compositional field representing stress tensor components is called ve_stress_yy.")); if (dim == 2) { - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_xy") == 2, + AssertThrow(stress_field_names[2] == "ve_stress_xy", ExcMessage("Rheology model Elasticity only works if the third " - "compositional field is called ve_stress_xy.")); + "compositional field representing stress tensor components is called ve_stress_xy.")); + + AssertThrow(stress_field_names[3] == "ve_stress_xx_old", + ExcMessage("Rheology model Elasticity only works if the fourth " + "compositional field representing stress tensor components is called ve_stress_xx_old.")); + AssertThrow(stress_field_names[4] == "ve_stress_yy_old", + ExcMessage("Rheology model Elasticity only works if the fifth " + "compositional field representing stress tensor components is called ve_stress_yy_old.")); + AssertThrow(stress_field_names[5] == "ve_stress_xy_old", + ExcMessage("Rheology model Elasticity only works if the sixth " + "compositional field representing stress tensor components is called ve_stress_xy_old.")); } else if (dim == 3) { - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_zz") == 2, + AssertThrow(stress_field_names[2] == "ve_stress_zz", ExcMessage("Rheology model Elasticity only works if the third " - "compositional field is called ve_stress_zz.")); - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_xy") == 3, + "compositional field representing stress tensor components is called ve_stress_zz.")); + AssertThrow(stress_field_names[3] == "ve_stress_xy", ExcMessage("Rheology model Elasticity only works if the fourth " - "compositional field is called ve_stress_xy.")); - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_xz") == 4, + "compositional field representing stress tensor components is called ve_stress_xy.")); + AssertThrow(stress_field_names[4] == "ve_stress_xz", ExcMessage("Rheology model Elasticity only works if the fifth " - "compositional field is called ve_stress_xz.")); - AssertThrow(this->introspection().compositional_index_for_name("ve_stress_yz") == 5, + "compositional field representing stress tensor components is called ve_stress_xz.")); + AssertThrow(stress_field_names[5] == "ve_stress_yz", ExcMessage("Rheology model Elasticity only works if the sixth " - "compositional field is called ve_stress_yz.")); + "compositional field representing stress tensor components is called ve_stress_yz.")); + + AssertThrow(stress_field_names[6] == "ve_stress_xx_old", + ExcMessage("Rheology model Elasticity only works if the seventh " + "compositional field representing stress tensor components is called ve_stress_xx_old.")); + AssertThrow(stress_field_names[7] == "ve_stress_yy_old", + ExcMessage("Rheology model Elasticity only works if the eighth " + "compositional field representing stress tensor components is called ve_stress_yy_old.")); + AssertThrow(stress_field_names[8] == "ve_stress_zz_old", + ExcMessage("Rheology model Elasticity only works if the ninth " + "compositional field representing stress tensor components is called ve_stress_zz_old.")); + AssertThrow(stress_field_names[9] == "ve_stress_xy_old", + ExcMessage("Rheology model Elasticity only works if the tenth " + "compositional field representing stress tensor components is called ve_stress_xy_old.")); + AssertThrow(stress_field_names[10] == "ve_stress_xz_old", + ExcMessage("Rheology model Elasticity only works if the eleventh " + "compositional field representing stress tensor components is called ve_stress_xz_old.")); + AssertThrow(stress_field_names[11] == "ve_stress_yz_old", + ExcMessage("Rheology model Elasticity only works if the twelfth " + "compositional field representing stress tensor components is called ve_stress_yz_old.")); } else AssertThrow(false, ExcNotImplemented()); - // Currently, it only makes sense to use this material model when the nonlinear solver - // scheme does a single Advection iteration and at minimum one Stokes iteration. More - // than one nonlinear Advection iteration will produce an unrealistic build-up of - // viscoelastic stress, which are tracked through compositional fields. + // We need to iterate over the Advection and Stokes equations. AssertThrow((this->get_parameters().nonlinear_solver == - Parameters::NonlinearSolver::single_Advection_single_Stokes - || + Parameters::NonlinearSolver::iterated_Advection_and_Stokes || this->get_parameters().nonlinear_solver == - Parameters::NonlinearSolver::single_Advection_iterated_Stokes - || + Parameters::NonlinearSolver::iterated_Advection_and_Newton_Stokes || this->get_parameters().nonlinear_solver == - Parameters::NonlinearSolver::single_Advection_iterated_Newton_Stokes - || - this->get_parameters().nonlinear_solver == - Parameters::NonlinearSolver::single_Advection_iterated_defect_correction_Stokes), + Parameters::NonlinearSolver::iterated_Advection_and_defect_correction_Stokes), ExcMessage("The material model will only work with the nonlinear " - "solver schemes 'single Advection, single Stokes', " - "'single Advection, iterated Stokes', " - "'single Advection, iterated Newton Stokes', and " - "'single Advection, iterated defect correction Stokes' ")); + "solver schemes 'iterated Advection and Stokes', " + "'iterated Advection and defect correction Stokes', " + "'iterated Advection and Newton Stokes'.")); // Functionality to average the additional RHS terms over the cell is not implemented. // Also, there is no option implemented in this rheology module to project to Q1 the viscosity @@ -233,14 +315,36 @@ namespace aspect template void - Elasticity::create_elastic_outputs (MaterialModel::MaterialModelOutputs &out) const + Elasticity::create_elastic_additional_outputs (MaterialModel::MaterialModelOutputs &out) const { + // Create the ElasticAdditionalOutputs that include the average shear modulus, elastic + // viscosity, timestep ratio and total deviatoric stress of the current timestep. if (out.template get_additional_output>() == nullptr) { const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } + + // We need to modify the shear heating outputs to correctly account for elastic stresses. + if (out.template get_additional_output>() == nullptr) + { + const unsigned int n_points = out.n_evaluation_points(); + out.additional_outputs.push_back( + std::make_unique> (n_points)); + } + + // Create the ReactionRateOutputs that are necessary for the operator splitting + // step (either on the fields or directly on the particles) + // that sets both sets of stresses to the total stress of the + // previous timestep. + if (out.template get_additional_output>() == nullptr && + (this->get_parameters().use_operator_splitting || (this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx")))) + { + const unsigned int n_points = out.n_evaluation_points(); + out.additional_outputs.push_back( + std::make_unique>(n_points, this->n_compositional_fields())); + } } @@ -251,13 +355,20 @@ namespace aspect const std::vector &average_elastic_shear_moduli, MaterialModel::MaterialModelOutputs &out) const { - // Create a reference to the structure for the elastic outputs + // Create a reference to the structure for the elastic outputs. + // The structure is created during the Stokes assembly. MaterialModel::ElasticOutputs *elastic_out = out.template get_additional_output>(); - if (elastic_out == nullptr) + // Create a reference to the structure for the prescribed shear heating outputs. + // The structure is created during the advection assembly. + HeatingModel::PrescribedShearHeatingOutputs + *heating_out = out.template get_additional_output>(); + + if (elastic_out == nullptr && heating_out == nullptr) return; + // TODO should a RHS term be a separate MaterialProperties? if (in.requests_property(MaterialProperties::additional_outputs)) { // The viscosity should be averaged if material averaging is applied. @@ -282,126 +393,296 @@ namespace aspect else effective_creep_viscosities = out.viscosities; + const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx"); + for (unsigned int i=0; i < in.n_evaluation_points(); ++i) { - // Get old stresses from compositional fields - const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][0], - &in.composition[i][0]+SymmetricTensor<2,dim>::n_independent_components)); - - elastic_out->elastic_force[i] = -effective_creep_viscosities[i] / calculate_elastic_viscosity(average_elastic_shear_moduli[i]) * stress_old; - // The viscoelastic strain rate is needed only when the Newton method is selected. - const typename Parameters::NonlinearSolver::Kind nonlinear_solver = this->get_parameters().nonlinear_solver; - if ((nonlinear_solver == Parameters::NonlinearSolver::iterated_Advection_and_Newton_Stokes) || - (nonlinear_solver == Parameters::NonlinearSolver::single_Advection_iterated_Newton_Stokes)) - elastic_out->viscoelastic_strain_rate[i] = calculate_viscoelastic_strain_rate(in.strain_rate[i], stress_old, average_elastic_shear_moduli[i]); + const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]); + + // Get stress from timestep $t$ rotated and advected into the current + // timestep $t+\Delta t_c$ from the compositional fields. + // This function is only evaluated during the assembly of the Stokes equations + // (the force term goes into the rhs of the momentum equation). + // This happens after the advection equations have been solved, and hence in.composition + // contains the rotated and advected stresses $tau^{0adv}$. + // Only at the beginning of the next timestep do we add the stress update of the + // current timestep to the stress stored in the compositional fields, giving + // $\tau{t+\Delta t_c}$ with $t+\Delta t_c$ being the current timestep. + const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], + &in.composition[i][stress_start_index]+n_independent_components)); + + // Get the old stress that is used to interpolate to timestep $t+\Delta t_c$. It is stored on the + // second set of n_independent_components fields, e.g. in 2D on field 3, 4 and 5. + const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], + &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); + + // Average effective creep viscosity + // Use the viscosity corresponding to the stresses selected above. + // out.viscosities is computed during the assembly of the Stokes equations + // based on the current_linearization_point. This means that it will be updated after every + // nonlinear Stokes iteration. + // The effective creep viscosity has already been scaled with the timestep ratio dtc/dte. + const double effective_creep_viscosity = effective_creep_viscosities[i]; + + // The force term is computed as: + // $\frac{-\eta_{effcreep} \tau_{0adv}}{\eta_{e}}$, where $\eta_{effcreep}$ is the + // current harmonic average of the viscous and elastic viscosity, or the yield stress + // divided by two times the second invariant of the deviatoric strain rate. + // In case the computational timestep differs from the elastic timestep, + // linearly interpolate between the two. + const double timestep_ratio = calculate_timestep_ratio(); + // The elastic viscosity has also already been scaled with the timestep ratio. + const double viscosity_ratio = effective_creep_viscosity / calculate_elastic_viscosity(average_elastic_shear_moduli[i]); + + if (elastic_out != nullptr) + { + elastic_out->elastic_force[i] = -1. * (viscosity_ratio * stress_0_advected + + (1. - timestep_ratio) * (1. - viscosity_ratio) * stress_old); + + // The viscoelastic strain rate is needed only when the Newton method is selected. + const typename Parameters::NonlinearSolver::Kind nonlinear_solver = this->get_parameters().nonlinear_solver; + if ((nonlinear_solver == Parameters::NonlinearSolver::iterated_Advection_and_Newton_Stokes) || + (nonlinear_solver == Parameters::NonlinearSolver::single_Advection_iterated_Newton_Stokes)) + elastic_out->viscoelastic_strain_rate[i] = calculate_viscoelastic_strain_rate( + in.strain_rate[i], stress_0_advected, stress_old, effective_creep_viscosity, average_elastic_shear_moduli[i]); + } + + // Apply the stress update to get the total stress of timestep t. + const SymmetricTensor<2, dim> stress = 2. * effective_creep_viscosity * deviatoric_strain_rate + viscosity_ratio * stress_0_advected + + (1. - timestep_ratio) * (1. - viscosity_ratio) * stress_old; + + // Obtain the computational timestep by multiplying the ratio between the computational + // and elastic timestep $\frac{\Delta t_c}{\Delta t_{el}}$ with the elastic timestep. + const double dtc = timestep_ratio * elastic_timestep(); + + // Assume incompressibility. + const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - ((stress - stress_0_advected) / (2. * dtc * average_elastic_shear_moduli[i])); + // If compressible, + // visco_plastic_strain_rate = visco_plastic_strain_rate - + // 1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor(); + + // The shear heating term needs to account for the elastic stress, but only the visco_plastic strain rate. + // This is best computed here, and stored for later use by the heating model. + if (heating_out != nullptr) + heating_out->prescribed_shear_heating_rates[i] = stress * visco_plastic_strain_rate; } + + } + } + + + + + template + void + Elasticity::fill_elastic_additional_outputs (const MaterialModel::MaterialModelInputs &in, + const std::vector &average_elastic_shear_moduli, + MaterialModel::MaterialModelOutputs &out) const + { + // Create a reference to the structure for the elastic additional outputs + MaterialModel::ElasticAdditionalOutputs + *elastic_additional_out = out.template get_additional_output>(); + + if (elastic_additional_out == nullptr || !in.requests_property(MaterialProperties::additional_outputs)) + return; + + const double timestep_ratio = calculate_timestep_ratio(); + const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx"); + + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) + { + const double eta = out.viscosities[i]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]); + const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], + &in.composition[i][stress_start_index]+n_independent_components)); + const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], + &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); + + const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_moduli[i]); + + // Apply the stress update to get the total deviatoric stress of timestep t. + elastic_additional_out->deviatoric_stress[i] = 2. * eta * deviatoric_strain_rate + eta / elastic_viscosity * stress_0_advected + (1. - timestep_ratio) * (1. - eta / elastic_viscosity) * stress_old; + elastic_additional_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; + elastic_additional_out->elastic_viscosity[i] = elastic_viscosity; } } + // Rotate the stresses of the previous timestep $t$ into the current timestep $t+dtc$. template void Elasticity::fill_reaction_outputs (const MaterialModel::MaterialModelInputs &in, - const std::vector &average_elastic_shear_moduli, + const std::vector &, MaterialModel::MaterialModelOutputs &out) const { - if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && in.requests_property(MaterialProperties::reaction_terms)) + if (in.current_cell.state() == IteratorState::valid + && in.requests_property(MaterialProperties::reaction_terms)) { - // Get old (previous time step) velocity gradients + // Get the velocity gradients of the current timestep $t+dtc$ + // at the requested location in in.position. std::vector> quadrature_positions(in.n_evaluation_points()); - for (unsigned int i=0; i < in.n_evaluation_points(); ++i) + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) quadrature_positions[i] = this->get_mapping().transform_real_to_unit_cell(in.current_cell, in.position[i]); + // Get the current velocity gradients, which get + // updated in each nonlinear iteration. + // This means we use the rotation tensor W^(t+dtc), not W^(t). std::vector solution_values(this->get_fe().dofs_per_cell); - in.current_cell->get_dof_values(this->get_old_solution(), + in.current_cell->get_dof_values(this->get_current_linearization_point(), solution_values.begin(), solution_values.end()); - // Only create the evaluator the first time we get here + // Only create the evaluator the first time we get here. if (!evaluator) evaluator = std::make_unique>(this->get_mapping(), this->get_fe(), update_gradients, this->introspection().component_indices.velocities[0]); - // Initialize the evaluator for the old velocity gradients + // Initialize the evaluator for the velocity gradients. evaluator->reinit(in.current_cell, quadrature_positions); evaluator->evaluate(solution_values, EvaluationFlags::gradients); - const double dte = elastic_timestep(); - const double dt = this->get_timestep(); - - // The viscosity should be averaged if material averaging is applied. - // Here the averaging scheme "project to Q1 (only viscosity)" is - // excluded, because there is no way to know the quadrature formula - // used for evaluation. - // TODO: find a way to include "project to Q1 (only viscosity)" as well. - std::vector effective_creep_viscosities; - if (this->get_parameters().material_averaging != MaterialAveraging::none && - this->get_parameters().material_averaging != MaterialAveraging::project_to_Q1 && - this->get_parameters().material_averaging != MaterialAveraging::project_to_Q1_only_viscosity) - { - MaterialModelOutputs out_copy(out.n_evaluation_points(), - this->introspection().n_compositional_fields); - out_copy.viscosities = out.viscosities; + // Get the fully updated stress of the previous timestep $t$. + const std::vector> stress_t = retrieve_stress_previous_timestep(in, quadrature_positions); - const MaterialAveraging::AveragingOperation averaging_operation_for_viscosity = - get_averaging_operation_for_viscosity(this->get_parameters().material_averaging); - MaterialAveraging::average(averaging_operation_for_viscosity, - in.current_cell, - Quadrature(quadrature_positions), - this->get_mapping(), - in.requested_properties, - out_copy); + Assert(out.reaction_terms.size() == in.n_evaluation_points(), ExcMessage("Out reaction terms not equal to n eval points.")); - effective_creep_viscosities = out_copy.viscosities; - } - else - effective_creep_viscosities = out.viscosities; + const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx"); - for (unsigned int i=0; i < in.n_evaluation_points(); ++i) + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) { - // Get old stresses from compositional fields - const SymmetricTensor<2,dim> stress_old(Utilities::Tensors::to_symmetric_tensor(&in.composition[i][0], - &in.composition[i][0]+SymmetricTensor<2,dim>::n_independent_components)); + Assert(out.reaction_terms[i].size() == this->n_compositional_fields(), ExcMessage("Out reaction terms i not equal to n fields.")); - // Calculate the rotated stresses // Rotation (vorticity) tensor (equation 25 in Moresi et al., 2003, J. Comp. Phys.) - const Tensor<2,dim> rotation = 0.5 * (evaluator->get_gradient(i) - transpose(evaluator->get_gradient(i))); - - // Calculate the current (new) stored elastic stress, which is a function of the material - // properties (viscoelastic viscosity, shear modulus), elastic time step size, strain rate, - // vorticity, prior (inherited) viscoelastic stresses and viscosity of the elastic damper. - // In the absence of the elastic damper, the expression for "stress_new" is identical - // to the one found in Moresi et al. (2003, J. Comp. Phys., equation 29). - const double damped_elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_moduli[i]); - - // stress_0 is the combination of the elastic stress tensor stored at the end of the last time step and the change in that stress generated by local rotation - const SymmetricTensor<2,dim> stress_0 = (stress_old + dte * ( symmetrize(rotation * Tensor<2,dim>(stress_old) ) - symmetrize(Tensor<2,dim>(stress_old) * rotation) ) ); + const Tensor<2, dim> rotation = 0.5 * (evaluator->get_gradient(i) - transpose(evaluator->get_gradient(i))); + + // stress_0 (i.e., $\tau^{0}$) is the sum of the stress tensor stored at the end of the last time step (stress_t) + // and the change in that stress generated by local rotation over the computational timestep $\Delta t_c$: + // stress_0 = stress_t + dtc * (symmetrize(rotation * stress_t) - symmetrize(stress_t * rotation)). + // The reaction terms should be filled with the change from the previous stress to the rotated stress_0. + // In case fields are used, in.composition holds the current_linearization_point, + // which for the first nonlinear iteration means the extrapolated solution from + // the last and previous to last timesteps. In later iterations, it holds the current solution. + // If we subtract the stress in in.composition from stress_0 to compute the stress change, we would alternately compute + // a reaction term that is correct and one that is zero. Therefore we subtract the old stress stress_t we retrieved + // above. This means the reaction terms are: + // stress_0 - stress_t = stress_t + dtc * (symmetrize(rotation * stress_t) - symmetrize(stress_t * rotation)) - stress_t + // = dtc * (symmetrize(rotation * stress_t) - symmetrize(stress_t * rotation)). + const SymmetricTensor<2, dim> stress_change = this->get_timestep() * (symmetrize(rotation * Tensor<2, dim>(stress_t[i])) - symmetrize(Tensor<2, dim>(stress_t[i]) * rotation)); + + Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_change, + &out.reaction_terms[i][stress_start_index], + &out.reaction_terms[i][stress_start_index]+n_independent_components); + } + } + } - // stress_creep is the stress experienced by the viscous and elastic components. - Assert(std::isfinite(in.strain_rate[i].norm()), - ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was " - "not filled by the caller.")); - const SymmetricTensor<2,dim> stress_creep = 2. * effective_creep_viscosities[i] * ( deviator(in.strain_rate[i]) + stress_0 / (2. * damped_elastic_viscosity ) ); + // The following function computes the reaction rates for the operator + // splitting step that at the beginning of the new timestep $t+dtc$ updates the + // stored compositions $tau^{0\mathrm{adv}}$ at time $t$ to $tau^{t}$. + // This update consists of the stress change resulting from system evolution, + // but does not advect or rotate the stress tensor. Advection is done by + // solving the advection equation and the stress tensor is rotated through + // the source term (reaction_terms) of that same equation. + template + void + Elasticity::fill_reaction_rates (const MaterialModel::MaterialModelInputs &in, + const std::vector &average_elastic_shear_moduli, + MaterialModel::MaterialModelOutputs &out) const + { + ReactionRateOutputs *reaction_rate_out = out.template get_additional_output>(); - // stress_new is the (new) stored elastic stress - SymmetricTensor<2,dim> stress_new = stress_creep * (1. - (elastic_damper_viscosity / damped_elastic_viscosity)) + elastic_damper_viscosity * stress_0 / damped_elastic_viscosity; + if (reaction_rate_out == nullptr) + return; - // Stress averaging scheme to account for difference between the elastic time step - // and the numerical time step (see equation 32 in Moresi et al., 2003, J. Comp. Phys.) - // Note that if there is no difference between the elastic timestep and the numerical - // timestep, then no averaging occurs as dt/dte = 1. - stress_new = ( ( 1. - ( dt / dte ) ) * stress_old ) + ( ( dt / dte ) * stress_new ) ; + // At the moment when the reaction rates are required (at the beginning of the timestep), + // the solution vector 'solution' holds the stress from the previous timestep, + // advected into the new position of the previous timestep, so $\tau^{t}_{0adv}$. + // This is the same as the vector 'old_solution' holds. At later moments during the current timestep, + // 'solution' will hold the current_linearization_point instead of the solution of the previous timestep. + // + // In case fields are used to track the stresses, MaterialModelInputs are based on 'solution' + // when calling the MaterialModel for the reaction rates. When particles are used, MaterialModelInputs + // for this function are filled with the old_solution (including for the strain rate), except for the + // compositions that represent the stress tensor components, these are taken directly from the + // particles. As the particles are restored to their pre-advection location at the beginning of + // each nonlinear iteration, their values and positions correspond to the old solution. + // This means that in both cases we can use 'in' to get to the $\tau^{t}_{0adv}$ and velocity/strain rate of the + // previous timestep. + if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && in.requests_property(MaterialProperties::reaction_rates)) + { + const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx"); - const SymmetricTensor<2,dim> stress_update = stress_new - stress_old; + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) + { + // Set all reaction rates to zero + for (unsigned int c = 0; c < in.composition[i].size(); ++c) + reaction_rate_out->reaction_rates[i][c] = 0.0; + + // Get $\tau^{0adv}$ of the previous timestep t from the compositional fields. + // This stress includes the rotation and advection of the previous timestep, + // i.e., the reaction term (which prescribes the change in stress due to rotation + // over the previous timestep) has already been applied during the previous timestep. + const SymmetricTensor<2, dim> stress_0_t (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], + &in.composition[i][stress_start_index]+n_independent_components)); + + // Get the old stress that is used to interpolate to timestep $t+\Delta t_c$. It is stored on the + // second set of n_independent_components fields, e.g. in 2D on field 3, 4 and 5. + // The old stress was advected into the previous timestep, but not rotated. + // Below we update it to full stress of the previous timestep, so that it can be + // advected into the current timestep. + const SymmetricTensor<2, dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], + &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); + + // $\eta^{t}_{effcreep}$. This viscosity has been calculated with the timestep_ratio dtc/dte. + const double effective_creep_viscosity = out.viscosities[i]; + + // $\eta_{el} = G \Delta t_c$. + // The elastic viscosity has also already been scaled with $\frac{\Delta t_c} / {\Delta t_{el}}$ + // in light of the linear interpolation between $t$ and $t+ \Delta t_{el}$ + // when $\Delta t_c$ and $t+\Delta t_el$ differ. + const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_moduli[i]); + + // The ratio between the computational and elastic timestep $\frac{\Delta t_c} / {\Delta t_{el}}$. + const double timestep_ratio = calculate_timestep_ratio(); + + // Compute the total stress at time t. + const SymmetricTensor<2, dim> + stress_t = 2. * effective_creep_viscosity * deviator(in.strain_rate[i]) + + effective_creep_viscosity / elastic_viscosity * stress_0_t + + (1. - timestep_ratio) * (1. - effective_creep_viscosity / elastic_viscosity) * stress_old; + + // Fill reaction rates. + // During this timestep, the reaction rates will be multiplied + // with the current timestep size to turn the rate of change into a change. + // However, this update belongs + // to the previous timestep. Therefore we divide by the + // current timestep and multiply with the previous one. + // When multiplied with the current timestep, this will give + // (rate * previous_dt / current_dt) * current_dt = rate * previous_dt = previous_change. + // previous_change = stress_t - stress_0_t. + // To compute the rate we should return to the operator splitting scheme, + // we therefore divide the change in stress by the current timestep current_dt (= dtc). + const double dtc = timestep_ratio * elastic_timestep(); + + const SymmetricTensor<2, dim> stress_update = (stress_t - stress_0_t) / dtc; - // Fill reaction terms Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_update, - &out.reaction_terms[i][0], - &out.reaction_terms[i][0]+SymmetricTensor<2,dim>::n_independent_components); + &reaction_rate_out->reaction_rates[i][stress_start_index], + &reaction_rate_out->reaction_rates[i][stress_start_index]+n_independent_components); + + // Also update the second set of stresses, stress_old, with the newly computed stress, + // which in the rest of the timestep will serve as the old stress advected but not rotated + // into the current timestep. This function fill_reaction_rates is only called at the + // beginning of the timestep, and so this update only happens once. + const SymmetricTensor<2, dim> stress_old_update = (stress_t - stress_old) / dtc; + + Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_old_update, + &reaction_rate_out->reaction_rates[i][stress_start_index+n_independent_components], + &reaction_rate_out->reaction_rates[i][stress_start_index+n_independent_components]+n_independent_components); } } } @@ -412,14 +693,17 @@ namespace aspect double Elasticity::elastic_timestep () const { - // The elastic time step (dte) is equal to the numerical time step if the time step number + // The elastic time step ($\Delta t_el$, dte) is equal to the numerical time step if the time step number // is greater than 0 and the parameter 'use_fixed_elastic_time_step' is set to false. - // On the first (0) time step the elastic time step is always equal to the value + // On the first (0) time step, the elastic time step is always equal to the value // specified in 'fixed_elastic_time_step', which is also used in all subsequent time // steps if 'use_fixed_elastic_time_step' is set to true. // // We also use this parameter when we are still *before* the first time step, // i.e., if the time step number is numbers::invalid_unsigned_int. + if (use_fixed_elastic_time_step && this->get_timestep_number() > 0 && this->simulator_is_past_initialization()) + AssertThrow(fixed_elastic_time_step >= this->get_timestep(), ExcMessage("The elastic timestep has to be equal to or bigger than the numerical timestep")); + const double dte = ( ( this->get_timestep_number() > 0 && this->simulator_is_past_initialization() && use_fixed_elastic_time_step == false ) @@ -432,6 +716,23 @@ namespace aspect + template + double + Elasticity::calculate_timestep_ratio() const + { + // Before the simulator is initialized, get_timestep() can return + // a NaN. During assembly in timestep 0, get_timestep() returns 0. Therefore we guess the + // timestep size using the maximum timestep parameters capped by the elastic timestep. + double dtc = this->get_timestep(); + if (!this->simulator_is_past_initialization() || + (this->get_timestep_number() == 0 && this->get_timestep() == 0)) + dtc = std::min(std::min(this->get_parameters().maximum_time_step, this->get_parameters().maximum_first_time_step), elastic_timestep()); + + return dtc / elastic_timestep(); + } + + + template const std::vector & Elasticity::get_elastic_shear_moduli () const @@ -446,7 +747,17 @@ namespace aspect Elasticity:: calculate_elastic_viscosity (const double shear_modulus) const { - return shear_modulus*elastic_timestep() + elastic_damper_viscosity; + // In case the computational timestep dtc differs from the elastic timestep, + // we need to scale the elastic viscosity with the timestep ratio dtc/dte: + // $\eta_{el} = \Delta t_{el} G$ + // $\eta_{el}^{c} = \Delta t_{el} G \frac{\Delta t_c}{\Delta t_{el}} = \Delta t_c G$. + // Since we already have a function that returns the timestep ratio $\frac{\Delta t_c}{\Delta t_{el}}$ + // and deals with the special cases where $\Delta t_c$ is not know yet, use that function. + // The damper is asserted to be zero when using this rheology at the moment. If the + // damper were to be applied, the user should input the damper viscosity over + // the computational timestep. + const double timestep_ratio = calculate_timestep_ratio(); + return shear_modulus*elastic_timestep()*timestep_ratio + elastic_damper_viscosity; } @@ -457,8 +768,17 @@ namespace aspect calculate_viscoelastic_viscosity (const double viscosity, const double shear_modulus) const { + // The elastic viscosity has been scaled with the timestep ratio $\frac{\Delta t_c}{\Delta t_{el}}$. + // The viscous viscosity has not been scaled yet, so we do it here. + // Scaling both viscosities with the timestep ratio before computing the effective + // viscoelastic viscosity equals scaling the viscoelastic viscosity computed from + // the unscaled elastic and viscous viscosity. + // If the computational timestep is small compared to the elastic timestep, + // the effective viscosity becomes small as well. This reduction is balanced by + // an increasing body force term in the right-hand-side of the momentum equation. + const double timestep_ratio = calculate_timestep_ratio(); const double elastic_viscosity = calculate_elastic_viscosity(shear_modulus); - return 1. / (1./elastic_viscosity + 1./viscosity); + return 1. / (1./elastic_viscosity + 1./(viscosity*timestep_ratio)); } @@ -467,19 +787,107 @@ namespace aspect SymmetricTensor<2,dim> Elasticity:: calculate_viscoelastic_strain_rate(const SymmetricTensor<2,dim> &strain_rate, - const SymmetricTensor<2,dim> &stored_stress, + const SymmetricTensor<2,dim> &stress_0_advected, + const SymmetricTensor<2,dim> &stress_old, + const double viscosity_pre_yield, const double shear_modulus) const { // The first term in the following expression is the deviator of the true strain rate // of one or more isostress rheological elements (in series). // One of these elements must be an elastic component (potentially damped). // The second term corresponds to a fictional strain rate arising from - // elastic stresses stored from the last time step. - // Note the parallels with the viscous part of the strain rate deviator, + // elastic stresses stored from the last time step. Note the parallels with the + // viscous part of the strain rate deviator, // which is equal to 0.5 * stress / viscosity. - return deviator(strain_rate) + 0.5 * deviator(stored_stress) / - calculate_elastic_viscosity(shear_modulus); + + // The elastic viscosity is already scaled with the timestep ratio. + const double elastic_viscosity = calculate_elastic_viscosity(shear_modulus); + // viscosity_pre_yield is also already scaled with the timestep ratio. + const double creep_viscosity = viscosity_pre_yield; + + // The ratio between the computational and elastic timestep. + const double timestep_ratio = calculate_timestep_ratio(); + + const SymmetricTensor<2, dim> + edot_deviator = deviator(strain_rate) + 0.5 * stress_0_advected / elastic_viscosity + + 0.5 * (1. - timestep_ratio) * (1. - creep_viscosity/elastic_viscosity) * stress_old / creep_viscosity; + + return edot_deviator; + } + + + + template + std::vector> + Elasticity:: + retrieve_stress_previous_timestep (const MaterialModel::MaterialModelInputs &in, + const std::vector> &quadrature_positions) const + { + std::vector> stress_t(in.n_evaluation_points(), SymmetricTensor<2, dim>()); + + const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx"); + + // Either get stress_t - the fully updated stress of the previous timestep - from the fields + // or from the particles. + if (this->get_parameters().compositional_field_methods[stress_start_index] == Parameters::AdvectionFieldMethod::fem_field) + { + // Get the compositional fields from the previous timestep $t$. + // The 'old_solution' has been updated to the full stress tensor + // of time $t$ by the operator splitting step at the beginning + // of the current timestep. + std::vector old_solution_values(this->get_fe().dofs_per_cell); + in.current_cell->get_dof_values(this->get_old_solution(), + old_solution_values.begin(), + old_solution_values.end()); + + // Only create the evaluator the first time we get here. + // With the field types, we can avoid specifying the stress + // as the first fields. However, the FEPointEvaluation evaluators + // can select a range of fields, but this has to be a consecutive + // range starting from a given index. To avoid having to evaluate + // all fields, we still request that the stress fields are listed + // in a consecutive order without interruption by other fields. + if (!evaluator_composition) + evaluator_composition.reset(new FEPointEvaluation(this->get_mapping(), + this->get_fe(), + update_values, + this->introspection().component_indices.compositional_fields[stress_start_index])); + + // Initialize the evaluator for the composition values. + evaluator_composition->reinit(in.current_cell, quadrature_positions); + evaluator_composition->evaluate(old_solution_values, + EvaluationFlags::values); + + // Get the composition values representing the viscoelastic stress field tensor components + // of the previous timestep from the evaluator. + // We assume (and assert in parse_parameters) that the 2*n_independent_components + // tensor components are in the correct order and consecutively listed. + // These stresses have not yet been rotated or advected to the current timestep. + std::vector> + stress_t(in.n_evaluation_points(), SymmetricTensor<2, dim>()); + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) + { + const Tensor<1,n_independent_components> composition_values = evaluator_composition->get_value(i); + stress_t[i] = Utilities::Tensors::to_symmetric_tensor(&composition_values[0], + &composition_values[0]+n_independent_components); + } + } + else + { + // If we use particles instead of fields to track the stresses, use the MaterialInputs, + // which include the stresses stored on the particles. + // The computation of the reaction terms (during which this function is called) happens + // after the particles have been restored to their position and values + // at the beginning of the timestep (i.e., the position and values of the previous timestep) and after they + // have been updated to the full stress of the previous timestep by the reaction rates. + for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) + stress_t[i] = Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], + &in.composition[i][stress_start_index]+n_independent_components); + } + + return stress_t; } + } } } diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc index 9a9b4da8522..aadc18a95c5 100644 --- a/source/material_model/rheology/visco_plastic.cc +++ b/source/material_model/rheology/visco_plastic.cc @@ -110,12 +110,32 @@ namespace aspect output_parameters.current_friction_angles.resize(volume_fractions.size(), numbers::signaling_nan()); output_parameters.current_cohesions.resize(volume_fractions.size(), numbers::signaling_nan()); - // Assemble stress tensor if elastic behavior is enabled - SymmetricTensor<2,dim> stress_old = numbers::signaling_nan>(); + // Assemble current and old stress tensor if elastic behavior is enabled + SymmetricTensor<2, dim> stress_0_advected = numbers::signaling_nan>(); + SymmetricTensor<2, dim> stress_old = numbers::signaling_nan>(); + double elastic_shear_modulus = numbers::signaling_nan(); if (this->get_parameters().enable_elasticity) { - for (unsigned int j=0; j < SymmetricTensor<2,dim>::n_independent_components; ++j) - stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j]; + // The first set of stresses in in.composition holds $\tau^{0adv}$ + // after the first advection nonlinear iteration, + // which is when the viscosity matters for the Stokes system. + const std::vector &stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress); + stress_0_advected = (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_field_indices[0]], + &in.composition[i][stress_field_indices[0]]+SymmetricTensor<2, dim>::n_independent_components)); + + // The old stresses are only changed in the operator splitting step and have been advected into + // the current timestep. They are stored in the second set of n_independent_components. + stress_old = (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]], + &in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]]+SymmetricTensor<2, dim>::n_independent_components)); + + + // Average the compositional contributions to elastic_shear_moduli here and use + // a volume-averaged shear modulus in the loop over the compositions below. + // Otherwise it is implied that each material is acting independently + // (different rotations, different stress changes), but this is inconsistent with storing only one stress tensor. + // The averaging used is the same as for the viscosity. + const std::vector &elastic_shear_moduli = elastic_rheology.get_elastic_shear_moduli(); + elastic_shear_modulus = MaterialUtilities::average_value(volume_fractions, elastic_shear_moduli, viscosity_averaging); } // Use a specified "reference" strain rate if the strain rate is not yet available @@ -188,7 +208,7 @@ namespace aspect : numbers::signaling_nan()); - // Step 1c: select which form of viscosity to use (diffusion, dislocation, fk, or composite), and apply + // Step 1c: select which form of viscosity to use (diffusion, dislocation, their minimum or composite, or fk), and apply // pre-exponential weakening, if required. switch (viscous_flow_law) { @@ -222,6 +242,16 @@ namespace aspect (scaled_viscosity_diffusion + scaled_viscosity_dislocation); break; } + case minimum_diffusion_dislocation: + { + const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \ + CompositionalViscosityPrefactors::ModifiedFlowLaws::diffusion); + const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \ + CompositionalViscosityPrefactors::ModifiedFlowLaws::dislocation); + + non_yielding_viscosity = std::min(scaled_viscosity_diffusion, scaled_viscosity_dislocation); + break; + } default: { AssertThrow(false, ExcNotImplemented()); @@ -260,7 +290,8 @@ namespace aspect if (this->get_parameters().enable_elasticity) { - const std::vector &elastic_shear_moduli = elastic_rheology.get_elastic_shear_moduli(); + // Step 3a: calculate viscoelastic (effective) viscosity + non_yielding_viscosity = elastic_rheology.calculate_viscoelastic_viscosity(non_yielding_viscosity, elastic_shear_modulus); if (use_reference_strainrate == false) { @@ -273,22 +304,20 @@ namespace aspect Assert(std::isfinite(in.strain_rate[i].norm()), ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was " "not filled by the caller.")); - const SymmetricTensor<2,dim> effective_strain_rate = - elastic_rheology.calculate_viscoelastic_strain_rate(in.strain_rate[i], + + const SymmetricTensor<2, dim> effective_strain_rate = elastic_rheology.calculate_viscoelastic_strain_rate(in.strain_rate[i], + stress_0_advected, stress_old, - elastic_shear_moduli[j]); + non_yielding_viscosity, + elastic_shear_modulus); effective_edot_ii = std::max(std::sqrt(std::max(-second_invariant(effective_strain_rate), 0.)), min_strain_rate); } - - // Step 3a: calculate the viscoelastic (effective) viscosity - non_yielding_viscosity = elastic_rheology.calculate_viscoelastic_viscosity(non_yielding_viscosity, - elastic_shear_moduli[j]); } // Step 3b: calculate non yielding (viscous or viscous + elastic) stress magnitude - double non_yielding_stress = 2. * non_yielding_viscosity * effective_edot_ii; + const double non_yielding_stress = 2. * non_yielding_viscosity * effective_edot_ii; // Step 4a: calculate the strain-weakened friction and cohesion const DruckerPragerParameters drucker_prager_parameters = drucker_prager_plasticity.compute_drucker_prager_parameters(j, @@ -351,7 +380,9 @@ namespace aspect // The following uses the effective_edot_ii // (which has been modified for elastic effects, above), // and calculates the effective viscosity over all active rheological elements - // assuming that the non-yielding viscosity is not strain rate dependent + // assuming that the non-yielding viscosity is not strain rate dependent. + // TODO When dislocation creep or nonlinear plastic strain weakening + // is used, we need to iterate on the yield stress. effective_viscosity = drucker_prager_plasticity.compute_viscosity(current_cohesion, current_friction, pressure_for_plasticity, @@ -386,6 +417,7 @@ namespace aspect ); output_parameters.composition_viscosities[j] = std::min(std::max(effective_viscosity, minimum_viscosity_for_composition), maximum_viscosity_for_composition); } + return output_parameters; } @@ -517,8 +549,11 @@ namespace aspect if (this->get_parameters().enable_elasticity) { - for (unsigned int i = 0; i < SymmetricTensor<2,dim>::n_independent_components ; ++i) - composition_mask.set(i,false); + // Set a mask (false) for the 2*n_independent_components fields representing the viscoelastic + // stress tensor components so that they are excluded from the volume fraction computation. + const std::vector &stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress); + for (unsigned int field_index: stress_field_indices) + composition_mask.set(field_index, false); } return composition_mask; @@ -560,10 +595,15 @@ namespace aspect "viscosity at that point. Select a weighted harmonic, arithmetic, " "geometric, or maximum composition."); prm.declare_entry ("Viscous flow law", "composite", - Patterns::Selection("diffusion|dislocation|frank kamenetskii|composite"), - "Select what type of viscosity law to use between diffusion, " - "dislocation, frank kamenetskii, and composite options. Soon there will be an option " - "to select a specific flow law for each assigned composition "); + Patterns::Selection("diffusion|dislocation|frank kamenetskii|composite|minimum diffusion dislocation"), + "Select what type of viscosity law to use between the options diffusion, " + "dislocation, frank kamenetskii, composite and the minimum of diffusion and dislocation. " + "Soon there will be an option " + "to select a specific flow law for each assigned composition, " + "When the full strain rate is used to compute each viscosity " + "instead of the properly partitioned diffusion and dislocation creep " + "strain rates, it is recommended to take the minimum of the creep " + "viscosities instead of their composite."); prm.declare_entry ("Yield mechanism", "drucker", Patterns::Selection("drucker|limiter"), "Select what type of yield mechanism to use between Drucker Prager " @@ -711,10 +751,12 @@ namespace aspect viscous_flow_law = dislocation; else if (prm.get ("Viscous flow law") == "frank kamenetskii") viscous_flow_law = frank_kamenetskii; + else if (prm.get ("Viscous flow law") == "minimum diffusion dislocation") + viscous_flow_law = minimum_diffusion_dislocation; else AssertThrow(false, ExcMessage("Not a valid viscous flow law")); - if (prm.get ("Yield mechanism") == "drucker") + if (prm.get("Yield mechanism") == "drucker") yield_mechanism = drucker_prager; else if (prm.get ("Yield mechanism") == "limiter") yield_mechanism = stress_limiter; @@ -836,14 +878,13 @@ namespace aspect // average over the volume volume fractions for (unsigned int j = 0; j < volume_fractions.size(); ++j) { - plastic_out->cohesions[i] += volume_fractions[j] * cohesions[j]; + plastic_out->cohesions[i] += volume_fractions[j] * cohesions[j]; // Also convert radians to degrees plastic_out->friction_angles[i] += constants::radians_to_degree * volume_fractions[j] * friction_angles_RAD[j]; plastic_out->yield_stresses[i] += volume_fractions[j] * drucker_prager_plasticity.compute_yield_stress(cohesions[j], friction_angles_RAD[j], pressure_for_plasticity, max_yield_stress); - } } } diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index f96fb8bc4b4..616cbbd5f2e 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -230,7 +230,7 @@ namespace aspect } } - // Now compute changes in the compositional fields (i.e. the accumulated strain). + // Now compute changes in the compositional fields (e.g., the accumulated strain). for (unsigned int c=0; cfill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities); if (this->get_parameters().enable_elasticity) @@ -250,12 +248,6 @@ namespace aspect average_elastic_shear_moduli[i] = MaterialUtilities::average_value(volume_fractions, rheology->elastic_rheology.get_elastic_shear_moduli(), rheology->viscosity_averaging); - - // Fill the material properties that are part of the elastic additional outputs - if (ElasticAdditionalOutputs *elastic_out = out.template get_additional_output>()) - { - elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; - } } } @@ -264,8 +256,18 @@ namespace aspect if (this->get_parameters().enable_elasticity) { + // Fill the elastic outputs with the body force term for the RHS, the viscoelastic strain rate + // and the viscous dissipation. rheology->elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out); + // Fill the elastic additional outputs with the shear modulus, elastic viscosity + // and deviatoric stress of the current timestep. + rheology->elastic_rheology.fill_elastic_additional_outputs(in, average_elastic_shear_moduli, out); + // Fill the reaction terms that account for the rotation of the stresses. rheology->elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out); + // Fill the reaction_rates that apply the stress update of the previous + // timestep to the advected and rotated stress computed in the previous timestep ($\tau^{0}$) + // to obtain $\tau^{t}$. + rheology->elastic_rheology.fill_reaction_rates(in, average_elastic_shear_moduli, out); } } @@ -398,7 +400,7 @@ namespace aspect rheology->create_plastic_outputs(out); if (this->get_parameters().enable_elasticity) - rheology->elastic_rheology.create_elastic_outputs(out); + rheology->elastic_rheology.create_elastic_additional_outputs(out); } } @@ -526,10 +528,13 @@ namespace aspect "compositional fields representing these components must be named " "and listed in a very specific format, which is designed to minimize " "mislabeling stress tensor components as distinct 'compositional " - "rock types' (or vice versa). For 2d models, the first three " - "compositional fields must be labeled 'stress\\_xx', 'stress\\_yy' and 'stress\\_xy'. " - "In 3d, the first six compositional fields must be labeled 'stress\\_xx', " - "'stress\\_yy', 'stress\\_zz', 'stress\\_xy', 'stress\\_xz', 'stress\\_yz'. " + "rock types' (or vice versa). For 2d models, three plus three consecutive " + "compositional fields must be labeled 'stress\\_xx', 'stress\\_yy', 'stress\\_xy', " + "'stress\\_xx\\_old', 'stress\\_yy\\_old', and 'stress\\_xy\\_old'. " + "In 3d, six plus six compositional fields must be labeled 'stress\\_xx', " + "'stress\\_yy', 'stress\\_zz', 'stress\\_xy', 'stress\\_xz', 'stress\\_yz', " + "'stress\\_xx\\_old', 'stress\\_yy\\_old', 'stress\\_zz\\_old', 'stress\\_xy\\_old', " + "'stress\\_xz\\_old', 'stress\\_yz\\_old'. " "\n\n " "Combining this viscoelasticity implementation with non-linear viscous flow " "and plasticity produces a constitutive relationship commonly referred to " @@ -545,10 +550,10 @@ namespace aspect "\n\n " "The overview below directly follows Moresi et al. (2003) eqns. 23-38. " "However, an important distinction between this material model and " - "the studies above is the use of compositional fields, rather than " + "the studies above is the option to use compositional fields, rather than " "particles, to track individual components of the viscoelastic stress " - "tensor. The material model will be updated when an option to track " - "and calculate viscoelastic stresses with particles is implemented. " + "tensor. Calculating viscoelastic stresses with particles is also implemented, " + "and can be switched on by using particles with the particle property 'elastic stress'. " "\n\n " "Moresi et al. (2003) begins (eqn. 23) by writing the deviatoric " "rate of deformation ($\\hat{D}$) as the sum of elastic " @@ -598,8 +603,9 @@ namespace aspect "viscosity is reduced relative to the initial viscosity. " "\n\n " "Elastic effects are introduced into the governing Stokes equations through " - "an elastic force term (eqn. 30) using stresses from the previous time step: " - "$F^{e,t} = -\\frac{\\eta_{eff}}{\\mu \\Delta t^{e}} \\tau^{t}$. " + "an elastic force term (eqn. 30 updated to the term in eqn. 5 in Farrington et al. 2014) " + "using stresses from the previous time step rotated and advected into the current time step: " + "$F^{e,t} = -\\frac{\\eta_{eff}}{\\mu \\Delta t^{e}} \\tau^{0adv}$. " "This force term is added onto the right-hand side force vector in the " "system of equations. " "\n\n " diff --git a/source/material_model/viscoelastic.cc b/source/material_model/viscoelastic.cc index a5981c907a7..85be37e32c7 100644 --- a/source/material_model/viscoelastic.cc +++ b/source/material_model/viscoelastic.cc @@ -36,13 +36,6 @@ namespace aspect { EquationOfStateOutputs eos_outputs (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1); - // Store which components to exclude during volume fraction computation. - ComponentMask composition_mask(this->n_compositional_fields(), true); - // assign compositional fields associated with viscoelastic stress a value of 0 - // assume these fields are listed first - for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i) - composition_mask.set(i,false); - std::vector average_elastic_shear_moduli (in.n_evaluation_points()); std::vector elastic_shear_moduli(elastic_rheology.get_elastic_shear_moduli()); @@ -73,26 +66,40 @@ namespace aspect for (unsigned int c=0; c *elastic_out = out.template get_additional_output>()) + // If we have multiple compositions, we need to first compute their respective their viscoelastic viscosities, + // based on their respective viscous viscosities and the averaged shear modulus, before averaging them + // into the final effective viscosity. + std::vector viscoelastic_viscosities(volume_fractions.size()); + for (unsigned int j=0; j < volume_fractions.size(); ++j) { - elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; + // The viscoelastic viscosity is scaled with the timestep ratio $\frac{\Delta t_c}{\Delta t_{el}}$ in the + // calculate_viscoelastic_viscosity function. + viscoelastic_viscosities[j] = elastic_rheology.calculate_viscoelastic_viscosity(viscosities[j], + average_elastic_shear_moduli[i]); } + + // Average viscoelastic (e.g., effective) viscosity (equation 28 in Moresi et al., 2003, J. Comp. Phys.). + out.viscosities[i] = MaterialUtilities::average_value(volume_fractions, viscoelastic_viscosities, viscosity_averaging); } + // Fill the body force term, viscoelastic strain rate and viscous dissipation. elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out); + // Fill the elastic additional outputs with the shear modulus, elastic viscosity + // and deviatoric stress of the current timestep. + elastic_rheology.fill_elastic_additional_outputs(in, average_elastic_shear_moduli, out); + // Fill the reaction terms to apply the rotation of the stresses into the current timestep. elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out); - + // Fill the reaction_rates that during operator splitting apply the stress update of the previous + // timestep to the advected and rotated stress computed in the previous timestep ($\tau^{0adv}$) + // to obtain $\tau^{t}$. + elastic_rheology.fill_reaction_rates(in, average_elastic_shear_moduli, out); } + + template bool Viscoelastic:: @@ -101,6 +108,8 @@ namespace aspect return equation_of_state.is_compressible(); } + + template void Viscoelastic::declare_parameters (ParameterHandler &prm) @@ -190,7 +199,7 @@ namespace aspect void Viscoelastic::create_additional_named_outputs (MaterialModel::MaterialModelOutputs &out) const { - elastic_rheology.create_elastic_outputs(out); + elastic_rheology.create_elastic_additional_outputs(out); } } } @@ -214,10 +223,13 @@ namespace aspect "compositional fields representing these components must be named " "and listed in a very specific format, which is designed to minimize " "mislabeling stress tensor components as distinct 'compositional " - "rock types' (or vice versa). For 2d models, the first three " - "compositional fields must be labeled 'stress\\_xx', 'stress\\_yy' and 'stress\\_xy'. " - "In 3d, the first six compositional fields must be labeled 'stress\\_xx', " - "'stress\\_yy', 'stress\\_zz', 'stress\\_xy', 'stress\\_xz', 'stress\\_yz'. " + "rock types' (or vice versa). For 2d models, the first six " + "compositional fields must be labeled 'stress\\_xx', 'stress\\_yy' and 'stress\\_xy', " + "'stress\\_xx\\_old', 'stress\\_yy\\_old' and 'stress\\_xy\\_old', " + "In 3d, the first twelve compositional fields must be labeled 'stress\\_xx', " + "'stress\\_yy', 'stress\\_zz', 'stress\\_xy', 'stress\\_xz', 'stress\\_yz', " + "'stress\\_xx\\_old', 'stress\\_yy\\_old', 'stress\\_zz\\_old', " + "'stress\\_xy\\_old', 'stress\\_xz\\_old', 'stress\\_yz\\_old'. " "\n\n " "Expanding the model to include non-linear viscous flow (e.g., " "diffusion/dislocation creep) and plasticity would produce a " @@ -263,8 +275,8 @@ namespace aspect "W^{t}\\tau^{t} + \\tau^{t}W^{t}$. " "In this material model, the size of the time step above ($\\Delta t^{e}$) " "can be specified as the numerical time step size or an independent fixed time " - "step. If the latter case is selected, the user has an option to apply a " - "stress averaging scheme to account for the differences between the numerical " + "step. If the latter case is selected, a linear interpolation will be applied" + "to account for the differences between the numerical " "and fixed elastic time step (eqn. 32). If one selects to use a fixed elastic time " "step throughout the model run, this can still be achieved by using CFL and " "maximum time step values that restrict the numerical time step to a specific time." @@ -282,8 +294,9 @@ namespace aspect "viscosity is reduced relative to the initial viscosity. " "\n\n " "Elastic effects are introduced into the governing Stokes equations through " - "an elastic force term (eqn. 30) using stresses from the previous time step: " - "$F^{e,t} = -\\frac{\\eta_\\text{eff}}{\\mu \\Delta t^{e}} \\tau^{t}$. " + "an elastic force term (eqn. 30 updated to the term in eqn. 5 in Farrington et al. 2014) " + "using stresses from the previous time step rotated and advected into the current time step: " + "$F^{e,t} = -\\frac{\\eta_\\text{eff}}{\\mu \\Delta t^{e}} \\tau^{0adv}$. " "This force term is added onto the right-hand side force vector in the " "system of equations. " "\n\n " diff --git a/source/particle/property/elastic_stress.cc b/source/particle/property/elastic_stress.cc index 203835a277f..fb8fe173a71 100644 --- a/source/particle/property/elastic_stress.cc +++ b/source/particle/property/elastic_stress.cc @@ -22,6 +22,7 @@ #include #include #include +#include namespace aspect { @@ -42,10 +43,6 @@ namespace aspect void ElasticStress::initialize () { - material_inputs = MaterialModel::MaterialModelInputs(1, this->n_compositional_fields()); - - material_outputs = MaterialModel::MaterialModelOutputs(1, this->n_compositional_fields()); - AssertThrow((Plugins::plugin_type_matches>(this->get_material_model()) || Plugins::plugin_type_matches>(this->get_material_model())), @@ -54,6 +51,187 @@ namespace aspect AssertThrow(this->get_parameters().enable_elasticity == true, ExcMessage ("This particle property should only be used if 'Enable elasticity' is set to true")); + const auto &manager = this->get_particle_manager(this->get_particle_manager_index()).get_property_manager(); + AssertThrow(!manager.plugin_name_exists("composition"), + ExcMessage("The 'elastic stress' plugin cannot be used in combination with the 'composition' plugin.")); + + material_inputs = MaterialModel::MaterialModelInputs(1, this->n_compositional_fields()); + + material_outputs = MaterialModel::MaterialModelOutputs(1, this->n_compositional_fields()); + + // The reaction rates are stored in additional outputs + this->get_material_model().create_additional_named_outputs(material_outputs); + + // Get the indices of those compositions that correspond to stress tensor elements. + stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress); + + // Get the indices of all compositions that do not correspond to stress tensor elements. + std::vector all_field_indices(this->n_compositional_fields()); + std::iota (std::begin(all_field_indices), std::end(all_field_indices), 0); + set_difference(all_field_indices.begin(), all_field_indices.end(), + stress_field_indices.begin(), stress_field_indices.end(), + std::inserter(non_stress_field_indices, non_stress_field_indices.begin())); + + // Connect to the signal after particles are restored at the beginning of + // a nonlinear iteration of iterative advection schemes. + this->get_signals().post_restore_particles.connect( + [&](typename Particle::Manager &particle_manager) + { + this->update_particles(particle_manager); + } + ); + } + + + + template + void + ElasticStress::update_particles(typename Particle::Manager &particle_manager) const + { + // There is no update of the stress to apply during the first (0th) timestep + if (this->simulator_is_past_initialization() == false || this->get_timestep_number() == 0) + return; + + // Determine the data position of the first stress tensor component + const unsigned int data_position = particle_manager.get_property_manager().get_data_info().get_position_by_field_name("ve_stress_xx"); + + // Get handler + Particle::ParticleHandler &particle_handler = particle_manager.get_particle_handler(); + + // Structure for the reaction rates + MaterialModel::ReactionRateOutputs *reaction_rate_outputs + = material_outputs.template get_additional_output>(); + + const std::vector update_flags = particle_manager.get_property_manager().get_update_flags(); + + // combine all update flags to a single flag, which is the required information + // for the mapping inside the solution evaluator + UpdateFlags mapping_flags = update_flags[0]; + for (unsigned int i=1; i> evaluators = construct_solution_evaluator(*this, + mapping_flags); + + // FEPointEvaluation uses different evaluation flags than the common UpdateFlags. + // Translate between the two. + std::vector evaluation_flags (update_flags.size(), EvaluationFlags::nothing); + + for (unsigned int i=0; i> positions; + + // Loop over all active and owned cells + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + if (cell->is_locally_owned()) + { + // Find which particles are in the current cell + typename ParticleHandler::particle_iterator_range + particles_in_cell = particle_handler.particles_in_cell(cell); + + // Only update particles, if there are any in this cell + if (particles_in_cell.begin() != particles_in_cell.end()) + { + const unsigned int n_particles_in_cell = particle_handler.n_particles_in_cell(cell); + + // Get the locations of the particles within the reference cell + positions.resize(n_particles_in_cell); + //for (auto particle = particles_in_cell.begin(); particle!=particles_in_cell.end(); ++particle) + // positions.push_back(particle->get_reference_location()); + unsigned int p = 0; + for (const auto &particle : particles_in_cell) + { + positions[p] = particle.get_reference_location(); + ++p; + } + + // Collect the values of the old solution restricted to the current cell's DOFs + small_vector old_solution_values(this->get_fe().dofs_per_cell); + cell->get_dof_values(this->get_old_solution(), + old_solution_values.begin(), + old_solution_values.end()); + + EvaluationFlags::EvaluationFlags evaluation_flags_union = EvaluationFlags::nothing; + for (unsigned int i=0; ireinit(cell, {positions.data(), positions.size()}); + + evaluators->evaluate({old_solution_values.data(),old_solution_values.size()}, + evaluation_flags); + } + + // To store the old solutions + std::vector> old_solution(n_particles_in_cell,small_vector(evaluators->n_components(), numbers::signaling_nan())); + + // To store the old gradients + std::vector,50>> old_gradients(n_particles_in_cell,small_vector,50>(evaluators->n_components(), numbers::signaling_nan>())); + + // Loop over all particles in the cell + auto particle = particles_in_cell.begin(); + for (unsigned int i = 0; particle!=particles_in_cell.end(); ++particle,++i) + { + // Evaluate the old solution, but only if it is requested in the update_flags + if (evaluation_flags_union & EvaluationFlags::values) + evaluators->get_solution(i, {&old_solution[i][0],old_solution[i].size()}, evaluation_flags); + + // Evaluate the old gradients, but only if they are requested in the update_flags + if (evaluation_flags_union & EvaluationFlags::gradients) + evaluators->get_gradients(i, {&old_gradients[i][0],old_gradients[i].size()}, evaluation_flags); + + // Fill material model input + // Get the real location of the particle + material_inputs.position[0] = particle->get_location(); + + material_inputs.current_cell = typename DoFHandler::active_cell_iterator(*particle->get_surrounding_cell(), + &(this->get_dof_handler())); + + material_inputs.temperature[0] = old_solution[i][this->introspection().component_indices.temperature]; + + material_inputs.pressure[0] = old_solution[i][this->introspection().component_indices.pressure]; + + for (unsigned int d = 0; d < dim; ++d) + material_inputs.velocity[0][d] = old_solution[i][this->introspection().component_indices.velocities[d]]; + + // Fill the non-stress composition inputs with the old_solution. + for (const unsigned int &n : non_stress_field_indices) + material_inputs.composition[0][n] = old_solution[i][this->introspection().component_indices.compositional_fields[n]]; + // Fill the ve_stress_* composition inputs with the values on the particles. + // After the particles have been restored, their properties have the values of the previous timestep. + for (unsigned int n = 0; n < 2*SymmetricTensor<2,dim>::n_independent_components; ++n) + material_inputs.composition[0][stress_field_indices[n]] = particle->get_properties()[data_position + n]; + + Tensor<2,dim> grad_u; + for (unsigned int d=0; dget_material_model().evaluate (material_inputs,material_outputs); + + // Add the reaction_rates * timestep = update to the corresponding stress + // tensor components + for (unsigned int p = 0; p < 2*SymmetricTensor<2,dim>::n_independent_components ; ++p) + particle->get_properties()[data_position + p] += reaction_rate_outputs->reaction_rates[0][stress_field_indices[p]] * this->get_timestep(); + } + } + } } @@ -71,6 +249,7 @@ namespace aspect if (dim == 2) { data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy"))); + } else if (dim == 3) { @@ -83,6 +262,24 @@ namespace aspect data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yz"))); } + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xx_old"))); + + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yy_old"))); + + if (dim == 2) + { + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old"))); + } + else if (dim == 3) + { + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_zz_old"))); + + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old"))); + + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xz_old"))); + + data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yz_old"))); + } } @@ -107,8 +304,12 @@ namespace aspect for (unsigned int d = 0; d < dim; ++d) material_inputs.velocity[0][d] = inputs.solution[p][this->introspection().component_indices.velocities[d]]; - for (unsigned int n = 0; n < this->n_compositional_fields(); ++n) + // Fill the non-stress composition inputs with the solution. + for (const unsigned int &n : non_stress_field_indices) material_inputs.composition[0][n] = inputs.solution[p][this->introspection().component_indices.compositional_fields[n]]; + // For the stress composition we use the ve_stress_* stored on the particles. + for (unsigned int n = 0; n < 2*SymmetricTensor<2,dim>::n_independent_components; ++n) + material_inputs.composition[0][stress_field_indices[n]] = particle.get_properties()[this->data_position + n]; Tensor<2,dim> grad_u; for (unsigned int d=0; dget_material_model().evaluate (material_inputs,material_outputs); + // Apply the stress rotation to the ve_stress_* fields, not the ve_stress_*_old fields. for (unsigned int i = 0; i < SymmetricTensor<2,dim>::n_independent_components ; ++i) - particle.get_properties()[this->data_position + i] += material_outputs.reaction_terms[0][i]; + particle.get_properties()[this->data_position + i] += material_outputs.reaction_terms[0][stress_field_indices[i]]; ++p; } @@ -153,31 +355,34 @@ namespace aspect { std::vector> property_information; - //Check which fields are used in model and make an output for each. - if (this->introspection().compositional_name_exists("ve_stress_xx")) - property_information.emplace_back("ve_stress_xx",1); - - if (this->introspection().compositional_name_exists("ve_stress_yy")) - property_information.emplace_back("ve_stress_yy",1); + property_information.emplace_back("ve_stress_xx",1); + property_information.emplace_back("ve_stress_yy",1); if (dim == 2) { - if (this->introspection().compositional_name_exists("ve_stress_xy")) - property_information.emplace_back("ve_stress_xy",1); + property_information.emplace_back("ve_stress_xy",1); } else if (dim == 3) { - if (this->introspection().compositional_name_exists("ve_stress_zz")) - property_information.emplace_back("ve_stress_zz",1); - - if (this->introspection().compositional_name_exists("ve_stress_xy")) - property_information.emplace_back("ve_stress_xy",1); + property_information.emplace_back("ve_stress_zz",1); + property_information.emplace_back("ve_stress_xy",1); + property_information.emplace_back("ve_stress_xz",1); + property_information.emplace_back("ve_stress_yz",1); + } - if (this->introspection().compositional_name_exists("ve_stress_xz")) - property_information.emplace_back("ve_stress_xz",1); + property_information.emplace_back("ve_stress_xx_old",1); + property_information.emplace_back("ve_stress_yy_old",1); - if (this->introspection().compositional_name_exists("ve_stress_yz")) - property_information.emplace_back("ve_stress_yz",1); + if (dim == 2) + { + property_information.emplace_back("ve_stress_xy_old",1); + } + else if (dim == 3) + { + property_information.emplace_back("ve_stress_zz_old",1); + property_information.emplace_back("ve_stress_xy_old",1); + property_information.emplace_back("ve_stress_xz_old",1); + property_information.emplace_back("ve_stress_yz_old",1); } return property_information; @@ -197,9 +402,13 @@ namespace aspect "elastic stress", "A plugin in which the particle property tensor is " "defined as the total elastic stress a particle has " - "accumulated. See the viscoelastic material model " + "accumulated. This plugin modifies the properties " + "with the name ve_stress_*. It first applies the stress " + "change resulting from system evolution during the previous " + "computational timestep, and then the rotation of those " + "stresses into the current timestep. " + "See the viscoelastic or visco_plastic material model " "documentation for more detailed information.") - } } } diff --git a/source/postprocess/visualization/maximum_horizontal_compressive_stress.cc b/source/postprocess/visualization/maximum_horizontal_compressive_stress.cc index 2ca4d413ded..3e2a2d50a9c 100644 --- a/source/postprocess/visualization/maximum_horizontal_compressive_stress.cc +++ b/source/postprocess/visualization/maximum_horizontal_compressive_stress.cc @@ -51,6 +51,7 @@ namespace aspect ExcInternalError()); Assert (input_data.solution_values[0].size() == this->introspection().n_components, ExcInternalError()); Assert (input_data.solution_gradients[0].size() == this->introspection().n_components, ExcInternalError()); + AssertThrow (this->get_parameters().enable_elasticity == false, ExcMessage("The maximum horizontal compressive stress plugin currently does not work with elasticity.")); MaterialModel::MaterialModelInputs in(input_data, this->introspection()); @@ -359,8 +360,10 @@ namespace aspect " \\label{fig:max-horizontal-compressive-stress}" "\\end{figure}" "\n\n" - "Physical units: \\si{\\pascal}." - ) + "Note that this plugin does not take into account elastic stresses " + "and therefore cannot be used when elasticity is switched on. " + "\n\n" + "Physical units: \\si{\\pascal}.") } } } diff --git a/source/postprocess/visualization/principal_stress.cc b/source/postprocess/visualization/principal_stress.cc index d5378f43ee6..2156ea8d762 100644 --- a/source/postprocess/visualization/principal_stress.cc +++ b/source/postprocess/visualization/principal_stress.cc @@ -21,6 +21,8 @@ #include #include +#include +#include namespace aspect { @@ -107,44 +109,39 @@ namespace aspect MaterialModel::MaterialModelOutputs out(n_quadrature_points, this->n_compositional_fields()); - // We do not need to compute anything but the viscosity - in.requested_properties = MaterialModel::MaterialProperties::viscosity; + // We do not need to compute anything but the viscosity and ElasticAdditionalOutputs + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; - // Compute the viscosity... + this->get_material_model().create_additional_named_outputs(out); + + // Compute the viscosity and additional outputs this->get_material_model().evaluate(in, out); for (unsigned int q=0; q stress; + const double eta = out.viscosities[q]; + + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? in.strain_rate[q] - 1. / 3. * trace(in.strain_rate[q]) * unit_symmetric_tensor() + : in.strain_rate[q]); + // Add elastic stresses if existent if (this->get_parameters().enable_elasticity == false) { - const SymmetricTensor<2,dim> deviatoric_strain_rate - = (this->get_material_model().is_compressible() - ? - in.strain_rate[q] - 1./3. * trace(in.strain_rate[q]) * unit_symmetric_tensor() - : - in.strain_rate[q]); - // Compressive stress is positive in geoscience applications - stress = -2. * out.viscosities[q] * deviatoric_strain_rate; + stress = -2. * eta * deviatoric_strain_rate; } else { - // Compressive stress is positive in geoscience applications - stress[0][0] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")]; - stress[1][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")]; - stress[0][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")]; - - if (dim == 3) - { - stress[2][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")]; - stress[0][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")]; - stress[1][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")]; - } - } + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'principal stress' postprocessor, but they have not been created.")); + + stress = -(elastic_additional_out->deviatoric_stress[q]); + } if (use_deviatoric_stress == false) stress += in.pressure[q] * unit_symmetric_tensor(); diff --git a/source/postprocess/visualization/shear_stress.cc b/source/postprocess/visualization/shear_stress.cc index 34d44e35532..19fe61a1ac5 100644 --- a/source/postprocess/visualization/shear_stress.cc +++ b/source/postprocess/visualization/shear_stress.cc @@ -20,8 +20,7 @@ #include - - +#include namespace aspect { @@ -58,13 +57,15 @@ namespace aspect MaterialModel::MaterialModelOutputs out(n_quadrature_points, this->n_compositional_fields()); - // We do not need to compute anything but the viscosity - in.requested_properties = MaterialModel::MaterialProperties::viscosity; + // We do not need to compute anything but the viscosity and ElasticAdditionalOutputs + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; - // Compute the viscosity... + this->get_material_model().create_additional_named_outputs(out); + + // Compute the viscosity and additional outputs this->get_material_model().evaluate(in, out); - // ...and use it to compute the stresses + // ...and use them to compute the stresses for (unsigned int q=0; q shear_stress; - // If elasticity is enabled, the deviatoric stress is stored + const double eta = out.viscosities[q]; + + const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() + : strain_rate); + + // If elasticity is enabled, the visco-elastic stress is stored // in compositional fields, otherwise the deviatoric stress - // can be obtained from the viscosity and strain rate. + // can be obtained from the viscosity and strain rate only. if (this->get_parameters().enable_elasticity == true) { - shear_stress[0][0] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")]; - shear_stress[1][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")]; - shear_stress[0][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")]; - - if (dim == 3) - { - shear_stress[2][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")]; - shear_stress[0][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")]; - shear_stress[1][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")]; - } + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); + + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'shear stress' postprocessor, but they have not been created.")); + + shear_stress = -(elastic_additional_out->deviatoric_stress[q]); } else { - const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q]; - const SymmetricTensor<2,dim> deviatoric_strain_rate - = (this->get_material_model().is_compressible() - ? - strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() - : - strain_rate); - - const double eta = out.viscosities[q]; shear_stress = -2. * eta * deviatoric_strain_rate; } diff --git a/source/postprocess/visualization/stress.cc b/source/postprocess/visualization/stress.cc index 4b85f1d3a49..bf1d6b6134e 100644 --- a/source/postprocess/visualization/stress.cc +++ b/source/postprocess/visualization/stress.cc @@ -20,8 +20,7 @@ #include - - +#include namespace aspect { @@ -58,13 +57,15 @@ namespace aspect MaterialModel::MaterialModelOutputs out(n_quadrature_points, this->n_compositional_fields()); - // We do not need to compute anything but the viscosity - in.requested_properties = MaterialModel::MaterialProperties::viscosity; + // We do not need to compute anything but the viscosity and the additional outputs + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; - // Compute the viscosity... + this->get_material_model().create_additional_named_outputs(out); + + // Compute the viscosity and additional outputs this->get_material_model().evaluate(in, out); - // ...and use it to compute the stresses + // ...and use them to compute the stresses for (unsigned int q=0; q stress = in.pressure[q] * unit_symmetric_tensor(); - // If elasticity is enabled, the deviatoric stress is stored - // in compositional fields, otherwise the deviatoric stress - // can be obtained from the viscosity and strain rate. + const double eta = out.viscosities[q]; + + const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() + : strain_rate); + + // If elasticity is enabled, the visco-elastic stress is stored + // in compositional fields and we can retrieve the deviatoric stress + // from the material model, otherwise the deviatoric stress + // can be computed from the viscosity and strain rate. if (this->get_parameters().enable_elasticity) { - stress[0][0] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")]; - stress[1][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")]; - stress[0][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")]; - - if (dim == 3) - { - stress[2][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")]; - stress[0][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")]; - stress[1][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")]; - } + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); + + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'principal stress' postprocessor, but they have not been created.")); + + stress -= elastic_additional_out->deviatoric_stress[q]; } else { - const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q]; - const SymmetricTensor<2,dim> deviatoric_strain_rate - = (this->get_material_model().is_compressible() - ? - strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() - : - strain_rate); - - const double eta = out.viscosities[q]; stress -= 2. * eta * deviatoric_strain_rate; } diff --git a/source/postprocess/visualization/stress_residual.cc b/source/postprocess/visualization/stress_residual.cc new file mode 100644 index 00000000000..65ebca8ee49 --- /dev/null +++ b/source/postprocess/visualization/stress_residual.cc @@ -0,0 +1,129 @@ +/* + Copyright (C) 2011 - 2021 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#include +#include + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + template + StressResidual:: + StressResidual () + : + DataPostprocessorScalar ("stress_residual", + update_values | update_gradients | update_quadrature_points) + {} + + + + template + void + StressResidual:: + evaluate_vector_field(const DataPostprocessorInputs::Vector &input_data, + std::vector> &computed_quantities) const + { + AssertThrow(Plugins::plugin_type_matches>(this->get_material_model()), + ExcMessage("This postprocessor only works with the viscoplastic material model. ")); + + const unsigned int n_quadrature_points = input_data.solution_values.size(); + Assert(computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert(computed_quantities[0].size() == 1, ExcInternalError()); + Assert(input_data.solution_values[0].size() == this->introspection().n_components, ExcInternalError()); + Assert(input_data.solution_gradients[0].size() == this->introspection().n_components, ExcInternalError()); + + MaterialModel::MaterialModelInputs in(input_data, + this->introspection()); + MaterialModel::MaterialModelOutputs out(n_quadrature_points, + this->n_compositional_fields()); + + this->get_material_model().create_additional_named_outputs(out); + + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; + in.current_cell = input_data.template get_cell(); + + // Get the viscosity + this->get_material_model().evaluate(in, out); + + for (unsigned int q = 0; q < n_quadrature_points; ++q) + { + const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() + : strain_rate); + + const double eta = out.viscosities[q]; + + // Compressive stress is positive in geoscience applications. + SymmetricTensor<2, dim> deviatoric_stress = -2. * eta * deviatoric_strain_rate; + + if (this->get_parameters().enable_elasticity == true) + { + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); + + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'principal stress' postprocessor, but they have not been created.")); + + deviatoric_stress = -(elastic_additional_out->deviatoric_stress[q]); + } + + // Compute the second moment invariant of the deviatoric stress + const double stress_invariant = std::sqrt(std::fabs(second_invariant(deviatoric_stress))); + + // Get the current yield_stress + const MaterialModel::PlasticAdditionalOutputs *plastic_output = + out.template get_additional_output>(); + const double yield_stress = plastic_output->yield_stresses[q]; + + // Compute the difference between the second stress invariant and the yield stress + computed_quantities[q](0) = stress_invariant - yield_stress; + } + + // average the values if requested + const auto &viz = this->get_postprocess_manager().template get_matching_postprocessor>(); + if (!viz.output_pointwise_stress_and_strain()) + average_quantities(computed_quantities); + } + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(StressResidual, + "stress residual", + "A visualization output object that generates output " + "for the difference between the second moment invariant " + "of the deviatoric stress tensor and the yield stress. " + "Note that this plugin currently only works when the " + "'visco plastic' material model is used. ") + } + } +} diff --git a/source/postprocess/visualization/stress_second_invariant.cc b/source/postprocess/visualization/stress_second_invariant.cc index 3195e256780..0d0ec199f4c 100644 --- a/source/postprocess/visualization/stress_second_invariant.cc +++ b/source/postprocess/visualization/stress_second_invariant.cc @@ -22,7 +22,7 @@ #include - +#include namespace aspect { @@ -58,47 +58,37 @@ namespace aspect MaterialModel::MaterialModelInputs in(input_data, this->introspection()); - in.requested_properties = MaterialModel::MaterialProperties::viscosity; + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; MaterialModel::MaterialModelOutputs out(n_quadrature_points, this->n_compositional_fields()); + this->get_material_model().create_additional_named_outputs(out); + this->get_material_model().evaluate(in, out); for (unsigned int q = 0; q < n_quadrature_points; ++q) { + const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() + : strain_rate); + + const double eta = out.viscosities[q]; + // Compressive stress is positive in geoscience applications. - SymmetricTensor<2, dim> stress = in.pressure[q] * unit_symmetric_tensor(); + SymmetricTensor<2, dim> deviatoric_stress = -2. * eta * deviatoric_strain_rate; - // Add elastic stresses if existent. if (this->get_parameters().enable_elasticity == true) { - stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")]; - stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")]; - stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")]; - - if (dim == 3) - { - stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")]; - stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")]; - stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")]; - } - } - else - { - const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; - const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() - ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() - : strain_rate); + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); - const double eta = out.viscosities[q]; + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'principal stress' postprocessor, but they have not been created.")); - stress += -2. * eta * deviatoric_strain_rate; + deviatoric_stress = -(elastic_additional_out->deviatoric_stress[q]); } - // Compute the deviatoric stress tensor after elastic stresses were added. - const SymmetricTensor<2, dim> deviatoric_stress = deviator(stress); - // Compute the second moment invariant of the deviatoric stress // in the same way as the second moment invariant of the deviatoric // strain rate is computed in the viscoplastic material model. diff --git a/source/postprocess/visualization/surface_stress.cc b/source/postprocess/visualization/surface_stress.cc index 56ebf1c6f3b..83dd7a34ea2 100644 --- a/source/postprocess/visualization/surface_stress.cc +++ b/source/postprocess/visualization/surface_stress.cc @@ -20,8 +20,7 @@ #include - - +#include namespace aspect { @@ -58,8 +57,10 @@ namespace aspect MaterialModel::MaterialModelOutputs out(n_quadrature_points, this->n_compositional_fields()); + this->get_material_model().create_additional_named_outputs(out); + // We do not need to compute anything but the viscosity - in.requested_properties = MaterialModel::MaterialProperties::viscosity; + in.requested_properties = MaterialModel::MaterialProperties::viscosity | MaterialModel::MaterialProperties::additional_outputs; // Compute the viscosity... this->get_material_model().evaluate(in, out); @@ -67,33 +68,36 @@ namespace aspect // ...and use it to compute the stresses for (unsigned int q=0; q strain_rate = in.strain_rate[q]; - const SymmetricTensor<2,dim> deviatoric_strain_rate - = (this->get_material_model().is_compressible() - ? - strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() - : - strain_rate); + // Compressive stress is negative by the sign convention + // used by the engineering community, and as input and used + // internally by ASPECT. + // Here, we change the sign of the stress to match the + // sign convention used by the geoscience community. + SymmetricTensor<2,dim> stress = in.pressure[q] * unit_symmetric_tensor(); const double eta = out.viscosities[q]; - // Compressive stress is positive in geoscience applications - SymmetricTensor<2,dim> stress = -2.*eta*deviatoric_strain_rate + - in.pressure[q] * unit_symmetric_tensor(); + const SymmetricTensor<2, dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible() + ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor() + : strain_rate); - // Add elastic stresses if existent - if (this->get_parameters().enable_elasticity == true) + // If elasticity is enabled, the visco-elastic stress is stored + // in compositional fields and we can retrieve the deviatoric stress + // from the material model, otherwise the deviatoric stress + // can be computed from the viscosity and strain rate. + if (this->get_parameters().enable_elasticity) + { + // Get the total deviatoric stress from the material model. + const MaterialModel::ElasticAdditionalOutputs *elastic_additional_out = out.template get_additional_output>(); + + Assert(elastic_additional_out != nullptr, ExcMessage("Elastic Additional Outputs are needed for the 'principal stress' postprocessor, but they have not been created.")); + + stress -= elastic_additional_out->deviatoric_stress[q]; + } + else { - stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")]; - stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")]; - stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")]; - - if (dim == 3) - { - stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")]; - stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")]; - stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")]; - } + stress -= 2. * eta * deviatoric_strain_rate; } for (unsigned int d=0; dcreate_additional_named_outputs(scratch.material_model_outputs); MaterialModel::ReactionRateOutputs *reaction_rate_outputs diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc index 5fed6560d1c..6b68fd89141 100644 --- a/source/simulator/solver_schemes.cc +++ b/source/simulator/solver_schemes.cc @@ -337,6 +337,56 @@ namespace aspect if (fields_advected_by_particles.size() > 0) interpolate_particle_properties(fields_advected_by_particles); + // If elasticity is switched on, the stress tensor can have components + // with very large values and components that are and stay practically zero. + // Therefore, we include all fields that represent stress tensor components + // of the same stress tensor in the computation of the initial residual + // for the fields that belong to that tensor. In other words, we compute an + // averaged initial residual using those fields that belong to the ve_stress tensor. + // The ve_stress_old values are not independent components of the solution vector, + // so we do not need to compute a residual for them and set their residual to zero. + // TODO Is this residual calculation representative of a second order tensor? + if (parameters.enable_elasticity == true) + { + double stress_initial_residual = 0.0; + std::vector stress_indices; + std::vector old_stress_indices; + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xx")); + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yy")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xx_old")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yy_old")); + if (dim == 2) + { + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy_old")); + } + else if (dim == 3) + { + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_zz")); + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy")); + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xz")); + stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yz")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_zz_old")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy_old")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xz_old")); + old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yz_old")); + } + + + if (residual) + { + const double n_stress_fields = stress_indices.size(); + for (auto &c : stress_indices) + stress_initial_residual += system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm() / n_stress_fields; + + // Overwrite the initial residual + for (auto &c : stress_indices) + (*residual)[c] = stress_initial_residual; + for (auto &c : old_stress_indices) + (*residual)[c] = 0.; + } + } + // for consistency we update the current linearization point only after we have solved // all fields, so that we use the same point in time for every field when solving @@ -346,7 +396,9 @@ namespace aspect = solution.block(introspection.block_indices.compositional_fields[c]); if ((initial_residual.size() > 0) && (initial_residual[c] > 0)) - current_residual[c] /= initial_residual[c]; + { + current_residual[c] /= initial_residual[c]; + } else current_residual[c] = 0.0; } @@ -948,6 +1000,15 @@ namespace aspect for (auto &particle_manager : particle_managers) particle_manager.restore_particles(); + // Apply a particle update if required by the particle properties. + // Apply the update even if nonlinear_iteration == 0, + // because the signal could be used, for example, to apply operator + // splitting on the particles, and this would need to be + // applied at the beginning of the timestep and after + // every restore_particles(). + for (auto &particle_manager : particle_managers) + signals.post_restore_particles(particle_manager); + const double relative_temperature_residual = assemble_and_solve_temperature(initial_temperature_residual, nonlinear_iteration == 0 ? &initial_temperature_residual : nullptr); @@ -1109,6 +1170,15 @@ namespace aspect for (auto &particle_manager : particle_managers) particle_manager.restore_particles(); + // Apply a particle update if required by the particle properties. + // Apply the update even if nonlinear_iteration == 0, + // because the signal could be used, for example, to apply operator + // splitting on the particles, and this would need to be + // applied at the beginning of the timestep and after + // every restore_particles(). + for (auto &particle_manager : particle_managers) + signals.post_restore_particles(particle_manager); + const double relative_temperature_residual = assemble_and_solve_temperature(initial_temperature_residual, nonlinear_iteration == 0 ? &initial_temperature_residual : nullptr); diff --git a/tests/boundary_traction_function_cartesian_free_surface.prm b/tests/boundary_traction_function_cartesian_free_surface.prm index bb1f55c6f78..64e66d7e7cf 100644 --- a/tests/boundary_traction_function_cartesian_free_surface.prm +++ b/tests/boundary_traction_function_cartesian_free_surface.prm @@ -5,12 +5,24 @@ set Dimension = 2 set Use years in output instead of seconds = true set Start time = 0 set End time = 5 -set Nonlinear solver scheme = single Advection, iterated Stokes +set Nonlinear solver scheme = iterated Advection and Stokes set CFL number = 0.5 set Maximum time step = 2.5 set Timing output frequency = 1 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 + subsection Geometry model set Model name = chunk @@ -41,6 +53,7 @@ subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 1 + set Use discontinuous composition discretization = true end subsection Formulation @@ -76,8 +89,9 @@ subsection Boundary traction model end 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 subsection Initial composition model @@ -86,7 +100,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 diff --git a/tests/boundary_traction_function_spherical_free_surface.prm b/tests/boundary_traction_function_spherical_free_surface.prm index c45db3508a8..1aa1559c1e4 100644 --- a/tests/boundary_traction_function_spherical_free_surface.prm +++ b/tests/boundary_traction_function_spherical_free_surface.prm @@ -5,12 +5,24 @@ set Dimension = 2 set Use years in output instead of seconds = true set Start time = 0 set End time = 5 -set Nonlinear solver scheme = single Advection, iterated Stokes +set Nonlinear solver scheme = iterated Advection and Stokes set CFL number = 0.5 set Maximum time step = 2.5 set Timing output frequency = 1 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 + subsection Geometry model set Model name = chunk @@ -41,6 +53,7 @@ subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 1 + set Use discontinuous composition discretization = true end subsection Formulation @@ -77,8 +90,9 @@ subsection Boundary traction model end 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 subsection Initial composition model @@ -87,7 +101,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 diff --git a/tests/damped_viscoelastic_bending_beam.prm b/tests/damped_viscoelastic_bending_beam.prm index 51ff4a8c6cc..fb690b9e5e3 100644 --- a/tests/damped_viscoelastic_bending_beam.prm +++ b/tests/damped_viscoelastic_bending_beam.prm @@ -26,9 +26,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, 2800, 2800, 2800, 2800, 2800, 2800, 3300 + set Viscosities = 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e18, 1.e24 + set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e10 set Fixed elastic time step = 1e3 set Elastic damper viscosity = 1e14 set Use fixed elastic time step = false diff --git a/tests/free_surface_timestep_repeat.prm b/tests/free_surface_timestep_repeat.prm index c5a69ff7dfe..fe68c47ca69 100644 --- a/tests/free_surface_timestep_repeat.prm +++ b/tests/free_surface_timestep_repeat.prm @@ -7,13 +7,24 @@ set Dimension = 2 set Start time = 0 set End time = 2 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.001 set Maximum time step = 1 set Output directory = output_free_surface_VE_cylinder_2D_loading-repeat set Timing output frequency = 10 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 + subsection Time stepping set List of model names = convection time step, repeat on cutback set Minimum time step size = 0.5 @@ -50,6 +61,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 @@ -91,8 +104,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 @@ -102,7 +116,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 @@ -136,7 +150,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 = 3 + set Fixed elastic time step = 1 set Viscosity averaging scheme = harmonic end end diff --git a/tests/newton_postprocess_nonlinear.prm b/tests/newton_postprocess_nonlinear.prm index b563800640c..30cf8d601c6 100644 --- a/tests/newton_postprocess_nonlinear.prm +++ b/tests/newton_postprocess_nonlinear.prm @@ -6,7 +6,7 @@ set Dimension = 2 set Start time = 0 set End time = 0 set Use years in output instead of seconds = true -set Nonlinear solver scheme = single Advection, iterated Newton Stokes +set Nonlinear solver scheme = iterated Advection and Newton Stokes set Nonlinear solver tolerance = 1e-4 set Max nonlinear iterations = 5 set CFL number = 0.5 @@ -14,6 +14,17 @@ set Maximum time step = 5 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 parameters subsection Solver parameters subsection Stokes solver parameters @@ -60,6 +71,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 @@ -81,9 +94,9 @@ 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 Types of fields = stress, stress, stress, chemical composition, chemical composition, chemical composition + 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 @@ -93,7 +106,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); @@ -140,16 +153,18 @@ subsection Material model set Reference strain rate = 1.e-13 set Maximum viscosity = 1.e23 set Minimum viscosity = 1.e17 - 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. set Elastic shear moduli = 1e50 set Use fixed elastic time step = false - set Fixed elastic time step = 5e3 + set Fixed elastic time step = 5 set Viscosity averaging scheme = harmonic - 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 diff --git a/tests/sea_level_postprocessor.prm b/tests/sea_level_postprocessor.prm index b07ad6a1198..94e38f3b3e7 100644 --- a/tests/sea_level_postprocessor.prm +++ b/tests/sea_level_postprocessor.prm @@ -46,6 +46,31 @@ subsection Nullspace removal set Remove nullspace = angular momentum end +subsection Compositional fields + # Although we set these stress tensor component fields, + # the material model chosen below does not implement elasticity. + 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 + +subsection Initial composition model + set Model name = function + + subsection Function + set Variable names = x,y,z + set Function constants = + set Function expression = 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; + end +end + +subsection Boundary composition model + set Allow fixed composition on outflow boundaries = true + set List of model names = initial composition +end + subsection Boundary temperature model set Fixed temperature boundary indicators = inner, outer set List of model names = initial temperature diff --git a/tests/skip_refinement.prm b/tests/skip_refinement.prm index 193151e2736..0c4f3d971e6 100644 --- a/tests/skip_refinement.prm +++ b/tests/skip_refinement.prm @@ -8,11 +8,23 @@ set Dimension = 2 set Start time = 0 set End time = 100 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 = 50 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 + # Solver parameters subsection Solver parameters subsection Stokes solver parameters @@ -54,6 +66,7 @@ subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 1 + set Use discontinuous composition discretization = true end # Formulation classification @@ -88,8 +101,9 @@ end # Number and name of compositional fields subsection Compositional fields - set Number of fields = 7 - set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, LM, TZ, UM, Lith + set Number of fields = 10 + 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, LM, TZ, UM, Lith + set Types of fields = stress, stress, stress, stress, stress, stress, chemical composition, chemical composition, chemical composition, chemical composition end # Spatial domain of different compositional fields @@ -99,7 +113,7 @@ subsection Initial composition model subsection Function set Variable names = x,y set Function constants = - set Function expression = 0; 0; 0; if(y<=2221e3,1,0); if(y<=2471e3 && y>2221e3,1,0); if(y<=2821e3 && y>2471e3,1,0); if(y>2821e3,1,0); + set Function expression = 0; 0; 0; 0; 0; 0; if(y<=2221e3,1,0); if(y<=2471e3 && y>2221e3,1,0); if(y<=2821e3 && y>2471e3,1,0); if(y>2821e3,1,0); end end @@ -130,11 +144,11 @@ subsection Material model set Material averaging = harmonic average only viscosity subsection Viscoelastic - set Densities = 4978, 4978, 4978, 4978, 4978, 3871, 3438, 3037 - set Viscosities = 2.e21, 2.e21, 2.e21, 2.e21, 2.e21, 1.e21, 1.e21, 1.e40 - set Elastic shear moduli = 2.28340e11, 2.28340e11, 2.28340e11, 2.28340e11, 2.28340e11, 1.05490e11, 0.70363e11, 0.50605e11 + set Densities = 4978, 4978, 3871, 3438, 3037 + set Viscosities = 2.e21, 2.e21, 1.e21, 1.e21, 1.e40 + set Elastic shear moduli = 2.28340e11, 2.28340e11, 1.05490e11, 0.70363e11, 0.50605e11 set Use fixed elastic time step = false - set Fixed elastic time step = 2.5 + set Fixed elastic time step = 50 set Viscosity averaging scheme = harmonic end end diff --git a/tests/vep_chemical_compositions.prm b/tests/vep_chemical_compositions.prm index 813cd99a568..3297f00eedc 100644 --- a/tests/vep_chemical_compositions.prm +++ b/tests/vep_chemical_compositions.prm @@ -6,14 +6,27 @@ include $ASPECT_SOURCE_DIR/tests/visco_plastic_yield_plastic_viscous_strain_weak # Global parameters set Dimension = 2 set End time = 1000 +set Nonlinear solver scheme = iterated Advection and Stokes set Max nonlinear iterations = 2 set Output directory = vep_chemical_compositions +# One operator splitting step to update the stresses +set Use operator splitting = true + +# Solver settings +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 + # 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, plastic_strain, viscous_strain, test - set Types of fields = stress, stress, stress, strain, strain, chemical composition + 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, plastic_strain, viscous_strain, test + set Types of fields = stress, stress, stress, stress, stress, stress, strain, strain, chemical composition end # We prescribe some initial strain at the center of the domain @@ -23,7 +36,7 @@ subsection Initial composition model subsection Function set Variable names = x,y - set Function expression = 0; 0; 0; \ + set Function expression = 0; 0; 0; 0; 0; 0; \ if(x>=45e3&x<=55e3&y>=45.3e3&y<=55.e3,0.2,0); \ if(x>=35e3&x<=65e3&y>=35.e3&y<=65.e3,0.8,0); 1 end @@ -39,6 +52,11 @@ subsection Formulation set Enable elasticity = true end +# Element types +subsection Discretization + set Use discontinuous composition discretization = true +end + # Material model (values for background material) subsection Material model set Model name = visco plastic @@ -50,19 +68,19 @@ subsection Material model subsection Visco Plastic set Reference strain rate = 1.e-16 set Viscous flow law = dislocation - set Prefactors for dislocation creep = 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23 - set Stress exponents for dislocation creep = 1, 0, 0, 0, 0, 0, 1 + set Prefactors for dislocation creep = 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23, 5e-23 + set Stress exponents for dislocation creep = 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 set Activation energies for dislocation creep = 0, 0 set Activation volumes for dislocation creep = 0 set Yield mechanism = drucker - set Angles of internal friction = 0,10,10,10,5,5,0 - set Cohesions = 1.e6,0,0,0,0,0,1e6 + set Angles of internal friction = 0,10,10,10,10,10,10,5,5,0 + set Cohesions = 1.e6,0,0,0,0,0,0,0,0,1e6 set Strain weakening mechanism = plastic weakening with plastic strain and viscous weakening with viscous strain set Start prefactor strain weakening intervals = 0,0 - set End prefactor strain weakening intervals = 1,1,1,1,1,1,1 + set End prefactor strain weakening intervals = 1,1,1,1,1,1,1,1,1,1 set Prefactor strain weakening factors = 0.5,0.5 set Start plasticity strain weakening intervals = 0.,0 - set End plasticity strain weakening intervals = 1,1,1,1,1,1,1 + set End plasticity strain weakening intervals = 1,1,1,1,1,1,1,1,1,1 set Cohesion strain weakening factors = 0.5,0.5 set Friction strain weakening factors = 0.5 set Elastic shear moduli = 1e10,1e10 diff --git a/tests/visco_elastic_top_beam.prm b/tests/visco_elastic_top_beam.prm index e27fc9f5fe3..3df7856d184 100644 --- a/tests/visco_elastic_top_beam.prm +++ b/tests/visco_elastic_top_beam.prm @@ -6,11 +6,22 @@ set Dimension = 2 set End time = 2e3 set Use years in output instead of seconds = true set Output directory = output_ve +set Nonlinear solver scheme = iterated Advection and Stokes +set Max nonlinear iterations = 5 + set CFL number = 0.5 set Maximum time step = 1e3 +# One operator splitting step to update the stresses +set Use operator splitting = true + # Solver settings 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 subsection Stokes solver parameters set Number of cheap Stokes solver steps = 2000 end @@ -44,8 +55,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 @@ -62,8 +73,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 @@ -73,7 +85,7 @@ subsection Initial composition model subsection Function set Variable names = x,y set Function constants = - set Function expression = 0; 0; 0; if (x>=3.5e3 && x<=4.0e3 && y>=3.0e3, 1, 0) + set Function expression = 0; 0; 0; 0; 0; 0; if (x>=3.5e3 && x<=4.0e3 && y>=3.0e3, 1, 0) end end @@ -110,13 +122,13 @@ subsection Material model set Model name = visco plastic subsection Visco Plastic - set Densities = 2800, 2800, 2800, 2800, 6000 + set Densities = 2800, 6000 set Viscous flow law = dislocation - set Prefactors for dislocation creep = 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-25 + set Prefactors for dislocation creep = 5.e-19, 5.e-25 set Stress exponents for dislocation creep = 1.0 set Activation energies for dislocation creep = 0. set Activation volumes for dislocation creep = 0. - set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e10 + set Elastic shear moduli = 1.e11, 1.e10 set Fixed elastic time step = 2e3 set Use fixed elastic time step = true set Viscosity averaging scheme = maximum composition diff --git a/tests/visco_plastic_vep_bending_beam.prm b/tests/visco_plastic_vep_bending_beam.prm index 28d8e3927b2..e12fc11787e 100644 --- a/tests/visco_plastic_vep_bending_beam.prm +++ b/tests/visco_plastic_vep_bending_beam.prm @@ -22,13 +22,13 @@ subsection Material model set Model name = visco plastic subsection Visco Plastic - set Densities = 2800, 2800, 2800, 2800, 3300 + set Densities = 2800, 2800, 2800, 2800, 2800, 2800, 2800, 3300 set Viscous flow law = dislocation - set Prefactors for dislocation creep = 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-25 + set Prefactors for dislocation creep = 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-25 set Stress exponents for dislocation creep = 1.0 set Activation energies for dislocation creep = 0. set Activation volumes for dislocation creep = 0. - set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e10 + set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 1.e11, 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/tests/visco_plastic_vep_brick_extension.prm b/tests/visco_plastic_vep_brick_extension.prm index dfe551c4d14..4dc87522084 100644 --- a/tests/visco_plastic_vep_brick_extension.prm +++ b/tests/visco_plastic_vep_brick_extension.prm @@ -37,7 +37,7 @@ subsection Initial composition model subsection Function set Variable names = x,y set Function constants = - set Function expression = 0; 0; 0; 0; \ + set Function expression = 0; 0; 0; 0; 0; 0; 0; \ if (y<=2.0e3 && x>=18.0e3 && x<=22.0e3, 1, 0); end end diff --git a/tests/visco_plastic_vep_plate_flexure.prm b/tests/visco_plastic_vep_plate_flexure.prm index eff50b0744c..33f9c8cc04b 100644 --- a/tests/visco_plastic_vep_plate_flexure.prm +++ b/tests/visco_plastic_vep_plate_flexure.prm @@ -27,16 +27,16 @@ subsection Material model subsection Visco Plastic set Viscous flow law = dislocation - set Prefactors for dislocation creep = 5.e-18, 5.e-18, 5.e-18, 5.e-18, 5.e-36, 5.e-36 + set Prefactors for dislocation creep = 5.e-18, 5.e-18, 5.e-18, 5.e-18, 5.e-18, 5.e-18, 5.e-18, 5.e-36, 5.e-36 set Stress exponents for dislocation creep = 1.0 set Activation energies for dislocation creep = 0. set Activation volumes for dislocation creep = 0 - set Elastic shear moduli = 1.e50, 1.e50, 1.e50, 1.e50, 3.e10, 3.e10 + set Elastic shear moduli = 1.e50, 1.e50, 1.e50, 1.e50, 1.e50, 1.e50, 1.e50, 3.e10, 3.e10 set Fixed elastic time step = 5 set Use fixed elastic time step = false set Viscosity averaging scheme = maximum composition set Cohesions = 1.e50 - set Densities = 2700, 2700, 2700, 2700, 2700, 1890 + set Densities = 2700, 2700, 2700, 2700, 2700, 2700, 2700, 2700, 1890 end end diff --git a/tests/visco_plastic_vep_stress_build-up_yield.prm b/tests/visco_plastic_vep_stress_build-up_yield.prm index 69185522072..f4274a38134 100644 --- a/tests/visco_plastic_vep_stress_build-up_yield.prm +++ b/tests/visco_plastic_vep_stress_build-up_yield.prm @@ -19,7 +19,7 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm -set Nonlinear solver scheme = single Advection, iterated Stokes +set Nonlinear solver scheme = iterated Advection and Stokes set Nonlinear solver tolerance = 1e-5 set Max nonlinear iterations = 2 diff --git a/tests/visco_plastic_vep_stress_build-up_yield_particles.prm b/tests/visco_plastic_vep_stress_build-up_yield_particles.prm index 635acf1d99c..6bec43a547c 100644 --- a/tests/visco_plastic_vep_stress_build-up_yield_particles.prm +++ b/tests/visco_plastic_vep_stress_build-up_yield_particles.prm @@ -19,7 +19,6 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm -set Nonlinear solver scheme = single Advection, iterated Stokes set Nonlinear solver tolerance = 1e-5 set Max nonlinear iterations = 2 set Output directory = visco_plastic_vep_stress_build-up_yield_particles @@ -27,6 +26,10 @@ set Output directory = visco_plastic_vep_stress_build-up_y # Global parameters set End time = 2e3 +# For particles, the operator splitting step +# is handled by the particle property elastic stress. +set Use operator splitting = false + # Model geometry (100x100 km, 5 km spacing) subsection Geometry model set Model name = box @@ -41,8 +44,9 @@ end # Track compositional fields with active particles subsection Compositional fields - set Compositional field methods = particles, particles, particles - set Mapped particle properties = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy + set Compositional field methods = particles, particles, particles, particles, particles, particles + 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 end # Material model diff --git a/tests/viscoelastic_bending_beam.prm b/tests/viscoelastic_bending_beam.prm index ce81b7b0196..c9dae49d935 100644 --- a/tests/viscoelastic_bending_beam.prm +++ b/tests/viscoelastic_bending_beam.prm @@ -31,21 +31,5 @@ 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 Types of fields = stress, stress, stress, chemical composition -end - -# Material model -subsection Material model - set Model name = viscoelastic - - subsection Viscoelastic - set Densities = 2800, 3300 - set Viscosities = 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 Viscosity averaging scheme = maximum composition - end + set Types of fields = stress, stress, stress, stress, stress, stress, chemical composition end diff --git a/tests/viscoelastic_bending_beam_particles.prm b/tests/viscoelastic_bending_beam_particles.prm index 876872cfe0d..d2c74cf65f3 100644 --- a/tests/viscoelastic_bending_beam_particles.prm +++ b/tests/viscoelastic_bending_beam_particles.prm @@ -5,6 +5,10 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm +# For particles, the stress update is done by the particle property elastic stress, +# so we do not use the operator splitting. +set Use operator splitting = false + set End time = 2e3 set Output directory = viscoelastic_bending_beam_particles @@ -20,10 +24,11 @@ end # Compositional field names and methods 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/tests/viscoelastic_numerical_elastic_time_step_particles.prm b/tests/viscoelastic_numerical_elastic_time_step_particles.prm index fd943255954..5fd28e67482 100644 --- a/tests/viscoelastic_numerical_elastic_time_step_particles.prm +++ b/tests/viscoelastic_numerical_elastic_time_step_particles.prm @@ -16,6 +16,10 @@ include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_ set End time = 2e3 set Output directory = viscoelastic_numerical_elastic_time_step_particles +# For particles, the stress update is done by the particle property elastic stress, +# so we do not need the operator splitting. +set Use operator splitting = false + # Model geometry (100x100 km, 5 km spacing) subsection Geometry model set Model name = box @@ -30,8 +34,9 @@ end # Track compositional fields with active particles subsection Compositional fields - set Compositional field methods = particles, particles, particles - set Mapped particle properties = ve_stress_xx:ve_stress_xx, ve_stress_yy:ve_stress_yy, ve_stress_xy:ve_stress_xy + set Compositional field methods = particles, particles, particles, particles, particles, particles + 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 end # Material model diff --git a/tests/viscoelastic_stress_build-up.prm b/tests/viscoelastic_stress_build-up.prm index 07b892908cd..496e8c0bc55 100644 --- a/tests/viscoelastic_stress_build-up.prm +++ b/tests/viscoelastic_stress_build-up.prm @@ -61,7 +61,7 @@ subsection Boundary velocity model end # Spatial domain of different compositional fields: -# ve_stress_xx, ve_stress_yy, ve_stress_xy +# ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old subsection Initial composition model set Model name = function @@ -69,6 +69,9 @@ subsection Initial composition model set Variable names = x,y set Function constants = eta = 1e22, mu = 1e10, eps = 1e-14, t = 250000, yinsec = 31557600 set Function expression = 2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ + -2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ + 0 ;\ + 2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ -2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ 0 end diff --git a/tests/viscoelastic_stress_build-up_gmg.prm b/tests/viscoelastic_stress_build-up_gmg.prm index 9c1f3235417..36d1bfad9ed 100644 --- a/tests/viscoelastic_stress_build-up_gmg.prm +++ b/tests/viscoelastic_stress_build-up_gmg.prm @@ -73,7 +73,7 @@ subsection Boundary velocity model end # Spatial domain of different compositional fields: -# ve_stress_xx, ve_stress_yy, ve_stress_xy +# ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old subsection Initial composition model set Model name = function @@ -81,6 +81,9 @@ subsection Initial composition model set Variable names = x,y set Function constants = eta = 1e22, mu = 1e10, eps = 1e-14, t = 250000, yinsec = 31557600 set Function expression = 2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ + -2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ + 0;\ + 2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ -2.*eta*eps*(1.-exp(-mu*t*yinsec/eta)); \ 0 end diff --git a/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm b/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm index 20f1706d585..77a8d9ec28e 100644 --- a/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm +++ b/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm @@ -4,13 +4,26 @@ include $ASPECT_SOURCE_DIR/tests/visco_plastic_yield_plastic_viscous_strain_weak # Global parameters set Dimension = 2 set End time = 1000 +set Nonlinear solver scheme = iterated Advection and Stokes set Max nonlinear iterations = 2 set Output directory = viscoelastoplastic_yield_plastic_viscous_strain_weakening +# 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 + # Number and name of compositional fields subsection Compositional fields - set Number of fields = 5 - set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, plastic_strain, viscous_strain + set Number of fields = 8 + 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, plastic_strain, viscous_strain + set Types of fields = stress, stress, stress, stress, stress, stress, strain, strain end # We prescribe some initial strain at the center of the domain @@ -20,7 +33,7 @@ subsection Initial composition model subsection Function set Variable names = x,y - set Function expression = 0; 0; 0; \ + set Function expression = 0; 0; 0; 0; 0; 0; \ if(x>=45e3&x<=55e3&y>=45.3e3&y<=55.e3,0.2,0); \ if(x>=35e3&x<=65e3&y>=35.e3&y<=65.e3,0.8,0) end @@ -36,6 +49,12 @@ subsection Formulation set Enable elasticity = true end +subsection Discretization + # DG for viscoelastic stresses + set Use discontinuous composition discretization = true +end + + # Material model (values for background material) subsection Material model set Model name = visco plastic diff --git a/tests/viscous_top_beam.prm b/tests/viscous_top_beam.prm index 63b291dda6f..0249bf6103a 100644 --- a/tests/viscous_top_beam.prm +++ b/tests/viscous_top_beam.prm @@ -69,6 +69,10 @@ subsection Boundary velocity model end # Number and name of compositional fields +# Note that although compositional fields are set with names corresponding to the +# stress tensor, elasticity is not actually enabled. Since no field types are +# specified, the type will be deduced from the field names as 'stress' and therefore +# these fields do not contribute to material properties. subsection Compositional fields set Number of fields = 4 set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, beam diff --git a/tests/visualization_shear_stress_viscoelastic.prm b/tests/visualization_shear_stress_viscoelastic.prm index 0e5a0b948c8..5ba0db5314c 100644 --- a/tests/visualization_shear_stress_viscoelastic.prm +++ b/tests/visualization_shear_stress_viscoelastic.prm @@ -18,12 +18,32 @@ include $ASPECT_SOURCE_DIR/tests/visualization_shear_stress.prm +set Nonlinear solver scheme = iterated Advection and Stokes +set Nonlinear solver tolerance = 1e-5 +set Max nonlinear iterations = 5 + set End time = 1 subsection Formulation set Enable elasticity = true end +# One operator splitting step to update the stresses +set Use operator splitting = true + +# Solver settings +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 Discretization + set Use discontinuous composition discretization = true +end + subsection Material model set Model name = viscoelastic @@ -37,8 +57,17 @@ subsection Material model end 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 + +subsection Initial composition model + set Model name = function + + subsection Function + set Function expression = 0; 0; 0; 0; 0; 0; + end end subsection Postprocess diff --git a/tests/visualization_stress_viscoelastic.prm b/tests/visualization_stress_viscoelastic.prm index b402535e662..0e4120b37da 100644 --- a/tests/visualization_stress_viscoelastic.prm +++ b/tests/visualization_stress_viscoelastic.prm @@ -18,12 +18,32 @@ include $ASPECT_SOURCE_DIR/tests/visualization_stress.prm +set Nonlinear solver scheme = iterated Advection and Stokes +set Nonlinear solver tolerance = 1e-5 +set Max nonlinear iterations = 5 + set End time = 1 subsection Formulation set Enable elasticity = true end +# One operator splitting step to update the stresses +set Use operator splitting = true + +# Solver settings +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 Discretization + set Use discontinuous composition discretization = true +end + subsection Material model set Model name = viscoelastic @@ -37,8 +57,17 @@ subsection Material model end 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 + +subsection Initial composition model + set Model name = function + + subsection Function + set Function expression = 0; 0; 0; 0; 0; 0; + end end subsection Postprocess