From f135f60d8a305cd79367254f5e19285ed96b4b60 Mon Sep 17 00:00:00 2001 From: naliboff Date: Sat, 13 May 2017 08:37:29 -0400 Subject: [PATCH 1/5] Add option to use full finite strain tensor for strain weakening --- include/aspect/material_model/visco_plastic.h | 1 + source/material_model/visco_plastic.cc | 129 +++++++++++++++--- ...ld_strain_weakening_full_strain_tensor.prm | 116 ++++++++++++++++ 3 files changed, 230 insertions(+), 16 deletions(-) create mode 100644 tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index 950e83a83a2..d46fb21c54b 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -186,6 +186,7 @@ namespace aspect const YieldScheme &yield_type) const; bool use_strain_weakening; + bool use_finite_strain_tensor; std::vector start_strain_weakening_intervals; std::vector end_strain_weakening_intervals; std::vector viscous_strain_weakening_factors; diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index c37564f22b1..ca1d586be1c 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -41,7 +41,17 @@ namespace aspect // assign compositional fields associated with strain a value of 0 if (use_strain_weakening == true) - x_comp[0] = 0.0; + { + if (use_finite_strain_tensor == false) + { + x_comp[0] = 0.0; + } + else + { + for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) + x_comp[i] = 0.; + } + } //sum the compositional fields for normalization purposes double sum_composition = 0.0; @@ -200,11 +210,25 @@ namespace aspect double coh = cohesions[j]; // Strain weakening + double strain_ii = 0.; if (use_strain_weakening == true) { - // Constrain the second strain invariant of the previous timestep - const double strain_ii = std::max(std::min(composition[0],end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]); + // Calculate and/or constrain the strain invariant of the previous timestep + if ( use_finite_strain_tensor == true ) + { + // Calculate second invarinat of left stretching tensor "L" + Tensor<2,dim> strain; + for (unsigned int q = 0; q < Tensor<2,dim>::n_independent_components ; ++q) + strain[Tensor<2,dim>::unrolled_to_component_indices(q)] = composition[q]; + const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) ); + strain_ii = std::fabs(second_invariant(L)); + } + else + { + // Here the compositional field already contains the finite strain invariant magnitude + strain_ii = std::max(std::min(composition[0],end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]); + } // Linear strain weakening of cohesion and internal friction angle between specified strain values const double strain_fraction = ( strain_ii - start_strain_weakening_intervals[j] ) / @@ -274,6 +298,10 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { + + + + // Loop through points for (unsigned int i=0; i < in.temperature.size(); ++i) { const double temperature = in.temperature[i]; @@ -347,14 +375,56 @@ namespace aspect out.reaction_terms[i][c] = 0.0; // If strain weakening is used, overwrite the first reaction term, // which represents the second invariant of the strain tensor - if (use_strain_weakening && this->get_timestep_number() > 0) + double edot_ii = 0.; + double e_ii = 0.; + if (use_strain_weakening == true && use_finite_strain_tensor == false && this->get_timestep_number() > 0) { - const double edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate); - - // Update reaction term - out.reaction_terms[i][0] = edot_ii*this->get_timestep(); + edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate); + e_ii = edot_ii*this->get_timestep(); + // Update reaction term + out.reaction_terms[i][0] = e_ii; } } + + // We need the velocity gradient for the finite strain (they are not included in material model inputs), + // so we get them from the finite element. + if (in.cell && use_finite_strain_tensor == true && this->get_timestep_number() > 0) + { + const QGauss quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+1); + FEValues fe_values (this->get_mapping(), + this->get_fe(), + quadrature_formula, + update_gradients); + + std::vector > velocity_gradients (quadrature_formula.size(), Tensor<2,dim>()); + + fe_values.reinit (*in.cell); + fe_values[this->introspection().extractors.velocities].get_function_gradients (this->get_solution(), + velocity_gradients); + + // Assign the strain components to the compositional fields reaction terms. + // If there are too many fields, we simply fill only the first fields with the + // existing strain rate tensor components. + for (unsigned int q=0; q < in.position.size(); ++q) + { + // Convert the compositional fields into the tensor quantity they represent. + Tensor<2,dim> strain; + for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) + { + strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i]; + } + + // Compute the strain accumulated in this timestep. + const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain); + + // Output the strain increment component-wise to its respective compositional field's reaction terms. + for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) + { + out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)]; + } + } + + } } template @@ -429,6 +499,10 @@ namespace aspect Patterns::Bool (), "Apply strain weakening to viscosity, cohesion and internal angle " "of friction based on accumulated finite strain. Units: None"); + prm.declare_entry ("Use finite strain tensor", "false", + Patterns::Bool (), + "Track and use the full finite strain tensor for strain weakening. " + "Units: None"); prm.declare_entry ("Start strain weakening intervals", "0.", Patterns::List(Patterns::Double(0)), "List of strain weakening interval initial strains " @@ -568,6 +642,9 @@ namespace aspect //increment by one for background: const unsigned int n_fields = this->n_compositional_fields() + 1; + // number of required compositional fields for full finite strain tensor + const unsigned int s = Tensor<2,dim>::n_independent_components; + prm.enter_subsection("Material model"); { prm.enter_subsection ("Visco Plastic"); @@ -603,7 +680,10 @@ namespace aspect if (use_strain_weakening) AssertThrow(this->n_compositional_fields() >= 1, ExcMessage("There must be at least one compositional field. ")); - + use_finite_strain_tensor = prm.get_bool ("Use finite strain tensor"); + if (use_finite_strain_tensor) + AssertThrow(this->n_compositional_fields() >= s, + ExcMessage("There must be enough compositional fields to track all components of the finite strain tensor (4 in 2D, 9 in 3D). ")); start_strain_weakening_intervals = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Start strain weakening intervals"))), n_fields, "Start strain weakening intervals"); @@ -758,12 +838,29 @@ namespace aspect "See, for example, Thieulot, C. (2011), PEPI 188, pp. 47-68. " "\n\n" "The user has the option to linearly reduce the cohesion and " - "internal friction angle as a function of the finite strain invariant. " - "The finite strain invariant may be calculated through particles or the compositional " - "field finite strain invariant plugin (see \\cite{buiter2008dissipation} benchmark). " - "In either case, the user specifies a compositional field for the finite " - "strain invariant, which must be the first listed compositional field " - "in the parameter file." + "internal friction angle as a function of the finite strain magnitude. " + "The finite strain invariant or full strain tensor is calculated through " + "compositional fields within the material model. This implementation is " + "identical to the compositional field finite strain plugin and cookbook " + "described in the manual (author: Gassmoeller). If the user selects to track " + "the finite strain invariant ($e_ii$), a single compositional field tracks " + "the value derived from $e_ii^t = (e_ii)^(t-1) + edot_ii*dt, where $t$ and $t-1$ " + "are the current and prior time steps, $edot_ii$ is the second invariant of the " + "strain rate tensor and $dt$ is the time step size. In the case of the " + "full strain tensor $F$, the finite strain magnitude is derived from the " + "second invariant of the symmetric stretching tensor $L$, where " + "$L = F * [F]^T$. In either case, the user specifies a single compositional " + "field for the finite strain invariant or multiple fields (4 in 2D, 9 in 3D) " + "for the finite strain tensor. These field(s) must be the first lised " + "compositional fields in the paramter file. Note that one or more of the finite strain " + "tensor components must be assigned a non-zero value intially. This value can be " + "be quite small (ex: 1.e-8), but still non-zero. While the option to track and use " + "the full finite strain tensor exists, tracking the associated compositional " + "is computationally expensive in 3D. Similarly, the finite strain magnitudes " + "may in fact decrease if the orientation of the deformation field switches " + "through time. Consequently, the ideal solution is track the finite strain " + "invariant (single compositional) field within the material and track " + "the full finite strain tensor through tracers." "" "\n\n" "Viscous stress may also be limited by a non-linear stress limiter " @@ -775,7 +872,7 @@ namespace aspect "reference strain rate, $\\varepsilon_{ii}$ is the strain rate " "and $n_y$ is the stress limiter exponent. The yield stress, " "$\\tau_y$, is defined through the Drucker Prager yield criterion " - "formulation. This method of limiting viscous stress has been used " + "formulation. This method of limiting viscous stress has been used " "in various forms within the geodynamic literature, including " "Christensen (1992), JGR, 97(B2), pp. 2015-2036; " "Cizkova and Bina (2013), EPSL, 379, pp. 95-103; " diff --git a/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm new file mode 100644 index 00000000000..76b10022b14 --- /dev/null +++ b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm @@ -0,0 +1,116 @@ +# Global parameters +set Dimension = 2 +set Start time = 0 +set End time = 1e2 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = iterated Stokes +set Max nonlinear iterations = 1 +set Number of cheap Stokes solver steps = 0 +set Output directory = visco_plastic_yield_strain_weakening_full_strain_tensor +set Timing output frequency = 1 + +# Model geometry (100x100 km, 10 km spacing) +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 + + +# Boundary classifications (fixed T boundaries, prescribed velocity) +subsection Model settings + set Fixed temperature boundary indicators = bottom, top, left, right + set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function +end + +# Velocity on boundaries characterized by functions +subsection Boundary velocity model + subsection Function + set Variable names = x,y + set Function constants = m=0.0005, year=1 + set Function expression = if (x<50e3 , -1*m/year, 1*m/year); if (y<50e3 , 1*m/year, -1*m/year); + end +end + +# Temperature boundary and initial conditions +subsection Boundary temperature model + set Model name = box + subsection Box + set Bottom temperature = 273 + set Left temperature = 273 + set Right temperature = 273 + set Top temperature = 273 + end +end +subsection Initial temperature model + set Model name = function + subsection Function + set Function expression = 273 + end +end + +# Compositional fields used to track finite strain +subsection Compositional fields + set Number of fields = 4 + set Names of fields = s11, s12, s21, s22 +end + +# All compositional fields have an initial value of 0 +subsection Initial composition model + set Model name = function + subsection Function + set Variable names = x,y + set Function expression = 0.5; 0; 0; 0.5; + end +end + +# Boundary composition specification +subsection Boundary composition model + set Model name = initial composition +end + +# Material model (values for background material and strain fields) +subsection Material model + set Model name = visco plastic + subsection Visco Plastic + set Reference strain rate = 1.e-16 + set Viscous flow law = dislocation + set Prefactors for dislocation creep = 5.e-23 + set Stress exponents for dislocation creep = 1.0 + set Activation energies for dislocation creep = 0. + set Activation volumes for dislocation creep = 0. + set Yield mechanism = drucker + set Angles of internal friction = 0. + set Cohesions = 1.e6 + set Use strain weakening = true + set Use finite strain tensor = true + set Start strain weakening intervals = 0. + set End strain weakening intervals = 1.0 + set Cohesion strain weakening factors = 0.5 + set Friction strain weakening factors = 0.5 + end +end + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 10.0 + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = velocity statistics, mass flux statistics +end From e29c9f55c63925529d259101f41992f8aa3a0087 Mon Sep 17 00:00:00 2001 From: naliboff Date: Sat, 13 May 2017 20:59:34 -0400 Subject: [PATCH 2/5] Address comments, fix astyle --- source/material_model/visco_plastic.cc | 22 +++++++++---------- ...ld_strain_weakening_full_strain_tensor.prm | 2 +- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index ca1d586be1c..d53cf4603fc 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -217,12 +217,12 @@ namespace aspect // Calculate and/or constrain the strain invariant of the previous timestep if ( use_finite_strain_tensor == true ) { - // Calculate second invarinat of left stretching tensor "L" + // Calculate second invariant of left stretching tensor "L" Tensor<2,dim> strain; for (unsigned int q = 0; q < Tensor<2,dim>::n_independent_components ; ++q) strain[Tensor<2,dim>::unrolled_to_component_indices(q)] = composition[q]; - const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) ); - strain_ii = std::fabs(second_invariant(L)); + const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) ); + strain_ii = std::max(std::min(std::fabs(second_invariant(L)),end_strain_weakening_intervals[j]),start_strain_weakening_intervals[j]); } else { @@ -379,16 +379,16 @@ namespace aspect double e_ii = 0.; if (use_strain_weakening == true && use_finite_strain_tensor == false && this->get_timestep_number() > 0) { - edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate); - e_ii = edot_ii*this->get_timestep(); - // Update reaction term + edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate); + e_ii = edot_ii*this->get_timestep(); + // Update reaction term out.reaction_terms[i][0] = e_ii; } } // We need the velocity gradient for the finite strain (they are not included in material model inputs), // so we get them from the finite element. - if (in.cell && use_finite_strain_tensor == true && this->get_timestep_number() > 0) + if (in.cell && use_strain_weakening == true && use_finite_strain_tensor == true && this->get_timestep_number() > 0) { const QGauss quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+1); FEValues fe_values (this->get_mapping(), @@ -416,7 +416,7 @@ namespace aspect // Compute the strain accumulated in this timestep. const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain); - + // Output the strain increment component-wise to its respective compositional field's reaction terms. for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) { @@ -424,7 +424,7 @@ namespace aspect } } - } + } } template @@ -849,10 +849,10 @@ namespace aspect "strain rate tensor and $dt$ is the time step size. In the case of the " "full strain tensor $F$, the finite strain magnitude is derived from the " "second invariant of the symmetric stretching tensor $L$, where " - "$L = F * [F]^T$. In either case, the user specifies a single compositional " + "$L = F * [F]^T$. The user must specify a single compositional " "field for the finite strain invariant or multiple fields (4 in 2D, 9 in 3D) " "for the finite strain tensor. These field(s) must be the first lised " - "compositional fields in the paramter file. Note that one or more of the finite strain " + "compositional fields in the parameter file. Note that one or more of the finite strain " "tensor components must be assigned a non-zero value intially. This value can be " "be quite small (ex: 1.e-8), but still non-zero. While the option to track and use " "the full finite strain tensor exists, tracking the associated compositional " diff --git a/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm index 76b10022b14..c18e90c53a5 100644 --- a/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm +++ b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor.prm @@ -66,7 +66,7 @@ subsection Compositional fields set Names of fields = s11, s12, s21, s22 end -# All compositional fields have an initial value of 0 +# The compositional fields have an initial value 0 or 0.5 subsection Initial composition model set Model name = function subsection Function From 2df6e70591846081ac871da08a6848bad6d3ce13 Mon Sep 17 00:00:00 2001 From: naliboff Date: Sun, 14 May 2017 22:08:52 -0400 Subject: [PATCH 3/5] Add fe_values access --- source/material_model/visco_plastic.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index d53cf4603fc..b5a8f67b929 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -20,6 +20,7 @@ #include #include +#include using namespace dealii; namespace aspect From 12d71cd3d3bc04570b3d98a354eb86b1da101774 Mon Sep 17 00:00:00 2001 From: naliboff Date: Tue, 16 May 2017 11:37:13 -0400 Subject: [PATCH 4/5] Add test output --- .../screen-output | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 tests/visco_plastic_yield_strain_weakening_full_strain_tensor/screen-output diff --git a/tests/visco_plastic_yield_strain_weakening_full_strain_tensor/screen-output b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor/screen-output new file mode 100644 index 00000000000..2fa23f2da36 --- /dev/null +++ b/tests/visco_plastic_yield_strain_weakening_full_strain_tensor/screen-output @@ -0,0 +1,52 @@ + +Number of active cells: 100 (on 1 levels) +Number of degrees of freedom: 3,208 (882+121+441+441+441+441+441) + +*** Timestep 0: t=0 years + Solving temperature system... 0 iterations. + Solving s11 system ... 0 iterations. + Skipping s12 composition solve because RHS is zero. + Skipping s21 composition solve because RHS is zero. + Solving s22 system ... 0 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+2 iterations. + Relative Stokes residual after nonlinear iteration 1: 1 + + + Postprocessing: + RMS, max velocity: 0.000408 m/year, 0.000691 m/year + Mass fluxes through boundary parts: 1.651e+05 kg/yr, 1.651e+05 kg/yr, -1.651e+05 kg/yr, -1.651e+05 kg/yr + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +*** Timestep 1: t=100 years + Solving temperature system... 0 iterations. + Solving s11 system ... 6 iterations. + Solving s12 system ... 12 iterations. + Solving s21 system ... 12 iterations. + Solving s22 system ... 6 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+0 iterations. + Relative Stokes residual after nonlinear iteration 1: 6.62146e-15 + + Postprocessing: + RMS, max velocity: 0.000408 m/year, 0.000691 m/year + Mass fluxes through boundary parts: 1.651e+05 kg/yr, 1.651e+05 kg/yr, -1.651e+05 kg/yr, -1.651e+05 kg/yr + + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + From 23cbf2a8ebe7476d6165b0225748649b1e235bdf Mon Sep 17 00:00:00 2001 From: naliboff Date: Tue, 16 May 2017 12:09:36 -0400 Subject: [PATCH 5/5] Add Dannberg as co-author for finite strain cookbook and plugin --- source/material_model/visco_plastic.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index b5a8f67b929..877f7f5181f 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -843,7 +843,7 @@ namespace aspect "The finite strain invariant or full strain tensor is calculated through " "compositional fields within the material model. This implementation is " "identical to the compositional field finite strain plugin and cookbook " - "described in the manual (author: Gassmoeller). If the user selects to track " + "described in the manual (author: Gassmoeller, Dannberg). If the user selects to track " "the finite strain invariant ($e_ii$), a single compositional field tracks " "the value derived from $e_ii^t = (e_ii)^(t-1) + edot_ii*dt, where $t$ and $t-1$ " "are the current and prior time steps, $edot_ii$ is the second invariant of the "