[cig-commits] [commit] master: Changed steinberger material model to new interface. Made the incompressible formulation consistent and added a compressible formulation. (9b515a0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed May 21 14:02:53 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/a7135c1f7697d39efff2f47a79ca1e1395cff504...73a71ba37f203bfed63bb8b602fdbd30ab99b1af
>---------------------------------------------------------------
commit 9b515a0ca20635ec478d64a312379c84fd2f8afa
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date: Mon May 19 16:05:12 2014 -0500
Changed steinberger material model to new interface. Made the incompressible formulation consistent and added a compressible formulation.
>---------------------------------------------------------------
9b515a0ca20635ec478d64a312379c84fd2f8afa
include/aspect/material_model/steinberger.h | 27 +++-
source/material_model/steinberger.cc | 185 ++++++++++++++++++++--------
2 files changed, 159 insertions(+), 53 deletions(-)
diff --git a/include/aspect/material_model/steinberger.h b/include/aspect/material_model/steinberger.h
index 97310ea..1ef0847 100644
--- a/include/aspect/material_model/steinberger.h
+++ b/include/aspect/material_model/steinberger.h
@@ -49,7 +49,7 @@ namespace aspect
* @ingroup MaterialModels
*/
template <int dim>
- class Steinberger: public MaterialModel::InterfaceCompatibility<dim>, public ::aspect::SimulatorAccess<dim>
+ class Steinberger: public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
@@ -184,6 +184,16 @@ namespace aspect
*/
/**
+ * Function to compute the material properties in @p out given the
+ * inputs in @p in. If MaterialModelInputs.strain_rate has the length
+ * 0, then the viscosity does not need to be computed.
+ */
+ virtual
+ void
+ evaluate(const typename Interface<dim>::MaterialModelInputs &in,
+ typename Interface<dim>::MaterialModelOutputs &out) const;
+
+ /**
* @name Functions used in dealing with run-time parameters
* @{
*/
@@ -214,7 +224,20 @@ namespace aspect
unsigned int n_material_data;
std::string radial_viscosity_file_name;
std::string lateral_viscosity_file_name;
- virtual double get_deltat (const Point<dim> &position) const;
+ virtual double get_corrected_temperature (const double temperature,
+ const double pressure,
+ const Point<dim> &position) const;
+ virtual double get_corrected_pressure (const double temperature,
+ const double pressure,
+ const Point<dim> &position) const;
+ virtual double get_density (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const;
+ virtual double get_corrected_density (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const;
/**
* Pointer to an object that reads and processes data we get from
diff --git a/source/material_model/steinberger.cc b/source/material_model/steinberger.cc
index 2bdb507..32f2923 100644
--- a/source/material_model/steinberger.cc
+++ b/source/material_model/steinberger.cc
@@ -470,18 +470,58 @@ namespace aspect
template <int dim>
double
Steinberger<dim>::
- get_deltat (const Point<dim> &position) const
+ get_corrected_temperature (const double temperature,
+ const double pressure,
+ const Point<dim> &position) const
{
- if (!(&this->get_adiabatic_conditions()))
- return 0.0;
+ if (!(&this->get_adiabatic_conditions())
+ || this->include_adiabatic_heating()
+ || compressible)
+ return temperature;
+
+ return temperature
+ + this->get_adiabatic_conditions().temperature(position)
+ - this->get_adiabatic_surface_temperature();
+ }
+
- // We do not need a temperature correction in the compressible case
- if (compressible)
- return 0.0;
- static const bool a = this->include_adiabatic_heating();
- return a ? 0.0 : (this->get_adiabatic_conditions().temperature(position)
- - this->get_adiabatic_surface_temperature());
+ template <int dim>
+ double
+ Steinberger<dim>::
+ get_corrected_pressure (const double temperature,
+ const double pressure,
+ const Point<dim> &position) const
+ {
+ if (!(&this->get_adiabatic_conditions())
+ || compressible)
+ return pressure;
+
+ return this->get_adiabatic_conditions().pressure(position);
+ }
+
+ template <int dim>
+ double
+ Steinberger<dim>::
+ get_corrected_density (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const
+ {
+ const double rho = get_density(temperature,pressure,compositional_fields,position);
+
+ const double depth_rho = get_density(this->get_adiabatic_conditions().temperature(position),
+ pressure,
+ compositional_fields,
+ position);
+
+ const Point<dim> surface_point = this->get_geometry_model().representative_point(0.0);
+ const double surface_temperature = this->get_adiabatic_surface_temperature();
+ const double surface_pressure = this->get_surface_pressure();
+ const double surface_rho = get_density(surface_temperature,surface_pressure,compositional_fields,surface_point);
+
+ //Return the density scaled to an incompressible profile
+ return (rho / depth_rho) * surface_rho;
}
@@ -525,28 +565,25 @@ namespace aspect
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
double cp = 0.0;
if (!latent_heat)
{
if (n_material_data == 1)
- cp = material_lookup[0]->specific_heat(temperature+get_deltat(position),pressure);
+ cp = material_lookup[0]->specific_heat(temperature,pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- cp += compositional_fields[i] * material_lookup[i]->specific_heat(temperature+get_deltat(position),pressure);
+ cp += compositional_fields[i] * material_lookup[i]->specific_heat(temperature,pressure);
}
}
else
{
if (n_material_data == 1)
- cp = material_lookup[0]->dHdT(temperature+get_deltat(position),pressure);
+ cp = material_lookup[0]->dHdT(temperature,pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- cp += compositional_fields[i] * material_lookup[i]->dHdT(temperature+get_deltat(position),pressure);
+ cp += compositional_fields[i] * material_lookup[i]->dHdT(temperature,pressure);
cp = std::max(std::min(cp,6000.0),500.0);
}
}
@@ -571,32 +608,40 @@ namespace aspect
template <int dim>
double
Steinberger<dim>::
- density (const double temperature,
+ get_density (const double temperature,
const double pressure,
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
double rho = 0.0;
if (n_material_data == 1)
{
- // TODO: This is wrong for the incompressible case: The pressure in this case is not
- // the real pressure necessary to look the density up in a material table (due to the
- // missing depth dependency of density). Should use the adiabatic pressure instead, or
- // even better the adiabatic pressure + the deviation of pressure from the lateral pressure
- rho = material_lookup[0]->density(temperature+get_deltat(position),pressure);
+ rho = material_lookup[0]->density(temperature,pressure);
}
else
{
for (unsigned i = 0; i < n_material_data; i++)
- // TODO: Wrong. See above.
- rho += compositional_fields[i] * material_lookup[i]->density(temperature+get_deltat(position),pressure);
+ rho += compositional_fields[i] * material_lookup[i]->density(temperature,pressure);
}
+
return rho;
}
+ template <int dim>
+ double
+ Steinberger<dim>::
+ density (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const
+ {
+ if (compressible
+ || !(&this->get_adiabatic_conditions()))
+ return get_density(temperature,pressure,compositional_fields,position);
+ else
+ return get_corrected_density(temperature,pressure,compositional_fields,position);
+ }
+
template <int dim>
@@ -607,31 +652,26 @@ namespace aspect
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
double alpha = 0.0;
if (!latent_heat)
{
if (n_material_data == 1)
- // TODO: Wrong. See above
- alpha = material_lookup[0]->thermal_expansivity(temperature+get_deltat(position),pressure);
+ alpha = material_lookup[0]->thermal_expansivity(temperature,pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- // TODO: Wrong. See above
- alpha += compositional_fields[i] * material_lookup[i]->thermal_expansivity(temperature+get_deltat(position),pressure);
+ alpha += compositional_fields[i] * material_lookup[i]->thermal_expansivity(temperature,pressure);
}
}
else
{
double dHdp = 0.0;
if (n_material_data == 1)
- dHdp += material_lookup[0]->dHdp(temperature+get_deltat(position),pressure);
+ dHdp += material_lookup[0]->dHdp(temperature,pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- dHdp += compositional_fields[i] * material_lookup[i]->dHdp(temperature+get_deltat(position),pressure);
+ dHdp += compositional_fields[i] * material_lookup[i]->dHdp(temperature,pressure);
}
alpha = (1 - density(temperature,pressure,compositional_fields,position) * dHdp) / temperature;
alpha = std::max(std::min(alpha,1e-3),1e-5);
@@ -649,16 +689,18 @@ namespace aspect
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
+ //this function is not called from evaluate() so we need to care about
+ //corrections for temperature and pressure
+ const double corrected_temperature = get_corrected_temperature(temperature,pressure,position);
+ const double corrected_pressure = get_corrected_pressure(temperature,pressure,position);
+
double vp = 0.0;
if (n_material_data == 1)
- vp += material_lookup[0]->seismic_Vp(temperature+get_deltat(position),pressure);
+ vp += material_lookup[0]->seismic_Vp(corrected_temperature,corrected_pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- vp += compositional_fields[i] * material_lookup[i]->seismic_Vp(temperature+get_deltat(position),pressure);
+ vp += compositional_fields[i] * material_lookup[i]->seismic_Vp(corrected_temperature,corrected_pressure);
}
return vp;
}
@@ -673,16 +715,19 @@ namespace aspect
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
+ //this function is not called from evaluate() so we need to care about
+ //corrections for temperature and pressure
+ const double corrected_temperature = get_corrected_temperature(temperature,pressure,position);
+ const double corrected_pressure = get_corrected_pressure(temperature,pressure,position);
+
+
double vs = 0.0;
if (n_material_data == 1)
- vs += material_lookup[0]->seismic_Vs(temperature+get_deltat(position),pressure);
+ vs += material_lookup[0]->seismic_Vs(corrected_temperature,corrected_pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- vs += compositional_fields[i] * material_lookup[i]->seismic_Vs(temperature+get_deltat(position),pressure);
+ vs += compositional_fields[i] * material_lookup[i]->seismic_Vs(corrected_temperature,corrected_pressure);
}
return vs;
}
@@ -697,16 +742,16 @@ namespace aspect
const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- Assert ((n_material_data <= compositional_fields.size()) || (n_material_data == 1),
- ExcMessage("There are more material files provided than compositional"
- " Fields. This can not be intended."));
+ if (!compressible)
+ return 0.0;
+
double dRhodp = 0.0;
if (n_material_data == 1)
- dRhodp += material_lookup[0]->dRhodp(temperature+get_deltat(position),pressure);
+ dRhodp += material_lookup[0]->dRhodp(temperature,pressure);
else
{
for (unsigned i = 0; i < n_material_data; i++)
- dRhodp += compositional_fields[i] * material_lookup[i]->dRhodp(temperature+get_deltat(position),pressure);
+ dRhodp += compositional_fields[i] * material_lookup[i]->dRhodp(temperature,pressure);
}
const double rho = density(temperature,pressure,compositional_fields,position);
return (1/rho)*dRhodp;
@@ -794,6 +839,44 @@ namespace aspect
return compressible;
}
+ template <int dim>
+ void
+ Steinberger<dim>::evaluate(const typename Interface<dim>::MaterialModelInputs &in,
+ typename Interface<dim>::MaterialModelOutputs &out) const
+ {
+
+ Assert ((n_material_data <= in.composition[0].size()) || (n_material_data == 1),
+ ExcMessage("There are more material files provided than compositional"
+ " Fields. This can not be intended."));
+
+ for (unsigned int i=0; i < in.temperature.size(); ++i)
+ {
+ const double temperature = get_corrected_temperature(in.temperature[i],
+ in.pressure[i],
+ in.position[i]);
+ const double pressure = get_corrected_pressure(in.temperature[i],
+ in.pressure[i],
+ in.position[i]);
+
+ /* 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.
+ */
+ if (&this->get_adiabatic_conditions() && in.strain_rate.size())
+ out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], 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]);
+ out.thermal_conductivities[i] = thermal_conductivity (temperature, pressure, in.composition[i], in.position[i]);
+ out.compressibilities[i] = compressibility (temperature, pressure, in.composition[i], in.position[i]);
+ out.entropy_derivative_pressure[i] = 0;
+ out.entropy_derivative_temperature[i] = 0;
+ for (unsigned int c=0; c<in.composition[i].size(); ++c)
+ out.reaction_terms[i][c] = 0;
+ }
+ }
template <int dim>
More information about the CIG-COMMITS
mailing list