[cig-commits] [commit] master: Made steinberger material viscosity calculation more flexible. (79c2eb6)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Jun 20 00:23:14 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/9398730971b5217314dc121bc16852d75951d0b4...646b028803359614a76a003795e6d97f0a5a0603
>---------------------------------------------------------------
commit 79c2eb602b1893283df064bc5b26140ae7fcfc29
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date: Thu Jun 19 17:41:40 2014 +0200
Made steinberger material viscosity calculation more flexible.
>---------------------------------------------------------------
79c2eb602b1893283df064bc5b26140ae7fcfc29
include/aspect/material_model/steinberger.h | 4 ++
source/material_model/steinberger.cc | 70 +++++++++++++++++++++++------
2 files changed, 60 insertions(+), 14 deletions(-)
diff --git a/include/aspect/material_model/steinberger.h b/include/aspect/material_model/steinberger.h
index 39c32f6..51a6394 100644
--- a/include/aspect/material_model/steinberger.h
+++ b/include/aspect/material_model/steinberger.h
@@ -220,8 +220,12 @@ namespace aspect
bool interpolation;
bool latent_heat;
bool compressible;
+ bool use_lateral_average_temperature;
double reference_eta;
std::vector<double> avg_temp;
+ double min_eta;
+ double max_eta;
+ double max_lateral_eta_variation;
std::string datadirectory;
std::vector<std::string> material_file_names;
unsigned int n_material_data;
diff --git a/source/material_model/steinberger.cc b/source/material_model/steinberger.cc
index 25c5411..a0fda3c 100644
--- a/source/material_model/steinberger.cc
+++ b/source/material_model/steinberger.cc
@@ -433,7 +433,8 @@ namespace aspect
Steinberger<dim>::
update()
{
- this->get_depth_average_temperature(avg_temp);
+ if (use_lateral_average_temperature)
+ this->get_depth_average_temperature(avg_temp);
}
@@ -448,16 +449,25 @@ namespace aspect
const Point<dim> &position) const
{
const double depth = this->get_geometry_model().depth(position);
+ const double adiabatic_temperature = this->get_adiabatic_conditions().temperature(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);
- const double vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),1e2),1e-2);
+ double delta_temperature;
+ if (use_lateral_average_temperature)
+ {
+ const unsigned int idx = static_cast<unsigned int>(avg_temp.size() * depth / this->get_geometry_model().maximal_depth());
+ delta_temperature = temperature-avg_temp[idx];
+ }
+ else
+ const double delta_temperature = temperature-adiabatic_temperature;
+
+ // For an explanation on this formula see the Steinberger & Calderwood 2006 paper
+ const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
+ // Limit the lateral viscosity variation to a reasonable interval
+ const double vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),max_lateral_eta_variation),1/max_lateral_eta_variation);
const double vis_radial = radial_viscosity_lookup->radial_viscosity(depth);
- return std::max(std::min(vis_lateral * vis_radial,1e23),1e19);
+ return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
}
@@ -860,11 +870,20 @@ namespace aspect
/* We are only asked to give viscosities if strain_rate.size() > 0
* and we can only calculate it if adiabatic_conditions are available.
* Note that the used viscosity formulation needs the not
- * corrected temperatures.
+ * corrected temperatures in case we compare it to the lateral
+ * temperature average.
*/
if (this->get_adiabatic_conditions().is_initialized() && in.strain_rate.size())
- out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
-
+ {
+ if (use_lateral_average_temperature)
+ {
+ out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
+ }
+ else
+ {
+ out.viscosities[i] = viscosity (temperature, pressure, in.composition[i], in.strain_rate[i], in.position[i]);
+ }
+ }
out.densities[i] = density (temperature, pressure, in.composition[i], in.position[i]);
out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (temperature, pressure, in.composition[i], in.position[i]);
out.specific_heat[i] = specific_heat (temperature, pressure, in.composition[i], in.position[i]);
@@ -905,6 +924,13 @@ namespace aspect
prm.declare_entry ("Lateral viscosity file name", "temp-viscosity-prefactor.txt",
Patterns::Anything (),
"The file name of the lateral viscosity data. ");
+ prm.declare_entry ("Use lateral average temperature for viscosity", "true",
+ Patterns::Bool (),
+ "Whether to use to use the laterally averaged temperature "
+ "instead of the adiabatic temperature for the viscosity "
+ "calculation. This ensures that the laterally averaged "
+ "viscosities remain more or less constant over the model "
+ "runtime. This behaviour might or might not be desired.");
prm.declare_entry ("Bilinear interpolation", "true",
Patterns::Bool (),
"Whether to use bilinear interpolation to compute "
@@ -921,6 +947,19 @@ namespace aspect
prm.declare_entry ("Reference viscosity", "1e23",
Patterns::Double(0),
"The reference viscosity that is used for pressure scaling. ");
+ prm.declare_entry ("Minimum viscosity", "1e19",
+ Patterns::Double(0),
+ "The minimum viscosity that is allowed in the viscosity "
+ "calculation. Smaller values will be cut off.");
+ prm.declare_entry ("Maximum viscosity", "1e23",
+ Patterns::Double(0),
+ "The maximum viscosity that is allowed in the viscosity "
+ "calculation. Larger values will be cut off.");
+ prm.declare_entry ("Maximum lateral viscosity variation", "1e2",
+ Patterns::Double(0),
+ "The relative cutoff value for lateral viscosity variations "
+ "caused by temperature deviations. The viscosity may vary "
+ "laterally by this factor squared.");
prm.leave_subsection();
}
prm.leave_subsection();
@@ -949,12 +988,15 @@ namespace aspect
material_file_names = Utilities::split_string_list
(prm.get ("Material file names"));
radial_viscosity_file_name = prm.get ("Radial viscosity file name");
- lateral_viscosity_file_name = prm.get ("Lateral viscosity file name");
+ lateral_viscosity_file_name = prm.get ("Lateral viscosity file name");
+ use_lateral_average_temperature = prm.get_bool ("Use lateral average temperature for viscosity");
interpolation = prm.get_bool ("Bilinear interpolation");
latent_heat = prm.get_bool ("Latent heat");
- compressible = prm.get_bool ("Compressible");
- reference_eta = prm.get_double ("Reference viscosity");
-
+ compressible = prm.get_bool ("Compressible");
+ reference_eta = prm.get_double ("Reference viscosity");
+ min_eta = prm.get_double ("Minimum viscosity");
+ max_eta = prm.get_double ("Maximum viscosity");
+ max_lateral_eta_variation = prm.get_double ("Maximum lateral viscosity variation");
prm.leave_subsection();
}
More information about the CIG-COMMITS
mailing list