Only require old stresses for elasticity if necessary

benchmarks/viscoelastic_relaxation/stress_xx_over_time.py

@@ -63,7 +63,7 @@
   # Read in the time and the minimum xx and yy components of the viscoelastic stress,
   # which are stored on the fields ve_stress_xx and ve_stress_yy.
   # The correct columns are selected with usecols.
-  time,stress_xx_min = np.genfromtxt(path, comments='#', usecols=(1,15), unpack=True)
+  time,stress_xx_min = np.genfromtxt(path, comments='#', usecols=(1,18), unpack=True)
 
   # Plot the stress elements in MPa against time in ky in
   # categorical batlow colors.

benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm

@@ -66,9 +66,9 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress, stress, stress, stress, stress, stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress, stress, stress
 end
 
 # Spatial domain of different compositional fields
@@ -77,7 +77,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 20e6; -20e6; 0; 20e6; -20e6; 0
+    set Function expression = 20e6; -20e6; 0
   end
 end
 

benchmarks/viscoelastic_stress_build_up/viscoelastic_stress_build_up.prm

@@ -114,9 +114,9 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress, stress, stress, stress, stress, stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress, stress, stress
 end
 
 # Spatial domain of different compositional fields
@@ -126,7 +126,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0
   end
 end
 

include/aspect/material_model/rheology/elasticity.h

@@ -233,6 +233,16 @@ namespace aspect
                                               const double viscosity_pre_yield,
                                               const double shear_modulus) const;
 
+          /**
+          * Whether the unrotated viscoelastic stress is required in order to update the
+          * stress values. This is the case if the elastic time step is not equal to the
+          * computational time step of ASPECT, either because a fixed elastic time step
+          * was requested, or because a fixed stabilization factor was chosen for the
+          * elastic time scale.
+          */
+          bool
+          require_unrotated_viscoelastic_stress() const;
+
           /**
            * Compute the elastic time step.
            */

source/material_model/rheology/visco_plastic.cc

@@ -126,9 +126,11 @@ namespace aspect
 
             // The old stresses are only changed in the operator splitting step and have been advected into
             // the current timestep. They are stored in the second set of n_independent_components.
-            stress_old = (Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]],
-                                                                       &in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]]+SymmetricTensor<2, dim>::n_independent_components));
-
+            if (stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components)
+              {
+                stress_old = (Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]],
+                                                                           &in.composition[i][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components]]+SymmetricTensor<2, dim>::n_independent_components));
+              }
 
             // Average the compositional contributions to elastic_shear_moduli here and use
             // a volume-averaged shear modulus in the loop over the compositions below.
@@ -625,7 +627,7 @@ namespace aspect
 
         if (this->get_parameters().enable_elasticity)
           {
-            // Set a mask (false) for the 2*n_independent_components fields representing the viscoelastic
+            // Set a mask (false) for the fields representing the viscoelastic
             // stress tensor components so that they are excluded from the volume fraction computation.
             const std::vector<unsigned int> &stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress);
             for (unsigned int field_index: stress_field_indices)

source/particle/property/elastic_stress.cc

@@ -65,8 +65,9 @@ namespace aspect
 
         // Get the indices of those compositions that correspond to stress tensor elements.
         stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress);
-        AssertThrow((stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components),
-                    ExcMessage("The number of stress tensor element fields in the 'elastic stress' plugin does not equal twice the number of independent components."));
+        AssertThrow((stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components ||
+                     stress_field_indices.size() == SymmetricTensor<2,dim>::n_independent_components),
+                    ExcMessage("The number of stress tensor element fields in the 'elastic stress' plugin does not equal the number of expected components."));
 
         // Get the indices of all compositions that do not correspond to stress tensor elements.
         std::vector<unsigned int> all_field_indices(this->n_compositional_fields());
@@ -283,23 +284,26 @@ namespace aspect
             data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yz")));
           }
 
-        data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xx_old")));
+        if (stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components)
+          {
+            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xx_old")));
 
-        data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yy_old")));
+            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yy_old")));
 
-        if (dim == 2)
-          {
-            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old")));
-          }
-        else if (dim == 3)
-          {
-            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_zz_old")));
+            if (dim == 2)
+              {
+                data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old")));
+              }
+            else if (dim == 3)
+              {
+                data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_zz_old")));
 
-            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old")));
+                data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xy_old")));
 
-            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xz_old")));
+                data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_xz_old")));
 
-            data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yz_old")));
+                data.push_back(this->get_initial_composition_manager().initial_composition(position,this->introspection().compositional_index_for_name("ve_stress_yz_old")));
+              }
           }
       }
 
@@ -310,6 +314,8 @@ namespace aspect
       ElasticStress<dim>::update_particle_properties(const ParticleUpdateInputs<dim> &inputs,
                                                      typename ParticleHandler<dim>::particle_iterator_range &particles) const
       {
+        const unsigned int n_total_stress_components = stress_field_indices.size();
+
         unsigned int p = 0;
         for (auto &particle: particles)
           {
@@ -329,7 +335,7 @@ namespace aspect
             for (const unsigned int &n : non_stress_field_indices)
               material_inputs.composition[0][n] = inputs.solution[p][this->introspection().component_indices.compositional_fields[n]];
             // For the stress composition we use the ve_stress_* stored on the particles.
-            for (unsigned int n = 0; n < 2*SymmetricTensor<2,dim>::n_independent_components; ++n)
+            for (unsigned int n = 0; n < n_total_stress_components; ++n)
               material_inputs.composition[0][stress_field_indices[n]] = particle.get_properties()[this->data_position + n];
 
             Tensor<2,dim> grad_u;
@@ -391,19 +397,22 @@ namespace aspect
             property_information.emplace_back("ve_stress_yz",1);
           }
 
-        property_information.emplace_back("ve_stress_xx_old",1);
-        property_information.emplace_back("ve_stress_yy_old",1);
-
-        if (dim == 2)
-          {
-            property_information.emplace_back("ve_stress_xy_old",1);
-          }
-        else if (dim == 3)
+        if  (stress_field_indices.size() == 2*SymmetricTensor<2,dim>::n_independent_components)
           {
-            property_information.emplace_back("ve_stress_zz_old",1);
-            property_information.emplace_back("ve_stress_xy_old",1);
-            property_information.emplace_back("ve_stress_xz_old",1);
-            property_information.emplace_back("ve_stress_yz_old",1);
+            property_information.emplace_back("ve_stress_xx_old",1);
+            property_information.emplace_back("ve_stress_yy_old",1);
+
+            if (dim == 2)
+              {
+                property_information.emplace_back("ve_stress_xy_old",1);
+              }
+            else if (dim == 3)
+              {
+                property_information.emplace_back("ve_stress_zz_old",1);
+                property_information.emplace_back("ve_stress_xy_old",1);
+                property_information.emplace_back("ve_stress_xz_old",1);
+                property_information.emplace_back("ve_stress_yz_old",1);
+              }
           }
 
         return property_information;

source/simulator/solver_schemes.cc

@@ -405,48 +405,27 @@ namespace aspect
     // of the same stress tensor in the computation of the initial residual
     // for the fields that belong to that tensor. In other words, we compute an
     // averaged initial residual using those fields that belong to the ve_stress tensor.
-    // The ve_stress_old values are not independent components of the solution vector,
-    // so we do not need to compute a residual for them and set their residual to zero.
     // TODO Is this residual calculation representative of a second order tensor?
-    if (parameters.enable_elasticity == true)
+    if (parameters.enable_elasticity == true && residual)
       {
+        const std::vector<unsigned int> &stress_field_indices = introspection.get_indices_for_fields_of_type(CompositionalFieldDescription::stress);
+        const double n_stress_fields = stress_field_indices.size();
+
+        AssertThrow((n_stress_fields == 2*SymmetricTensor<2,dim>::n_independent_components ||
+                     n_stress_fields == SymmetricTensor<2,dim>::n_independent_components),
+                    ExcMessage("The number of stress tensor element fields does not equal the number of expected components."));
+
         double stress_initial_residual = 0.0;
-        std::vector<unsigned int> stress_indices;
-        std::vector<unsigned int> old_stress_indices;
-        stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xx"));
-        stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yy"));
-        old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xx_old"));
-        old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yy_old"));
-        if (dim == 2)
-          {
-            stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy"));
-            old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy_old"));
-          }
-        else if (dim == 3)
-          {
-            stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_zz"));
-            stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy"));
-            stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xz"));
-            stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yz"));
-            old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_zz_old"));
-            old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xy_old"));
-            old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_xz_old"));
-            old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yz_old"));
-          }
 
+        for (unsigned int c=0; c<SymmetricTensor<2,dim>::n_independent_components; ++c)
+          stress_initial_residual += system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm();
 
-        if (residual)
-          {
-            const double n_stress_fields = stress_indices.size();
-            for (auto &c : stress_indices)
-              stress_initial_residual += system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm() / n_stress_fields;
-
-            // Overwrite the initial residual
-            for (auto &c : stress_indices)
-              (*residual)[c] = stress_initial_residual;
-            for (auto &c : old_stress_indices)
-              (*residual)[c] = 0.;
-          }
+        stress_initial_residual /= static_cast<double>(SymmetricTensor<2,dim>::n_independent_components);
+
+        // The ve_stress_old values if they exist are not independent components of the solution vector,
+        // so we do not compute a residual for them and set their residual to zero.
+        for (unsigned int c=0; c<n_stress_fields; ++c)
+          (*residual)[stress_field_indices[c]] = (c<SymmetricTensor<2,dim>::n_independent_components) ? stress_initial_residual : 0.;
       }
 
 

tests/boundary_traction_function_cartesian_free_surface.prm

@@ -90,9 +90,9 @@ subsection Boundary traction model
 end
 
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress,       stress,       stress,       stress,           stress,           stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress,       stress,       stress
 end
 
 subsection Initial composition model
@@ -101,7 +101,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0
   end
 end
 

tests/boundary_traction_function_cartesian_free_surface_mesh_deformation_degree.prm

@@ -94,9 +94,9 @@ subsection Boundary traction model
 end
 
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress,       stress,       stress,       stress,           stress,           stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress,       stress,       stress
 end
 
 subsection Initial composition model
@@ -105,7 +105,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0
   end
 end
 

tests/boundary_traction_function_spherical_free_surface.prm

@@ -91,9 +91,9 @@ subsection Boundary traction model
 end
 
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress,       stress,       stress,       stress,           stress,           stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress,       stress,       stress
 end
 
 subsection Initial composition model
@@ -102,7 +102,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0
   end
 end
 

tests/free_surface_timestep_repeat.prm

@@ -108,9 +108,9 @@ end
 
 # Number and name of compositional fields
 subsection Compositional fields
-  set Number of fields = 6
-  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old
-  set Types of fields  = stress,       stress,       stress,       stress,           stress,           stress
+  set Number of fields = 3
+  set Names of fields  = ve_stress_xx, ve_stress_yy, ve_stress_xy
+  set Types of fields  = stress,       stress,       stress
 end
 
 # Spatial domain of different compositional fields
@@ -120,7 +120,7 @@ subsection Initial composition model
   subsection Function
     set Variable names      = x,y
     set Function constants  =
-    set Function expression = 0; 0; 0; 0; 0; 0;
+    set Function expression = 0; 0; 0
   end
 end