Simplify order ve viscosity calculation

geodynamics · Jul 26, 2022 · f02d6d6 · f02d6d6
1 parent 4f68443
commit f02d6d6
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 9 deletions.
diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc
@@ -643,15 +643,14 @@ namespace aspect
         // which is equal to 0.5 * stress / viscosity.
         // in.composition contains $\tau^{0}$, so that we can compute
         // $\dot\varepsilon_T + \frac{\tau^{0adv}}{2 \eta_{el}}$.
-        // The input parameter viscosity_pre_yield has not yet been scaled with the timestep ratio dtc/dte.
         // During assembly in timestep 0, get_timestep() returns 0. Therefore we have to make an estimated guess
         // using the maximum timestep parameters capped by the elastic timestep.
         double dtc = this->get_timestep();
         if (this->get_timestep_number() == 0 && this->get_timestep() == 0)
           dtc = std::min(std::min(this->get_parameters().maximum_time_step, this->get_parameters().maximum_first_time_step), elastic_timestep());
         const double timestep_ratio = dtc / elastic_timestep();
         const double elastic_viscosity = timestep_ratio * calculate_elastic_viscosity(shear_modulus);
-        const double creep_viscosity = timestep_ratio *  calculate_viscoelastic_viscosity(viscosity_pre_yield, shear_modulus);
+        const double creep_viscosity = viscosity_pre_yield; //timestep_ratio *  calculate_viscoelastic_viscosity(viscosity_pre_yield, shear_modulus);
 
         const SymmetricTensor<2, dim>
         edot_deviator = deviator(strain_rate) + 0.5 * stress_0_advected / elastic_viscosity

diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc
@@ -256,6 +256,14 @@ namespace aspect
 
             if (use_elasticity)
               {
+
+                // Step 3a: calculate viscoelastic (effective) viscosity
+                // Estimate the timestep size when in timestep 0
+                double dtc = this->get_timestep();
+                if (this->get_timestep_number() == 0 && this->get_timestep() == 0)
+                  dtc = std::min(std::min(this->get_parameters().maximum_time_step, this->get_parameters().maximum_first_time_step), elastic_rheology.elastic_timestep());
+                viscosity_pre_yield = dtc / elastic_rheology.elastic_timestep() * elastic_rheology.calculate_viscoelastic_viscosity(viscosity_pre_yield, elastic_shear_modulus);
+
                 if (use_reference_strainrate == true)
                   current_edot_ii = ref_strain_rate;
                 else
@@ -274,13 +282,6 @@ namespace aspect
                                                min_strain_rate);
                   }
 
-                // Step 3a: calculate viscoelastic (effective) viscosity
-                // Estimate the timestep size when in timestep 0
-                double dtc = this->get_timestep();
-                if (this->get_timestep_number() == 0 && this->get_timestep() == 0)
-                  dtc = std::min(std::min(this->get_parameters().maximum_time_step, this->get_parameters().maximum_first_time_step), elastic_rheology.elastic_timestep());
-                viscosity_pre_yield = dtc / elastic_rheology.elastic_timestep() * elastic_rheology.calculate_viscoelastic_viscosity(viscosity_pre_yield,
-                                      elastic_shear_modulus);
               }
 
             // Step 3b: calculate current (viscous or viscous + elastic) stress magnitude