Pass shear modulus to fill_elastic_outputs

anne-glerum · Mar 11, 2024 · 3d50613 · 3d50613
1 parent 6828a2b
commit 3d50613
Show file tree

Hide file tree

Showing 4 changed files with 6 additions and 9 deletions.
diff --git a/include/aspect/material_model/rheology/elasticity.h b/include/aspect/material_model/rheology/elasticity.h
@@ -124,6 +124,7 @@ namespace aspect
            */
           void fill_elastic_outputs(const unsigned int point_index,
                                     const MaterialModel::MaterialModelInputs<dim> &in,
+                                    const double average_elastic_shear_modulus,
                                     MaterialModel::MaterialModelOutputs<dim> &out) const;
 
           /**

diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc
@@ -590,16 +590,13 @@ namespace aspect
       Elasticity<dim>::
       fill_elastic_outputs(const unsigned int q,
                            const MaterialModel::MaterialModelInputs<dim> &in,
+                           const double average_elastic_shear_modulus,
                            MaterialModel::MaterialModelOutputs<dim> &out) const
       {
         MaterialModel::ElasticOutputs<dim> *elastic_out = out.template get_additional_output<MaterialModel::ElasticOutputs<dim>>();
 
         if (elastic_out != nullptr)
           {
-
-            ElasticAdditionalOutputs<dim> *elastic_additional_out = out.template get_additional_output<ElasticAdditionalOutputs<dim>>();
-            AssertThrow (elastic_additional_out != nullptr, ExcMessage("To compute the ElasticOutputs viscous_dissipation, the ElasticAdditionalOutputs are required."));
-
             const SymmetricTensor<2, dim> deviatoric_strain_rate = in.strain_rate[q];
 
             SymmetricTensor<2, dim> stress_0;
@@ -615,8 +612,7 @@ namespace aspect
             // Retrieve the timestep ratio dtc/dte and the elastic viscosity,
             // only two material models support elasticity.
             const double timestep_ratio = calculate_timestep_ratio();
-            const double shear_modulus = elastic_additional_out->elastic_shear_moduli[q];
-            const double elastic_viscosity = calculate_elastic_viscosity(shear_modulus);
+            const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_modulus);
 
             // Apply the stress update to get the total stress of timestep t.
             // Both the viscous and the elastic viscosity have been scaled with the timestep ratio.
@@ -630,7 +626,7 @@ namespace aspect
             // Assume incompressibility. If compressible,
             // visco_plastic_strain_rate = visco_plastic_strain_rate -
             //                             1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor<dim>();
-            const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - ((stress - stress_0) / (2. * dtc * shear_modulus));
+            const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - ((stress - stress_0) / (2. * dtc * average_elastic_shear_modulus));
 
             elastic_out->viscous_dissipation[q] = stress * visco_plastic_strain_rate;
           }

diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
@@ -304,7 +304,7 @@ namespace aspect
               // At the moment, that means compute the viscous dissipation based on the viscoplastic strain rate.
               if (ElasticOutputs<dim> *elastic_out = out.template get_additional_output<ElasticOutputs<dim>>())
                 {
-                  rheology->elastic_rheology.fill_elastic_outputs(i, in, out);
+                  rheology->elastic_rheology.fill_elastic_outputs(i, in, average_elastic_shear_moduli[i], out);
                 }
             }
         }

diff --git a/source/material_model/viscoelastic.cc b/source/material_model/viscoelastic.cc
@@ -88,7 +88,7 @@ namespace aspect
           if (ElasticAdditionalOutputs<dim> *elastic_out = out.template get_additional_output<ElasticAdditionalOutputs<dim>>())
             {
               elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i];
-              elastic_rheology.fill_elastic_outputs(i, in, out);
+              elastic_rheology.fill_elastic_outputs(i, in, average_elastic_shear_moduli[i], out);
             }
         }