Merge pull request #2170 from naliboff/fix_material_model_density

[WIP] Fix density used for thermal conductivity in visco_plastic and diffusion_dislocation models
geodynamics · Apr 27, 2018 · 89e76ee · 89e76ee
2 parents d99edf9 + 2cbd9be
commit 89e76ee
Show file tree

Hide file tree

Showing 5 changed files with 45 additions and 9 deletions.
diff --git a/source/initial_temperature/adiabatic.cc b/source/initial_temperature/adiabatic.cc
@@ -89,7 +89,14 @@ namespace aspect
       in.strain_rate.resize(0); // adiabat has strain=0.
       this->get_material_model().evaluate(in, out);
 
-      const double kappa = out.thermal_conductivities[0] / (out.densities[0] * out.specific_heat[0]);
+      const double kappa = ( (this->get_parameters().formulation_temperature_equation ==
+                              Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
+                             ?
+                             out.thermal_conductivities[0] /
+                             (this->get_adiabatic_conditions().density(in.position[0]) * out.specific_heat[0])
+                             :
+                             out.thermal_conductivities[0] / (out.densities[0] * out.specific_heat[0])
+                           );
 
       // analytical solution for the thermal boundary layer from half-space cooling model
       const double surface_cooling_temperature = age_top > 0.0 ?

diff --git a/source/material_model/diffusion_dislocation.cc b/source/material_model/diffusion_dislocation.cc
@@ -20,7 +20,7 @@
 
 #include <aspect/material_model/diffusion_dislocation.h>
 #include <aspect/utilities.h>
-
+#include <aspect/adiabatic_conditions/interface.h>
 
 namespace aspect
 {
@@ -305,8 +305,15 @@ namespace aspect
           out.thermal_expansion_coefficients[i] = thermal_expansivity;
           // Specific heat at the given positions.
           out.specific_heat[i] = heat_capacity;
-          // Thermal conductivity at the given positions.
-          out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity * density;
+          // Thermal conductivity at the given positions. If the temperature equation uses
+          // the reference density profile formulation, use the reference density to
+          // calculate thermal conductivity. Otherwise, use the real density.
+          if (this->get_parameters().formulation_temperature_equation ==
+              Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
+            out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity *
+                                            this->get_adiabatic_conditions().density(in.position[i]);
+          else
+            out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity * density;
           // Compressibility at the given positions.
           // The compressibility is given as
           // $\frac 1\rho \frac{\partial\rho}{\partial p}$.

diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
@@ -23,6 +23,7 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/base/signaling_nan.h>
 #include <aspect/newton.h>
+#include <aspect/adiabatic_conditions/interface.h>
 
 namespace aspect
 {
@@ -508,8 +509,15 @@ namespace aspect
           out.thermal_expansion_coefficients[i] = thermal_expansivity;
           // Specific heat at the given positions.
           out.specific_heat[i] = heat_capacity;
-          // Thermal conductivity at the given positions.
-          out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity * density;
+          // Thermal conductivity at the given positions. If the temperature equation uses
+          // the reference density profile formulation, use the reference density to
+          // calculate thermal conductivity. Otherwise, use the real density.
+          if (this->get_parameters().formulation_temperature_equation ==
+              Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
+            out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity *
+                                            this->get_adiabatic_conditions().density(in.position[i]);
+          else
+            out.thermal_conductivities[i] = thermal_diffusivity * heat_capacity * density;
           // Compressibility at the given positions.
           // The compressibility is given as
           // $\frac 1\rho \frac{\partial\rho}{\partial p}$.

diff --git a/source/postprocess/basic_statistics.cc b/source/postprocess/basic_statistics.cc
@@ -72,8 +72,14 @@ namespace aspect
 
           this->get_material_model().evaluate(in, out);
 
-          const double thermal_diffusivity = out.thermal_conductivities[0] /
-                                             (out.densities[0] * out.specific_heat[0]);
+          const double thermal_diffusivity = ( (this->get_parameters().formulation_temperature_equation ==
+                                                Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
+                                               ?
+                                               out.thermal_conductivities[0] /
+                                               (this->get_adiabatic_conditions().density(in.position[0]) * out.specific_heat[0])
+                                               :
+                                               out.thermal_conductivities[0] / (out.densities[0] * out.specific_heat[0])
+                                             );
 
           // check whether diffusivity is set to 0 (in case of backward advection)
           const double Ra = (thermal_diffusivity != 0) ?

diff --git a/source/postprocess/visualization/thermal_diffusivity.cc b/source/postprocess/visualization/thermal_diffusivity.cc
@@ -20,6 +20,7 @@
 
 
 #include <aspect/postprocess/visualization/thermal_diffusivity.h>
+#include <aspect/adiabatic_conditions/interface.h>
 
 
 
@@ -59,7 +60,14 @@ namespace aspect
         this->get_material_model().evaluate(in, out);
 
         for (unsigned int q=0; q<n_quadrature_points; ++q)
-          computed_quantities[q](0) = out.thermal_conductivities[q] / (out.densities[q] * out.specific_heat[q]);
+
+          if (this->get_parameters().formulation_temperature_equation ==
+              Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
+            computed_quantities[q](0) = out.thermal_conductivities[q] /
+                                        (this->get_adiabatic_conditions().density(in.position[q]) * out.specific_heat[q]);
+          else
+            computed_quantities[q](0) = out.thermal_conductivities[q] / (out.densities[q] * out.specific_heat[q]);
+
       }
     }
   }