[cig-commits] [commit] master: Update steinberger material model to run with the new material model interface (calling also viscosity for the creation of the adiabatic conditions). (fa5f128)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/aspect

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/a7135c1f7697d39efff2f47a79ca1e1395cff504...73a71ba37f203bfed63bb8b602fdbd30ab99b1af

>---------------------------------------------------------------

commit fa5f1281ab749971fbaa036adcb81e93092af6be
Author: Rene Gassmoeller <R.Gassmoeller at gmx.de>
Date:   Thu May 15 14:28:16 2014 -0500

    Update steinberger material model to run with the new material model interface (calling also viscosity for the creation of the adiabatic conditions).


>---------------------------------------------------------------

fa5f1281ab749971fbaa036adcb81e93092af6be
 source/material_model/steinberger.cc | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/source/material_model/steinberger.cc b/source/material_model/steinberger.cc
index 0a3ef2f..9c67e05 100644
--- a/source/material_model/steinberger.cc
+++ b/source/material_model/steinberger.cc
@@ -419,6 +419,7 @@ namespace aspect
                                   (new internal::MaterialLookup(datadirectory+material_file_names[i],interpolation)));
       lateral_viscosity_lookup.reset(new internal::LateralViscosityLookup(datadirectory+lateral_viscosity_file_name));
       radial_viscosity_lookup.reset(new internal::RadialViscosityLookup(datadirectory+radial_viscosity_file_name));
+      avg_temp.resize(1000);
     }
 
 
@@ -443,13 +444,22 @@ namespace aspect
                const Point<dim> &position) const
     {
       const double depth = this->get_geometry_model().depth(position);
-      const unsigned int idx = static_cast<unsigned int>(avg_temp.size() * depth / this->get_geometry_model().maximal_depth());
-      const double delta_temp = temperature-avg_temp[idx];
-      const double adia_temp = this->get_adiabatic_conditions().temperature(position);
 
-      const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temp/(temperature*adia_temp);
+      double vis_lateral;
+
+      if (&this->get_adiabatic_conditions())
+        {
+          const unsigned int idx = static_cast<unsigned int>(avg_temp.size() * depth / this->get_geometry_model().maximal_depth());
+          const double delta_temp = temperature-avg_temp[idx];
+          const double adia_temp = this->get_adiabatic_conditions().temperature(position);
+          const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temp/(temperature*adia_temp);
+          vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),1e2),1e-2);
+        }
+      else
+        // The adiabatic conditions are not there yet. This means we are currently creating them.
+        // Assume we are at an adiabatic temperature.
+        vis_lateral = 1;
 
-      const double vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),1e2),1e-2);
       const double vis_radial = radial_viscosity_lookup->radial_viscosity(depth);
 
       return std::max(std::min(vis_lateral * vis_radial,1e23),1e19);