diff --git a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm new file mode 100644 index 00000000000..fd4baf6d198 --- /dev/null +++ b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm @@ -0,0 +1,184 @@ +# This parameter file reproduces a 2D benchmark for viscoelastic deformation +# published in "Numerical modelling of magma dynamics coupled to tectonic +# deformation of lithosphere and crust" by Keller, May and Kaus, 2013, +# Geophys. J. Int., v. 195, p. 1406-1422. The model setup and results are +# described in detail in section B.2.3 and figure B5. +# +# The benchmark examines bending and unbending (e.g., recovery) of a +# viscoelastic beam surrounded by a less dense and viscous (inelastic) fluid. +# Gravitational forces drive initial bending (elastic strain) of the beam for +# 50 Kyr, which then recovers its shape over ~ 500 Kyr if gravity is turned +# off. The recovery is driven by the accumulated elastic stresses and thus +# provides a basic test for the viscoelasticity implementation. +# +# Compositional fields are used to track the viscoelastic stresses and material +# representing the beam. To improve the accuracy of tracking the beam interface, +# a discontinuous discretization and bound preserving limiter is used for the +# compositional fields. Significantly, the time step is limited to 1 Kyr through +# the "Maximum time step" parameter and a relatively high (0.5) CFL value. Using +# a constant time step ensures the effective (viscoelastic) viscosity of the +# viscoelastic beam and viscous fluid remains constants throughout the model run. +# +# As currently constructed, the model will run for 50 Kyr with a gravitational +# acceleration of 10 m/s^2. To produce unbending of the beam after the model +# has finished, change the parameter "End time" to 500 Kyr and set the +# gravitational acceleration to 0 m/s^2. A restart file is written every 10 Kyr +# and the parameter "Resume computation = auto" specifies that a model should +# restart from a checkpoint file if one is present. + +# Global parameters +set Dimension = 2 +set Start time = 0 +set End time = 50e3 +set Use years in output instead of seconds = true +set Resume computation = auto +set CFL number = 0.5 +set Maximum time step = 1e3 +set Output directory = output +set Pressure normalization = surface +set Surface pressure = 0. + +# Solver settings +subsection Solver parameters + subsection Stokes solver parameters + set Number of cheap Stokes solver steps = 2000 + end +end + +# Model geometry (7.5x5 km, 0.1 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 75 + set Y repetitions = 50 + set X extent = 7.5e3 + set Y extent = 5e3 + end +end + +# Mesh refinement specifications +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 + set Time steps between mesh refinement = 0 +end + +# Element types +subsection Discretization + set Composition polynomial degree = 2 + set Stokes velocity polynomial degree = 2 + set Temperature polynomial degree = 1 + set Use locally conservative discretization = false + set Use discontinuous temperature discretization = false + set Use discontinuous composition discretization = true + subsection Stabilization parameters + set Use limiter for discontinuous composition solution = true + set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0 + set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0 + end +end + +# Formulation classification +subsection Formulation + set Enable elasticity = true +end + +# Velocity boundary conditions +subsection Boundary velocity model + set Zero velocity boundary indicators = left + set Tangential velocity boundary indicators = bottom, top, right +end + +# Number and name of compositional fields +subsection Compositional fields + set Number of fields = 4 + set Names of fields = stress_xx, stress_yy, stress_xy, beam +end + +# Spatial domain of different compositional fields +subsection Initial composition model + set Model name = function + subsection Function + set Variable names = x,y + set Function constants = + set Function expression = 0; 0; 0; if (x<=4.5e3 && y>=2.5e3 && y<=3.0e3, 1, 0) + end +end + +# Composition boundary conditions +subsection Boundary composition model + set Fixed composition boundary indicators = bottom, top, right + set List of model names = initial composition +end + +# Temperature boundary conditions +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top, left, right + set List of model names = box + subsection Box + set Bottom temperature = 293 + set Left temperature = 293 + set Right temperature = 293 + set Top temperature = 293 + end +end + +# Temperature initial conditions +subsection Initial temperature model + set Model name = function + subsection Function + set Function expression = 293 + end +end + +# Material model +subsection Material model + + set Model name = viscoelastic + + subsection Viscoelastic + set Densities = 2800, 2800, 2800, 2800, 3300 + set Viscosities = 1.e18, 1.e18, 1.e18, 1.e18, 1.e24 + set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e10 + set Fixed elastic time step = 1e3 + set Use fixed elastic time step = false + set Use stress averaging = false + set Viscosity averaging scheme = maximum composition + end + +end + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 10. + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization + subsection Visualization + set List of output variables = material properties, strain rate + + subsection Material properties + set List of material properties = density, viscosity + end + + set Time between graphical output = 0 + set Interpolate output = true + end + +end + +# Termination criteria +subsection Termination criteria + set Termination criteria = end step + set End step = 500 +end + +subsection Checkpointing + set Steps between checkpoint = 10 +end diff --git a/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm new file mode 100644 index 00000000000..740578c967c --- /dev/null +++ b/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm @@ -0,0 +1,192 @@ +# This parameter file reproduces a 2D analytical benchmark for the build-up +# of stress in an initially unstressed viscoelastic medium subject to a +# constant strain-rate. This benchmark is described in "Robust +# characteristics method for modelling multiphase visco-elasto-plastic +# thermo-mechanical problems" by Gerya and Yuen, 2007, Phys. Earth. Planet. +# Inter., v. 163, p. 83-105. Full details of the benchmark are located in +# section 3.1 and figure 3 of this manuscript. +# +# The analytical solution for viscoelastic stress build-up in +# an incompressible medium with a constant strain-rate is: +# simga_ij = 2 * edot_ii * eta * (1 - e^(-mu*t/eta)), +# where sigma_ij is the elastic deviatoric stress tensor, edot_ij is +# the deviatoric strain-rate, eta is the background viscosity, mu is the +# shear modulus and t is time. +# +# Following the conditions described in section 3.1 and figure 3 from +# Gerya and Yuen (2007), a 100x100 km body is subject to a constant +# strain-rate of 1.e-14 s^-1 in both the horizontal and vertical directions. +# Constant deformation is driven by inflow and outflow, respectively, +# on the right and bottom walls. The top and left walls are free-slip. +# The material has a viscosity of 1e22 Pa s and a shear modulus of 1e10 Pa. +# With these values, the analytical solution predicts a horizontal +# or vertical viscoelastic stress magnitude of ~ 200 MPa after 250 Kyr +# of deformation. Significantly, the effective "viscoelastic" viscosity +# is enforced to be constant through time by using a constant time step. +# This is achieved by setting a maximum time step (1000 years) much lower +# than the time step size given by the CFL number of 0.5. +# +# This result can be observed by viewing the compositional field values +# representing the horizontal (stress_xx) or vertical (stress_yy) components +# of the viscoselastic stress tensor. Significantly, the stress near the +# model boundaries is incorrect due to the compositional field boundary +# conditions, which are based on the initial compositional values (e.g., zero). +# This leads to oscillations in the stress field near the boundaries, which +# decay towards a constant value in the model interior. This phenomen highlights +# one advantage to tracking viscoelastic stress with particles, which do not require +# defining compositional boundary conditions. + +# Global parameters +set Dimension = 2 +set Start time = 0 +set End time = 250e3 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = single Advection, single Stokes +set Nonlinear solver tolerance = 1e-5 +set Max nonlinear iterations = 1 +set CFL number = 0.5 +set Maximum time step = 1000 +set Output directory = output +set Timing output frequency = 1 +set Pressure normalization = surface +set Surface pressure = 0. + +# Solver settings +subsection Solver parameters + subsection Stokes solver parameters + set Use direct solver for Stokes system = false + set Linear solver tolerance = 1e-7 + set Number of cheap Stokes solver steps = 2000 + end +end + +# Model geometry (100x100 km, 2 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 50 + set Y repetitions = 50 + set X extent = 100e3 + set Y extent = 100e3 + end +end + +# Mesh refinement specifications +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 + set Time steps between mesh refinement = 0 +end + +# Element types +subsection Discretization + set Composition polynomial degree = 2 + set Stokes velocity polynomial degree = 2 + set Temperature polynomial degree = 1 +end + +# Formulation classification +subsection Formulation + set Enable elasticity = true +end + +# Velocity boundary conditions +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, left + set Prescribed velocity boundary indicators = bottom y:function, right x:function + subsection Function + set Variable names = x,y + set Function constants = cm=0.01, year=1, vel=3.154 + set Function expression = if (x>50e3 , vel*cm/year, 0.); if (y<50e3 , vel*cm/year, 0.); + end +end + +# Number and name of compositional fields +subsection Compositional fields + set Number of fields = 3 + set Names of fields = stress_xx, stress_yy, stress_xy +end + +# Spatial domain of different compositional fields +subsection Initial composition model + set Model name = function + subsection Function + set Variable names = x,y + set Function constants = + set Function expression = 0; 0; 0; + end +end + +# Composition boundary conditions +subsection Boundary composition model + set Fixed composition boundary indicators = bottom, top, left, right + set List of model names = initial composition +end + +# Temperature boundary conditions +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top, left, right + set List of model names = box + subsection Box + set Bottom temperature = 293 + set Left temperature = 293 + set Right temperature = 293 + set Top temperature = 293 + end +end + +# Temperature initial conditions +subsection Initial temperature model + set Model name = function + subsection Function + set Function expression = 293 + end +end + +# Material model +subsection Material model + + set Model name = viscoelastic + + subsection Viscoelastic + set Densities = 2800 + set Viscosities = 1.e22 + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = false + set Fixed elastic time step = 1e3 + set Use stress averaging = false + set Viscosity averaging scheme = harmonic + end + +end + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 0. + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = basic statistics, composition statistics, temperature statistics, velocity statistics, visualization + + subsection Visualization + set List of output variables = material properties, strain rate + + subsection Material properties + set List of material properties = density, viscosity + end + + set Time between graphical output = 10e3 + set Interpolate output = true + end + +end + +# Termination criteria +subsection Termination criteria + set Termination criteria = end time +end diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h index 32645019c63..7f82847c6a3 100644 --- a/include/aspect/material_model/interface.h +++ b/include/aspect/material_model/interface.h @@ -32,6 +32,7 @@ #include #include #include +#include namespace aspect { @@ -778,6 +779,38 @@ namespace aspect std::vector rhs_melt_pc; }; + /** + * A class for an elastic force term to be added to the RHS of the + * Stokes system, which can be attached to the + * MaterialModel::MaterialModelOutputs structure and filled in the + * MaterialModel::Interface::evaluate() function. + */ + template + class ElasticOutputs: public AdditionalMaterialOutputs + { + public: + ElasticOutputs(const unsigned int n_points) + : elastic_force(n_points, numbers::signaling_nan >() ) + {} + + virtual ~ElasticOutputs() + {} + + virtual void average (const MaterialAveraging::AveragingOperation operation, + const FullMatrix &/*projection_matrix*/, + const FullMatrix &/*expansion_matrix*/) + { + AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented()); + return; + } + + /** + * Force tensor (elastic terms) on the right-hand side for the conservation of + * momentum equation (first part of the Stokes equation) in each + * quadrature point. + */ + std::vector > elastic_force; + }; /** diff --git a/include/aspect/material_model/viscoelastic.h b/include/aspect/material_model/viscoelastic.h new file mode 100644 index 00000000000..80f5d355df3 --- /dev/null +++ b/include/aspect/material_model/viscoelastic.h @@ -0,0 +1,334 @@ +/* + Copyright (C) 2014 - 2018 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 doc/COPYING. If not see + . +*/ + +#ifndef _aspect_material_model_viscoelastic_h +#define _aspect_material_model_viscoelastic_h + +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + 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. + */ + template + class ElasticAdditionalOutputs : public NamedAdditionalMaterialOutputs + { + public: + ElasticAdditionalOutputs(const unsigned int n_points); + + virtual std::vector get_nth_output(const unsigned int idx) const; + + /** + * Elastic shear moduli at the evaluation points passed to + * the instance of MaterialModel::Interface::evaluate() that fills + * the current object. + */ + std::vector elastic_shear_moduli; + }; + + /** + * An implementation of a simple linear viscoelastic rheology that only + * includes the deviatoric components of elasticity. Specifically, the + * viscoelastic rheology only takes into account the elastic shear + * strength (e.g., shear modulus), while the tensile and volumetric + * strength (e.g., Young's and bulk modulus) are not considered. The model + * is incompressible and allows specifying an arbitrary number of + * compositional fields, where each field represents a different rock type + * or component of the viscoelastic stress tensor. The stress tensor in 2D + * and 3D, respectively, contains 3 or 6 components. The 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. + * + * Expanding the model to include non-linear viscous flow (e.g., + * diffusion/dislocation creep) and plasticity would produce a constitutive + * relationship commonly referred to as partial elastoviscoplastic + * (e.g., pEVP) in the geodynamics community. While extensively discussed + * and applied within the geodynamics literature, notable references include: + * Moresi et al. (2003), J. Comp. Phys., v. 184, p. 476-497. + * Gerya and Yuen (2007), Phys. Earth. Planet. Inter., v. 163, p. 83-105. + * Gerya (2010), Introduction to Numerical Geodynamic Modeling. + * Kaus (2010), Tectonophysics, v. 484, p. 36-47. + * Choi et al. (2013), J. Geophys. Res., v. 118, p. 2429-2444. + * Keller et al. (2013), Geophys. J. Int., v. 195, p. 1406-1442. + * + * 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 + * tracers, to track individual components of the viscoelastic stress + * tensor. The material model will be udpated when an option to track + * and calculate viscoelastic stresses with tracers is implemented. + * + * 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: + * $\hat{D} = \hat{D_{e}} + \hat{D_{v}}$. + * These terms further decompose into + * $\hat{D_{v}} = \frac{\tau}{2\eta}$ and + * $\hat{D_{e}} = \frac{\overset{\triangledown}{\tau}}{2\mu}$, where + * $\tau$ is the viscous deviatoric stress, $\eta$ is the shear viscosity, + * $\mu$ is the shear modulus and $\overset{\triangledown}{\tau}$ is the + * Jaumann corotational stress rate. This later term (eqn. 24) contains the + * time derivative of the deviatoric stress ($\dot{\tau}$) and terms that + * account for material spin (e.g., rotation) due to advection: + * $\overset{\triangledown}{\tau} = \dot{\tau} + {\tau}W -W\tau$. + * Above, $W$ is the material spin tensor (eqn. 25): + * $W_{ij} = \frac{1}{2} \left (\frac{\partial V_{i}}{\partial x_{j}} - + * \frac{\partial V_{j}}{\partial x_{i}} \right )$. + * + * The Jaumann stress-rate can also be approximated using terms from the time + * at the previous time step ($t$) and current time step ($t + \Delta t_^{e}$): + * $\smash[t]{\overset{\triangledown}{\tau}}^{t + \Delta t^{e}} \approx + * \frac{\tau^{t + \Delta t^{e} - \tau^{t}}}{\Delta t^{e}} - + * 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. + * + * The formulation above allows rewriting the total rate of deformation (eqn. 29) as + * $\tau^{t + \Delta t^{e}} = \eta_{eff} \left ( + * 2\\hat{D}^{t + \\triangle t^{e}} + \\frac{\\tau^{t}}{\\mu \\Delta t^{e}} + + * \\frac{W^{t}\\tau^{t} - \\tau^{t}W^{t}}{\\mu} \\right )$. + * + * The effective viscosity (eqn. 28) is a function of the viscosity ($\eta$), + * elastic time step size ($\Delta t^{e}$) and shear relaxation time + * ($ \alpha = \frac{\eta}{\mu} $): + * $\eta_{eff} = \eta \frac{\Delta t^{e}}{\Delta t^{e} + \alpha}$ + * 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}$. + * This force term is added onto the right-hand side force vector in the + * system of equations. + * + * The value of each compositional field representing distinct + * rock types at a point is interpreted to be a volume fraction of that + * rock type. If the sum of the compositional field volume fractions is + * less than one, then the remainder of the volume is assumed to be + * 'background material'. + * + * Several model parameters (densities, elastic shear moduli, + * thermal expansivities, thermal conductivies, specific heats) can + * be defined per-compositional field. For each material parameter the + * 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 + * 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, + * except density, which varies linearly with temperature according to the + * thermal expansivity. + * + * When more than one compositional field is present at a point, they are + * averaged arithmetically. An exception is viscosity, which may be averaged + * arithmetically, harmonically, geometrically, or by selecting the + * viscosity of the composition with the greatest volume fraction. + * + * @ingroup MaterialModels + */ + template + class Viscoelastic : public MaterialModel::Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Function to compute the material properties in @p out given the + * inputs in @p in. If MaterialModelInputs.strain_rate has the length + * 0, then the viscosity does not need to be computed. + */ + virtual void evaluate(const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const; + + /** + * @name Qualitative properties one can ask a material model + * @{ + */ + + /** + * This model is not compressible, so this returns false. + */ + virtual bool is_compressible () const; + /** + * @} + */ + + /** + * @name Reference quantities + * @{ + */ + virtual double reference_viscosity () const; + /** + * @} + */ + + /** + * @name Functions used in dealing with run-time parameters + * @{ + */ + + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + /** + * @} + */ + + virtual + void + create_additional_named_outputs (MaterialModel::MaterialModelOutputs &out) const; + + + private: + /** + * The first 3 (2D) or 6 (3D) compositional fields are assumed + * to be components of the viscoelastic stress tensor and + * assigned volume fractions of zero. + */ + const std::vector compute_volume_fractions( + const std::vector &compositional_fields) const; + + /** + * Reference temperature for thermal expansion. All components use + * the same reference_T. + */ + double reference_T; + + /** + * Enumeration for selecting which averaging scheme to use. Select + * between harmonic, arithmetic, geometric, and maximum_composition. + * The max composition scheme simply uses the parameter of whichever + * field has the highest volume fraction. + */ + enum AveragingScheme + { + harmonic, + arithmetic, + geometric, + maximum_composition + }; + + + AveragingScheme viscosity_averaging; + + double average_value (const std::vector &composition, + const std::vector ¶meter_values, + const enum AveragingScheme &average_type) const; + + + /** + * Used for calculating average elastic shear modulus and viscosity + */ + double calculate_average_vector (const std::vector &composition, + const std::vector ¶meter_values, + const enum AveragingScheme &average_type) const; + + + double calculate_average_viscoelastic_viscosity (const double average_viscosity, + const double average_elastic_shear_modulus, + const double dte) const; + + + /** + * Vector for field densities, read from parameter file. + */ + std::vector densities; + + /** + * Vector for field viscosities, read from parameter file. + */ + std::vector viscosities; + + /** + * Vector for field thermal expnsivities, read from parameter file. + */ + std::vector thermal_expansivities; + + /** + * Vector for field thermal conductivities, read from parameter file. + */ + std::vector thermal_conductivities; + + /** + * Vector for field specific heats, read from parameter file. + */ + std::vector specific_heats; + + /** + * Vector for field elastic shear moduli, read from parameter file. + */ + std::vector elastic_shear_moduli; + + /** + * Bool indicating whether to use a fixed material time scale in the + * viscoelastic rheology for all time steps (if true) or to use the + * actual (variable) advection time step of the model (if false). Read + * from parameter file. + */ + bool use_fixed_elastic_time_step; + + /** + * Bool indicating whether to use a stress averaging scheme to account + * for differences between the numerical and fixed elastic time step + * (if true). When set to false, the viscoelastic stresses are not + * modified to account for differences between the viscoelastic time + * step and the numerical time step. Read from parameter file. + */ + bool use_stress_averaging; + + /** + * Double for fixed elastic time step value, read from parameter file + */ + double fixed_elastic_time_step; + + }; + + } +} + +#endif diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h index 6c7d6ad4bfa..60b0134cfdd 100644 --- a/include/aspect/parameters.h +++ b/include/aspect/parameters.h @@ -372,6 +372,12 @@ namespace aspect */ typename Formulation::TemperatureEquation::Kind formulation_temperature_equation; + /** + * This variable determines whether additional terms related to elastic forces + * are added to the Stokes equation. + */ + bool enable_elasticity; + /** * @} */ diff --git a/source/material_model/viscoelastic.cc b/source/material_model/viscoelastic.cc new file mode 100644 index 00000000000..067714dbbdd --- /dev/null +++ b/source/material_model/viscoelastic.cc @@ -0,0 +1,711 @@ +/* + Copyright (C) 2011 - 2018 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 doc/COPYING. If not see + . +*/ + + +#include +#include +#include +#include +#include +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + + namespace + { + std::vector make_elastic_additional_outputs_names() + { + std::vector names; + names.push_back("elastic_shear_modulus"); + return names; + } + } + + template + ElasticAdditionalOutputs::ElasticAdditionalOutputs (const unsigned int n_points) + : + NamedAdditionalMaterialOutputs(make_elastic_additional_outputs_names()), + elastic_shear_moduli(n_points, numbers::signaling_nan()) + {} + + template + std::vector + ElasticAdditionalOutputs::get_nth_output(const unsigned int idx) const + { + AssertIndexRange (idx, 1); + switch (idx) + { + case 0: + return elastic_shear_moduli; + + default: + AssertThrow(false, ExcInternalError()); + } + // we will never get here, so just return something + return elastic_shear_moduli; + } + + template + const std::vector + Viscoelastic:: + compute_volume_fractions(const std::vector &compositional_fields) const + { + std::vector volume_fractions(compositional_fields.size()+1); + + //clip the compositional fields so they are between zero and one + std::vector x_comp = compositional_fields; + for (unsigned int i=0; i < x_comp.size(); ++i) + x_comp[i] = std::min(std::max(x_comp[i], 0.0), 1.0); + + // assign compositional fields associated with viscoelastic stress a value of 0 + for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i) + x_comp[i] = 0.0; + + //sum the compositional fields for normalization purposes + double sum_composition = 0.0; + for (unsigned int i=0; i < x_comp.size(); ++i) + sum_composition += x_comp[i]; + + if (sum_composition >= 1.0) + { + volume_fractions[0] = 0.0; //background mantle + for (unsigned int i=1; i <= x_comp.size(); ++i) + volume_fractions[i] = x_comp[i-1]/sum_composition; + } + else + { + volume_fractions[0] = 1.0 - sum_composition; //background mantle + for (unsigned int i=1; i <= x_comp.size(); ++i) + volume_fractions[i] = x_comp[i-1]; + } + return volume_fractions; + } + + template + double + Viscoelastic:: + average_value (const std::vector &volume_fractions, + const std::vector ¶meter_values, + const enum AveragingScheme &average_type) const + { + double averaged_parameter = 0.0; + + switch (average_type) + { + case arithmetic: + { + for (unsigned int i=0; i< volume_fractions.size(); ++i) + averaged_parameter += volume_fractions[i]*parameter_values[i]; + break; + } + case harmonic: + { + for (unsigned int i=0; i< volume_fractions.size(); ++i) + { + AssertThrow(parameter_values[i] != 0, + ExcMessage ("All values must be greater than 0 during harmonic averaging")); + averaged_parameter += volume_fractions[i]/(parameter_values[i]); + } + averaged_parameter = 1.0/averaged_parameter; + break; + } + case geometric: + { + for (unsigned int i=0; i < volume_fractions.size(); ++i) + averaged_parameter += volume_fractions[i]*std::log(parameter_values[i]); + averaged_parameter = std::exp(averaged_parameter); + break; + } + case maximum_composition: + { + const unsigned int i = (unsigned int)(std::max_element( volume_fractions.begin(), + volume_fractions.end() ) + - volume_fractions.begin()); + averaged_parameter = parameter_values[i]; + break; + } + default: + { + AssertThrow( false, ExcNotImplemented() ); + break; + } + } + return averaged_parameter; + } + + template + double + Viscoelastic:: + calculate_average_vector (const std::vector &composition, + const std::vector ¶meter_values, + const enum AveragingScheme &average_type) const + { + const std::vector volume_fractions = compute_volume_fractions(composition); + const double averaged_vector = average_value(volume_fractions, parameter_values, average_type); + return averaged_vector; + } + + template + double + Viscoelastic:: + calculate_average_viscoelastic_viscosity (const double average_viscosity, + const double average_elastic_shear_modulus, + const double dte) const + { + return ( average_viscosity * dte ) / ( dte + ( average_viscosity / average_elastic_shear_modulus ) ); + } + + + template + void + Viscoelastic:: + evaluate(const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const + { + + // Create the structure for the elastic force terms that are needed to compute the + // right-hand side of the Stokes system + MaterialModel::ElasticOutputs + *force_out = out.template get_additional_output >(); + + + // The elastic time step (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 + // 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. + const double dte = ( ( this->get_timestep_number() > 0 && use_fixed_elastic_time_step == false ) + ? + this->get_timestep() + : + fixed_elastic_time_step); + + for (unsigned int i=0; i < in.temperature.size(); ++i) + { + const double temperature = in.temperature[i]; + const std::vector composition = in.composition[i]; + const std::vector volume_fractions = compute_volume_fractions(composition); + + out.specific_heat[i] = average_value(volume_fractions, specific_heats, arithmetic); + + // Arithmetic averaging of thermal conductivities + // This may not be strictly the most reasonable thing, but for most Earth materials we hope + // that they do not vary so much that it is a big problem. + out.thermal_conductivities[i] = average_value(volume_fractions, thermal_conductivities, arithmetic); + + double density = 0.0; + for (unsigned int j=0; j < volume_fractions.size(); ++j) + { + // not strictly correct if thermal expansivities are different, since we are interpreting + // these compositions as volume fractions, but the error introduced should not be too bad. + const double temperature_factor= (1.0 - thermal_expansivities[j] * (temperature - reference_T)); + density += volume_fractions[j] * densities[j] * temperature_factor; + } + out.densities[i] = density; + + out.thermal_expansion_coefficients[i] = average_value(volume_fractions, thermal_expansivities, arithmetic); + + // Compressibility at the given positions. + // The compressibility is given as + // $\frac 1\rho \frac{\partial\rho}{\partial p}$. + // (here we use an incompressible medium) + out.compressibilities[i] = 0.0; + // Pressure derivative of entropy at the given positions. + out.entropy_derivative_pressure[i] = 0.0; + // Temperature derivative of entropy at the given positions. + out.entropy_derivative_temperature[i] = 0.0; + // Change in composition due to chemical reactions at the + // given positions. The term reaction_terms[i][c] is the + // change in compositional field c at point i. + for (unsigned int c=0; c *elastic_out = out.template get_additional_output >()) + { + elastic_out->elastic_shear_moduli[i] = average_elastic_shear_modulus; + } + + // Fill elastic force outputs (assumed to be zero during intial time step) + if (force_out) + { + force_out->elastic_force[i] = 0.; + } + } + + // Viscoelasticity section + if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && in.strain_rate.size() > 0) + { + // Get old (previous time step) velocity gradients + std::vector > quadrature_positions(in.position.size()); + for (unsigned int i=0; i < in.position.size(); ++i) + quadrature_positions[i] = this->get_mapping().transform_real_to_unit_cell(in.current_cell, in.position[i]); + + // FEValues requires a quadrature and we provide the default quadrature + // as we only need to evaluate the solution and gradients. + FEValues fe_values (this->get_mapping(), + this->get_fe(), + Quadrature(quadrature_positions), + update_gradients); + + fe_values.reinit (in.current_cell); + std::vector > old_velocity_gradients (quadrature_positions.size(), Tensor<2,dim>()); + fe_values[this->introspection().extractors.velocities].get_function_gradients (this->get_old_solution(), + old_velocity_gradients); + + MaterialModel::ElasticOutputs + *force_out = out.template get_additional_output >(); + + for (unsigned int i=0; i < in.position.size(); ++i) + { + // Get old stresses from compositional fields + SymmetricTensor<2,dim> stress_old; + 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]; + + // Calculate the rotated stresses + // Rotation (vorticity) tensor (equation 25 in Moresi et al., 2003, J. Comp. Phys.) + const Tensor<2,dim> rotation = 0.5 * ( old_velocity_gradients[i] - transpose(old_velocity_gradients[i]) ); + + // Recalculate average elastic shear modulus + const std::vector composition = in.composition[i]; + const double average_elastic_shear_modulus = calculate_average_vector(composition, + elastic_shear_moduli, + viscosity_averaging); + + // Average viscoelastic viscosity + const double average_viscoelastic_viscosity = out.viscosities[i]; + + // Calculate the current (new) viscoelastic stress, which is a function of the material + // properties (viscoelastic viscosity, shear modulus), elastic time step size, strain rate, + // vorticity and prior (inherited) viscoelastic stresses (see equation 29 in Moresi et al., + // 2003, J. Comp. Phys.) + SymmetricTensor<2,dim> stress_new = ( 2. * average_viscoelastic_viscosity * deviator(in.strain_rate[i]) ) + + ( ( average_viscoelastic_viscosity / ( average_elastic_shear_modulus * dte ) ) * stress_old ) + + ( ( average_viscoelastic_viscosity / average_elastic_shear_modulus ) * + ( symmetrize(rotation * Tensor<2,dim>(stress_old) ) - symmetrize(Tensor<2,dim>(stress_old) * rotation) ) ); + + // Stress averaging scheme to account for difference betweed fixed elastic time step + // and numerical time step (see equation 32 in Moresi et al., 2003, J. Comp. Phys.) + const double dt = this->get_timestep(); + if (use_fixed_elastic_time_step == true && use_stress_averaging == true) + { + stress_new = ( ( 1. - ( dt / dte ) ) * stress_old ) + ( ( dt / dte ) * stress_new ) ; + } + + // Fill reaction terms + for (unsigned int j = 0; j < SymmetricTensor<2,dim>::n_independent_components ; ++j) + out.reaction_terms[i][j] = -stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] + + stress_new[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)]; + + // Fill elastic force outputs (See equation 30 in Moresi et al., 2003, J. Comp. Phys.) + if (force_out) + { + force_out->elastic_force[i] = -1. * ( ( average_viscoelastic_viscosity / ( average_elastic_shear_modulus * dte ) ) * stress_old ); + } + + } + } + + } + + template + double + Viscoelastic:: + reference_viscosity () const + { + return viscosities[0]; //background + } + + template + bool + Viscoelastic:: + is_compressible () const + { + return false; + } + + template + void + Viscoelastic::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Material model"); + { + prm.enter_subsection("Viscoelastic"); + { + prm.declare_entry ("Reference temperature", "293", + Patterns::Double (0), + "The reference temperature $T_0$. Units: $K$."); + prm.declare_entry ("Densities", "3300.", + Patterns::List(Patterns::Double(0)), + "List of densities for background mantle and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "If only one value is given, then all use the same value. Units: $kg / m^3$"); + prm.declare_entry ("Viscosities", "1.e21", + Patterns::List(Patterns::Double(0)), + "List of viscosities for background mantle and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "If only one value is given, then all use the same value. Units: $Pa s$"); + prm.declare_entry ("Thermal expansivities", "4.e-5", + Patterns::List(Patterns::Double(0)), + "List of thermal expansivities for background mantle and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "If only one value is given, then all use the same value. Units: $1/K$"); + prm.declare_entry ("Specific heats", "1250.", + Patterns::List(Patterns::Double(0)), + "List of specific heats $C_p$ for background mantle and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "If only one value is given, then all use the same value. Units: $J /kg /K$"); + prm.declare_entry ("Thermal conductivities", "4.7", + Patterns::List(Patterns::Double(0)), + "List of thermal conductivities for background mantle and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "If only one value is given, then all use the same value. Units: $W/m/K$ "); + prm.declare_entry ("Viscosity averaging scheme", "harmonic", + Patterns::Selection("arithmetic|harmonic|geometric|maximum composition "), + "When more than one compositional field is present at a point " + "with different viscosities, we need to come up with an average " + "viscosity at that point. Select a weighted harmonic, arithmetic, " + "geometric, or maximum composition."); + prm.declare_entry ("Elastic shear moduli", "75.0e9", + Patterns::List(Patterns::Double(0)), + "List of elastic shear moduli, $G$, " + "for background material and compositional fields, " + "for a total of N+1 values, where N is the number of compositional fields. " + "The default value of 75 GPa is representative of mantle rocks. Units: Pa."); + prm.declare_entry ("Use fixed elastic time step", "unspecified", + Patterns::Selection("true|false|unspecified"), + "Select whether the material time scale in the viscoelastic constitutive " + "relationship uses the regular numerical time step or a separate fixed " + "elastic time step throughout the model run. The fixed elastic time step " + "is always used during the initial time step. If a fixed elastic time " + "step is used throughout the model run, a stress averaging scheme can be " + "applied to account for differences with the numerical time step. An " + "alternative approach is to limit the maximum time step size so that it " + "is equal to the elastic time step. The default value of this parameter is " + "'unspecified', which throws an exception during runtime. In order for " + "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; " + "seconds otherwise."); + prm.declare_entry ("Use stress averaging","false", + Patterns::Bool (), + "Whether to apply a stress averaging scheme to account for differences " + "between the fixed elastic time step and numerical time step. "); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + template + void + Viscoelastic::parse_parameters (ParameterHandler &prm) + { + + // Get the number of fields for composition-dependent material properties + const unsigned int n_fields = this->n_compositional_fields() + 1; + + prm.enter_subsection("Material model"); + { + prm.enter_subsection("Viscoelastic"); + { + reference_T = prm.get_double ("Reference temperature"); + + if (prm.get ("Viscosity averaging scheme") == "harmonic") + viscosity_averaging = harmonic; + else if (prm.get ("Viscosity averaging scheme") == "arithmetic") + viscosity_averaging = arithmetic; + else if (prm.get ("Viscosity averaging scheme") == "geometric") + viscosity_averaging = geometric; + else if (prm.get ("Viscosity averaging scheme") == "maximum composition") + viscosity_averaging = maximum_composition; + else + AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme")); + + // Parse viscoelastic properties + densities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Densities"))), + n_fields, + "Densities"); + viscosities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Viscosities"))), + n_fields, + "Viscosities"); + thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivities"))), + n_fields, + "Thermal conductivities"); + thermal_expansivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal expansivities"))), + n_fields, + "Thermal expansivities"); + specific_heats = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Specific heats"))), + n_fields, + "Specific heats"); + elastic_shear_moduli = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Elastic shear moduli"))), + n_fields, + "Elastic shear moduli"); + + if (prm.get ("Use fixed elastic time step") == "true") + use_fixed_elastic_time_step = true; + else if (prm.get ("Use fixed elastic time step") == "false") + use_fixed_elastic_time_step = false; + else + AssertThrow(false, ExcMessage("'Use fixed elastic time step' must be set to 'true' or 'false'")); + + use_stress_averaging = prm.get_bool ("Use stress averaging"); + if (use_stress_averaging) + AssertThrow(use_fixed_elastic_time_step == true, + ExcMessage("Stress averaging can only be used if 'Use fixed elastic time step' is set to true'")); + + fixed_elastic_time_step = prm.get_double ("Fixed elastic time step"); + AssertThrow(fixed_elastic_time_step > 0, + ExcMessage("The fixed elastic time step must be greater than zero")); + + 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. + if (dim == 2) + { + AssertThrow(this->introspection().compositional_index_for_name("stress_xx") == 0, + ExcMessage("Material model Viscoelastic only works if the first " + "compositional field is called stress_xx.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_yy") == 1, + ExcMessage("Material model Viscoelastic only works if the second " + "compositional field is called stress_yy.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_xy") == 2, + ExcMessage("Material model Viscoelastic only works if the third " + "compositional field is called stress_xy.")); + } + else if (dim == 3) + { + AssertThrow(this->introspection().compositional_index_for_name("stress_xx") == 0, + ExcMessage("Material model Viscoelastic only works if the first " + "compositional field is called stress_xx.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_yy") == 1, + ExcMessage("Material model Viscoelastic only works if the second " + "compositional field is called stress_yy.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_zz") == 2, + ExcMessage("Material model Viscoelastic only works if the third " + "compositional field is called stress_zz.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_xy") == 3, + ExcMessage("Material model Viscoelastic only works if the fourth " + "compositional field is called stress_xy.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_xz") == 4, + ExcMessage("Material model Viscoelastic only works if the fifth " + "compositional field is called stress_xz.")); + AssertThrow(this->introspection().compositional_index_for_name("stress_yz") == 5, + ExcMessage("Material model Viscoelastic only works if the sixth " + "compositional field is called stress_yz.")); + } + else + 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. + AssertThrow((this->get_parameters().nonlinear_solver == + Parameters::NonlinearSolver::single_Advection_single_Stokes + || + this->get_parameters().nonlinear_solver == + Parameters::NonlinearSolver::single_Advection_iterated_Stokes), + ExcMessage("The material model will only work with the nonlinear " + "solver schemes 'single Advection, single Stokes' and " + "'single Advection, iterated Stokes'")); + + // Functionality to average the additional RHS terms over the cell is not implemented. + // This enforces that the variable 'Material averaging' is set to 'none'. + AssertThrow((this->get_parameters().material_averaging == + MaterialModel::MaterialAveraging::none), + ExcMessage("The viscoelastic material model cannot be used with " + "material averaging. The variable 'Material averaging' " + "in the 'Material model' subsection must be set to 'none'.")); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + + + // Declare dependencies on solution variables + this->model_dependence.viscosity = NonlinearDependence::compositional_fields; + this->model_dependence.density = NonlinearDependence::temperature | NonlinearDependence::compositional_fields; + this->model_dependence.compressibility = NonlinearDependence::none; + this->model_dependence.specific_heat = NonlinearDependence::compositional_fields; + this->model_dependence.thermal_conductivity = NonlinearDependence::compositional_fields; + } + + + + template + void + Viscoelastic::create_additional_named_outputs (MaterialModel::MaterialModelOutputs &out) const + { + if (out.template get_additional_output >() == NULL) + { + const unsigned int n_points = out.viscosities.size(); + out.additional_outputs.push_back( + std_cxx11::shared_ptr > + (new MaterialModel::ElasticAdditionalOutputs (n_points))); + } + } + } +} + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { + ASPECT_REGISTER_MATERIAL_MODEL(Viscoelastic, + "viscoelastic", + "An implementation of a simple linear viscoelastic rheology that " + "only includes the deviatoric components of elasticity. Specifically, " + "the viscoelastic rheology only takes into account the elastic shear " + "strength (e.g., shear modulus), while the tensile and volumetric " + "strength (e.g., Young's and bulk modulus) are not considered. The " + "model is incompressible and allows specifying an arbitrary number " + "of compositional fields, where each field represents a different " + "rock type or component of the viscoelastic stress tensor. The stress " + "tensor in 2D and 3D, respectively, contains 3 or 6 components. The " + "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'. " + "\n\n " + "Expanding the model to include non-linear viscous flow (e.g., " + "diffusion/dislocation creep) and plasticity would produce a " + "constitutive relationship commonly referred to as partial " + "elastoviscoplastic (e.g., pEVP) in the geodynamics community. " + "While extensively discussed and applied within the geodynamics " + "literature, notable references include: " + "Moresi et al. (2003), J. Comp. Phys., v. 184, p. 476-497. " + "Gerya and Yuen (2007), Phys. Earth. Planet. Inter., v. 163, p. 83-105. " + "Gerya (2010), Introduction to Numerical Geodynamic Modeling. " + "Kaus (2010), Tectonophysics, v. 484, p. 36-47. " + "Choi et al. (2013), J. Geophys. Res., v. 118, p. 2429-2444. " + "Keller et al. (2013), Geophys. J. Int., v. 195, p. 1406-1442. " + "\n\n " + "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 " + "tracers, to track individual components of the viscoelastic stress " + "tensor. The material model will be udpated when an option to track " + "and calculate viscoelastic stresses with tracers is implemented. " + "\n\n " + "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: " + "$\\hat{D} = \\hat{D_{e}} + \\hat{D_{v}}$. " + "These terms further decompose into " + "$\\hat{D_{v}} = \\frac{\\tau}{2\\eta}$ and " + "$\\hat{D_{e}} = \\frac{\\overset{\\triangledown}{\\tau}}{2\\mu}$, where " + "$\\tau$ is the viscous deviatoric stress, $\\eta$ is the shear viscosity, " + "$\\mu$ is the shear modulus and $\\overset{\\triangledown}{\\tau}$ is the " + "Jaumann corotational stress rate. This later term (eqn. 24) contains the " + "time derivative of the deviatoric stress ($\\dot{\\tau}$) and terms that " + "account for material spin (e.g., rotation) due to advection: " + "$\\overset{\\triangledown}{\\tau} = \\dot{\\tau} + {\\tau}W -W\\tau$. " + "Above, $W$ is the material spin tensor (eqn. 25): " + "$W_{ij} = \\frac{1}{2} \\left (\\frac{\\partial V_{i}}{\\partial x_{j}} - " + "\\frac{\\partial V_{j}}{\\partial x_{i}} \\right )$. " + "\n\n " + "The Jaumann stress-rate can also be approximated using terms from the time " + "at the previous time step ($t$) and current time step ($t + \\Delta t_^{e}$): " + "$\\smash[t]{\\overset{\\triangledown}{\\tau}}^{t + \\Delta t^{e}} \\approx " + "\\frac{\\tau^{t + \\Delta t^{e} - \\tau^{t}}}{\\Delta t^{e}} - " + "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." + "\n\n " + "The formulation above allows rewriting the total rate of deformation (eqn. 29) as " + "$\\tau^{t + \\Delta t^{e}} = \\eta_{eff} \\left ( " + "2\\hat{D}^{t + \\triangle t^{e}} + \\frac{\\tau^{t}}{\\mu \\Delta t^{e}} + " + "\\frac{W^{t}\\tau^{t} - \\tau^{t}W^{t}}{\\mu} \\right )$. " + "\n\n " + "The effective viscosity (eqn. 28) is a function of the viscosity ($\\eta$), " + "elastic time step size ($\\Delta t^{e}$) and shear relaxation time " + "($ \\alpha = \\frac{\\eta}{\\mu} $): " + "$\\eta_{eff} = \\eta \\frac{\\Delta t^{e}}{\\Delta t^{e} + \\alpha}$ " + "The magnitude of the shear modulus thus controls how much the effective " + "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}$. " + "This force term is added onto the right-hand side force vector in the " + "system of equations. " + "\n\n " + "The value of each compositional field representing distinck rock types at a " + "point is interpreted to be a volume fraction of that rock type. If the sum of " + "the compositional field volume fractions is less than one, then the remainder " + "of the volume is assumed to be 'background material'." + "\n\n " + "Several model parameters (densities, elastic shear moduli, thermal expansivities, " + "thermal conductivies, specific heats) can be defined per-compositional field. " + "For each material parameter the 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 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, except density, which varies linearly with temperature according to the " + "thermal expansivity. " + "\n\n " + "When more than one compositional field is present at a point, they are averaged " + "arithmetically. An exception is viscosity, which may be averaged arithmetically, " + "harmonically, geometrically, or by selecting the viscosity of the composition field " + "with the greatest volume fraction.") + } +} diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 1a1ebb0da0e..a380d67db56 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -170,6 +170,9 @@ namespace aspect const MaterialModel::AdditionalMaterialOutputsStokesRHS *force = scratch.material_model_outputs.template get_additional_output >(); + const MaterialModel::ElasticOutputs + *elastic_outputs = scratch.material_model_outputs.template get_additional_output >(); + for (unsigned int q=0; qget_parameters().enable_additional_stokes_rhs) data.local_rhs(i) += (force->rhs_u[q] * scratch.phi_u[i] + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]) * JxW; + if (elastic_outputs != NULL && this->get_parameters().enable_elasticity) + data.local_rhs(i) += (double_contract(elastic_outputs->elastic_force[q],Tensor<2,dim>(scratch.grads_phi_u[i]))) + * JxW; + if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; j > (new MaterialModel::AdditionalMaterialOutputsStokesRHS (n_points))); } + Assert(!this->get_parameters().enable_additional_stokes_rhs || outputs.template get_additional_output >()->rhs_u.size() == n_points, ExcInternalError()); + + if ((this->get_parameters().enable_elasticity) && + outputs.template get_additional_output >() == NULL) + { + outputs.additional_outputs.push_back( + std_cxx11::shared_ptr > + (new MaterialModel::ElasticOutputs (n_points))); + } + + Assert(!this->get_parameters().enable_elasticity + || + outputs.template get_additional_output >()->elastic_force.size() + == n_points, ExcInternalError()); } diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 09156a2027e..cd05df7e13b 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -482,6 +482,10 @@ namespace aspect "of the Stokes equation. This feature is likely only used when implementing force " "vectors for manufactured solution problems and requires filling additional outputs " "of type AdditionalMaterialOutputsStokesRHS."); + prm.declare_entry ("Enable elasticity", "false", + Patterns::Bool (), + "Whether to include the additional elastic terms on the right-hand side of " + "the Stokes equation."); } prm.leave_subsection(); @@ -1163,6 +1167,7 @@ namespace aspect else AssertThrow(false, ExcNotImplemented()); enable_additional_stokes_rhs = prm.get_bool ("Enable additional Stokes RHS"); + enable_elasticity = prm.get_bool("Enable elasticity"); } prm.leave_subsection (); diff --git a/tests/viscoelastic_fixed_elastic_time_step.prm b/tests/viscoelastic_fixed_elastic_time_step.prm new file mode 100644 index 00000000000..4eda54b2d14 --- /dev/null +++ b/tests/viscoelastic_fixed_elastic_time_step.prm @@ -0,0 +1,57 @@ +# This test checks viscoelastic stress buildup with a +# constant strain rate and an elastic time step equal +# to a user-specified constant elastic time step. This +# fixed elastic time step is double the numerical time step. +# Viscoelastic stresses (sigma_xx, sigma_yy, sigma_xy) +# are tracked with compositional fields. The parameter +# file is derived from the 2D viscoelastic_stress_build-up +# benchmark with modifications to the resolution, end time, +# material properties and postprocessors. + +include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm + +# Global parameters +set End time = 1e3 +set Output directory = viscoelastic_fixed_elastic_time_step + +# Model geometry (100x100 km, 5 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 20 + set Y repetitions = 20 + set X extent = 100e3 + set Y extent = 100e3 + end +end + + +# Material model +subsection Material model + + set Model name = viscoelastic + + subsection Viscoelastic + set Densities = 2800 + set Viscosities = 1.e21 + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = true + set Fixed elastic time step = 2e3 + set Use stress averaging = false + set Viscosity averaging scheme = harmonic + end + +end + +# Post processing +subsection Postprocess + set List of postprocessors = velocity statistics, depth average + subsection Depth average + set Time between graphical output = 0 + set Number of zones = 5 + set List of output variables = composition + set Output format = txt + end +end + diff --git a/tests/viscoelastic_fixed_elastic_time_step/depth_average.txt b/tests/viscoelastic_fixed_elastic_time_step/depth_average.txt new file mode 100644 index 00000000000..e9b2f319d73 --- /dev/null +++ b/tests/viscoelastic_fixed_elastic_time_step/depth_average.txt @@ -0,0 +1,11 @@ +# time depth C_0 C_1 C_2 + 0 10000 0 0 0 + 0 30000 0 0 0 + 0 50000 0 0 0 + 0 70000 0 0 0 + 0 90000 0 0 0 + 1000 10000 7.42553e+06 -7.42553e+06 5.70638 + 1000 30000 7.64595e+06 -7.64595e+06 4.0183 + 1000 50000 7.646e+06 -7.646e+06 6.12554 + 1000 70000 7.64569e+06 -7.64569e+06 4.544 + 1000 90000 7.41715e+06 -7.41715e+06 0.23685 diff --git a/tests/viscoelastic_fixed_elastic_time_step/screen-output b/tests/viscoelastic_fixed_elastic_time_step/screen-output new file mode 100644 index 00000000000..fc76a8cdef3 --- /dev/null +++ b/tests/viscoelastic_fixed_elastic_time_step/screen-output @@ -0,0 +1,47 @@ + +Number of active cells: 400 (on 1 levels) +Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681) + +*** Timestep 0: t=0 years + Solving temperature system... 0 iterations. + Skipping stress_xx composition solve because RHS is zero. + Skipping stress_yy composition solve because RHS is zero. + Skipping stress_xy composition solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 27+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_fixed_elastic_time_step/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +*** Timestep 1: t=1000 years + Solving temperature system... 0 iterations. + Solving stress_xx system ... 11 iterations. + Solving stress_yy system ... 11 iterations. + Solving stress_xy system ... 12 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 17+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_fixed_elastic_time_step/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + diff --git a/tests/viscoelastic_fixed_elastic_time_step/statistics b/tests/viscoelastic_fixed_elastic_time_step/statistics new file mode 100644 index 00000000000..7d27fa0f55f --- /dev/null +++ b/tests/viscoelastic_fixed_elastic_time_step/statistics @@ -0,0 +1,18 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for composition solver 2 +# 11: Iterations for composition solver 3 +# 12: Iterations for Stokes solver +# 13: Velocity iterations in Stokes preconditioner +# 14: Schur complement iterations in Stokes preconditioner +# 15: RMS velocity (m/year) +# 16: Max. velocity (m/year) +0 0.000000000000e+00 0.000000000000e+00 400 3803 441 5043 0 0 0 0 27 28 28 2.57523020e-02 4.43529485e-02 +1 1.000000000000e+03 1.000000000000e+03 400 3803 441 5043 0 11 11 12 17 18 18 2.57523047e-02 4.43618279e-02 diff --git a/tests/viscoelastic_numerical_elastic_time_step.prm b/tests/viscoelastic_numerical_elastic_time_step.prm new file mode 100644 index 00000000000..fd5af6921fc --- /dev/null +++ b/tests/viscoelastic_numerical_elastic_time_step.prm @@ -0,0 +1,57 @@ +# This test checks viscoelastic stress buildup with a +# constant strain rate and an elastic time step equal +# to the numerical time step. The fixed elastic time +# step used during the first time step is equal to the +# numerical time step. Viscoelastic stresses (sigma_xx, +# sigma_yy, sigma_xy) are tracked with compositional +# fields. The parameter file is derived from the 2D +# viscoelastic_stress_build-up benchmark with modifications +# to the resolution, end time, material properties and +# postprocessors. + +include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm + +# Global parameters +set End time = 1e3 +set Output directory = viscoelastic_numerical_elastic_time_step + +# Model geometry (100x100 km, 5 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 20 + set Y repetitions = 20 + set X extent = 100e3 + set Y extent = 100e3 + end +end + + +# Material model +subsection Material model + + set Model name = viscoelastic + + subsection Viscoelastic + set Densities = 2800 + set Viscosities = 1.e21 + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = false + set Fixed elastic time step = 1e3 + set Use stress averaging = false + set Viscosity averaging scheme = harmonic + end + +end + +# Post processing +subsection Postprocess + set List of postprocessors = velocity statistics, depth average + subsection Depth average + set Time between graphical output = 0 + set Number of zones = 5 + set List of output variables = composition + set Output format = txt + end +end diff --git a/tests/viscoelastic_numerical_elastic_time_step/depth_average.txt b/tests/viscoelastic_numerical_elastic_time_step/depth_average.txt new file mode 100644 index 00000000000..4d3ea132c05 --- /dev/null +++ b/tests/viscoelastic_numerical_elastic_time_step/depth_average.txt @@ -0,0 +1,11 @@ +# time depth C_0 C_1 C_2 + 0 10000 0 0 0 + 0 30000 0 0 0 + 0 50000 0 0 0 + 0 70000 0 0 0 + 0 90000 0 0 0 + 1000 10000 4.60336e+06 -4.60336e+06 2.3362 + 1000 30000 4.74e+06 -4.74e+06 3.19024 + 1000 50000 4.74003e+06 -4.74003e+06 4.43544 + 1000 70000 4.73984e+06 -4.73984e+06 3.2894 + 1000 90000 4.59816e+06 -4.59816e+06 0.286356 diff --git a/tests/viscoelastic_numerical_elastic_time_step/screen-output b/tests/viscoelastic_numerical_elastic_time_step/screen-output new file mode 100644 index 00000000000..9a016bdb674 --- /dev/null +++ b/tests/viscoelastic_numerical_elastic_time_step/screen-output @@ -0,0 +1,47 @@ + +Number of active cells: 400 (on 1 levels) +Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681) + +*** Timestep 0: t=0 years + Solving temperature system... 0 iterations. + Skipping stress_xx composition solve because RHS is zero. + Skipping stress_yy composition solve because RHS is zero. + Skipping stress_xy composition solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 27+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_numerical_elastic_time_step/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +*** Timestep 1: t=1000 years + Solving temperature system... 0 iterations. + Solving stress_xx system ... 11 iterations. + Solving stress_yy system ... 11 iterations. + Solving stress_xy system ... 12 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 17+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_numerical_elastic_time_step/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + diff --git a/tests/viscoelastic_numerical_elastic_time_step/statistics b/tests/viscoelastic_numerical_elastic_time_step/statistics new file mode 100644 index 00000000000..2c4cfccab24 --- /dev/null +++ b/tests/viscoelastic_numerical_elastic_time_step/statistics @@ -0,0 +1,18 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for composition solver 2 +# 11: Iterations for composition solver 3 +# 12: Iterations for Stokes solver +# 13: Velocity iterations in Stokes preconditioner +# 14: Schur complement iterations in Stokes preconditioner +# 15: RMS velocity (m/year) +# 16: Max. velocity (m/year) +0 0.000000000000e+00 0.000000000000e+00 400 3803 441 5043 0 0 0 0 27 28 28 2.57523019e-02 4.43529487e-02 +1 1.000000000000e+03 1.000000000000e+03 400 3803 441 5043 0 11 11 12 17 18 18 2.57523061e-02 4.43639581e-02 diff --git a/tests/viscoelastic_stress_averaging.prm b/tests/viscoelastic_stress_averaging.prm new file mode 100644 index 00000000000..fd71aa82222 --- /dev/null +++ b/tests/viscoelastic_stress_averaging.prm @@ -0,0 +1,59 @@ +# This test checks viscoelastic stress buildup with a +# constant strain rate and an elastic time step equal +# to a user-specified constant elastic time step. This +# fixed elastic time step is double the numerical time step. +# Stress averaging is used to account for the difference +# between the numerical and elastic time step. Viscoelastic +# stresses (sigma_xx, sigma_yy, sigma_xy) are tracked +# with compositional fields. The parameter file is derived +# from the 2D viscoelastic_stress_build-up benchmark with +# modifications to the resolution, end time, material +# properties and postprocessors. + +include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm + +# Global parameters +set End time = 1e3 +set Output directory = viscoelastic_stress_averaging + +# Model geometry (100x100 km, 5 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 20 + set Y repetitions = 20 + set X extent = 100e3 + set Y extent = 100e3 + end +end + + +# Material model +subsection Material model + + set Model name = viscoelastic + + subsection Viscoelastic + set Densities = 2800 + set Viscosities = 1.e21 + set Elastic shear moduli = 1.e10 + set Use fixed elastic time step = true + set Fixed elastic time step = 2e3 + set Use stress averaging = true + set Viscosity averaging scheme = harmonic + end + +end + +# Post processing +subsection Postprocess + set List of postprocessors = velocity statistics, depth average + subsection Depth average + set Time between graphical output = 0 + set Number of zones = 5 + set List of output variables = composition + set Output format = txt + end +end + diff --git a/tests/viscoelastic_stress_averaging/depth_average.txt b/tests/viscoelastic_stress_averaging/depth_average.txt new file mode 100644 index 00000000000..954f6509b0e --- /dev/null +++ b/tests/viscoelastic_stress_averaging/depth_average.txt @@ -0,0 +1,11 @@ +# time depth C_0 C_1 C_2 + 0 10000 0 0 0 + 0 30000 0 0 0 + 0 50000 0 0 0 + 0 70000 0 0 0 + 0 90000 0 0 0 + 1000 10000 3.71277e+06 -3.71277e+06 2.85319 + 1000 30000 3.82297e+06 -3.82297e+06 2.00915 + 1000 50000 3.823e+06 -3.823e+06 3.06277 + 1000 70000 3.82285e+06 -3.82285e+06 2.272 + 1000 90000 3.70858e+06 -3.70858e+06 0.118425 diff --git a/tests/viscoelastic_stress_averaging/screen-output b/tests/viscoelastic_stress_averaging/screen-output new file mode 100644 index 00000000000..c0306780d50 --- /dev/null +++ b/tests/viscoelastic_stress_averaging/screen-output @@ -0,0 +1,47 @@ + +Number of active cells: 400 (on 1 levels) +Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681) + +*** Timestep 0: t=0 years + Solving temperature system... 0 iterations. + Skipping stress_xx composition solve because RHS is zero. + Skipping stress_yy composition solve because RHS is zero. + Skipping stress_xy composition solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 27+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_stress_averaging/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +*** Timestep 1: t=1000 years + Solving temperature system... 0 iterations. + Solving stress_xx system ... 11 iterations. + Solving stress_yy system ... 11 iterations. + Solving stress_xy system ... 12 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 16+0 iterations. + + Postprocessing: + RMS, max velocity: 0.0258 m/year, 0.0444 m/year + Writing depth average: output-viscoelastic_stress_averaging/depth_average.txt + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + diff --git a/tests/viscoelastic_stress_averaging/statistics b/tests/viscoelastic_stress_averaging/statistics new file mode 100644 index 00000000000..c0de1c73443 --- /dev/null +++ b/tests/viscoelastic_stress_averaging/statistics @@ -0,0 +1,18 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for composition solver 2 +# 11: Iterations for composition solver 3 +# 12: Iterations for Stokes solver +# 13: Velocity iterations in Stokes preconditioner +# 14: Schur complement iterations in Stokes preconditioner +# 15: RMS velocity (m/year) +# 16: Max. velocity (m/year) +0 0.000000000000e+00 0.000000000000e+00 400 3803 441 5043 0 0 0 0 27 28 28 2.57523020e-02 4.43529485e-02 +1 1.000000000000e+03 1.000000000000e+03 400 3803 441 5043 0 11 11 12 16 17 17 2.57523025e-02 4.43573871e-02