[cig-commits] [commit] master: add a new material model and cookbook with reactions between compositional fields (dc663a1)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri May 16 18:30:50 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/fec9fff95d1f00eabc29ab290d8294518b63b49a...7b658caa05d489ae3f8d61b8049b015e7be94f77
>---------------------------------------------------------------
commit dc663a1e6eaaa8b74746daa2c4391b170ac61d0b
Author: Juliane Dannberg <dannberg at gfz-potsdam.de>
Date: Fri May 16 20:35:27 2014 +0200
add a new material model and cookbook with reactions between compositional fields
>---------------------------------------------------------------
dc663a1e6eaaa8b74746daa2c4391b170ac61d0b
...elike-boundary.prm => composition-reaction.prm} | 28 +++-
.../{simple.h => composition_reaction.h} | 24 ++-
.../{simple.cc => composition_reaction.cc} | 165 ++++++++++++++-------
3 files changed, 154 insertions(+), 63 deletions(-)
diff --git a/cookbooks/platelike-boundary.prm b/cookbooks/composition-reaction.prm
similarity index 80%
copy from cookbooks/platelike-boundary.prm
copy to cookbooks/composition-reaction.prm
index c2e8684..7684ab7 100644
--- a/cookbooks/platelike-boundary.prm
+++ b/cookbooks/composition-reaction.prm
@@ -37,6 +37,11 @@ subsection Model settings
end
+subsection Compositional fields
+ set Number of fields = 2
+end
+
+
# We then set the temperature to one at the bottom and zero
# at the top:
subsection Boundary temperature model
@@ -48,6 +53,10 @@ subsection Boundary temperature model
end
end
+subsection Boundary composition model
+ set Model name = box
+end
+
# The velocity along the top boundary models a spreading
# center that is moving left and right:
@@ -78,14 +87,26 @@ subsection Initial conditions
end
end
+subsection Compositional initial conditions
+ set Model name = function
+
+ subsection Function
+ set Variable names = x,z
+ set Function expression = if(z<0.5,1,0);0
+ end
+end
+
subsection Material model
- set Model name = simple
+ set Model name = composition reaction
- subsection Simple model
+ subsection Composition reaction model
set Thermal conductivity = 1e-6
set Thermal expansion coefficient = 1e-4
set Viscosity = 1
+ set Density differential for compositional field 1 = -5
+ set Density differential for compositional field 2 = 5
+ set Reaction depth = 0.2
end
end
@@ -94,7 +115,7 @@ end
# mesh is refined and what to do with the solution once computed
subsection Mesh refinement
set Initial adaptive refinement = 0
- set Initial global refinement = 5
+ set Initial global refinement = 4
set Time steps between mesh refinement = 0
end
@@ -104,6 +125,7 @@ subsection Postprocess
subsection Visualization
set Time between graphical output = 0.1
+ set List of output variables = density, viscosity
end
end
diff --git a/include/aspect/material_model/simple.h b/include/aspect/material_model/composition_reaction.h
similarity index 88%
copy from include/aspect/material_model/simple.h
copy to include/aspect/material_model/composition_reaction.h
index 5432d10..13e7ffc 100644
--- a/include/aspect/material_model/simple.h
+++ b/include/aspect/material_model/composition_reaction.h
@@ -20,8 +20,8 @@
/* $Id$ */
-#ifndef __aspect__model_simple_h
-#define __aspect__model_simple_h
+#ifndef __aspect__model_composition_reaction_h
+#define __aspect__model_composition_reaction_h
#include <aspect/material_model/interface.h>
#include <aspect/simulator_access.h>
@@ -44,7 +44,7 @@ namespace aspect
* @ingroup MaterialModels
*/
template <int dim>
- class Simple : public MaterialModel::InterfaceCompatibility<dim>, public ::aspect::SimulatorAccess<dim>
+ class CompositionReaction : public MaterialModel::InterfaceCompatibility<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
/**
@@ -81,6 +81,12 @@ namespace aspect
const double pressure,
const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
+
+ virtual double reaction_term (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position,
+ const unsigned int compositional_variable) const;
/**
* @}
*/
@@ -185,7 +191,8 @@ namespace aspect
double reference_rho;
double reference_T;
double eta;
- double composition_viscosity_prefactor;
+ double composition_viscosity_prefactor_1;
+ double composition_viscosity_prefactor_2;
double thermal_viscosity_exponent;
double thermal_alpha;
double reference_specific_heat;
@@ -195,7 +202,14 @@ namespace aspect
*/
double k_value;
- double compositional_delta_rho;
+ double compositional_delta_rho_1;
+ double compositional_delta_rho_2;
+
+ /**
+ * Above this depth the compositional fields react:
+ * The first field gets converted to the second field.
+ */
+ double reaction_depth;
};
}
diff --git a/source/material_model/simple.cc b/source/material_model/composition_reaction.cc
similarity index 60%
copy from source/material_model/simple.cc
copy to source/material_model/composition_reaction.cc
index d535708..8ae408d 100644
--- a/source/material_model/simple.cc
+++ b/source/material_model/composition_reaction.cc
@@ -20,7 +20,7 @@
/* $Id$ */
-#include <aspect/material_model/simple.h>
+#include <aspect/material_model/composition_reaction.h>
#include <deal.II/base/parameter_handler.h>
using namespace dealii;
@@ -32,7 +32,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
viscosity (const double temperature,
const double,
const std::vector<double> &composition, /*composition*/
@@ -46,20 +46,25 @@ namespace aspect
temperature_dependence = 1.0;
double composition_dependence = 1.0;
- if ((composition_viscosity_prefactor != 1.0) && (composition.size() > 0))
- {
+ switch(composition.size())
+ {
+ case 0:
+ return composition_dependence * temperature_dependence * eta;
+ case 1:
//geometric interpolation
return (pow(10, ((1-composition[0]) * log10(eta*temperature_dependence)
- + composition[0] * log10(eta*composition_viscosity_prefactor*temperature_dependence))));
- }
-
- return composition_dependence * temperature_dependence * eta;
+ + composition[0] * log10(eta*composition_viscosity_prefactor_1*temperature_dependence))));
+ default:
+ return (pow(10, ((1 - 0.5*composition[0] - 0.5*composition[1]) * log10(eta*temperature_dependence)
+ + 0.5 * composition[0] * log10(eta*composition_viscosity_prefactor_1*temperature_dependence)
+ + 0.5 * composition[1] * log10(eta*composition_viscosity_prefactor_2*temperature_dependence))));
+ }
}
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
reference_viscosity () const
{
return eta;
@@ -67,7 +72,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
reference_density () const
{
return reference_rho;
@@ -75,7 +80,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
reference_thermal_expansion_coefficient () const
{
return thermal_alpha;
@@ -83,7 +88,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
specific_heat (const double,
const double,
const std::vector<double> &, /*composition*/
@@ -94,7 +99,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
reference_cp () const
{
return reference_specific_heat;
@@ -102,7 +107,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
thermal_conductivity (const double,
const double,
const std::vector<double> &, /*composition*/
@@ -113,7 +118,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
reference_thermal_diffusivity () const
{
return k_value/(reference_rho*reference_specific_heat);
@@ -121,24 +126,28 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
density (const double temperature,
const double,
const std::vector<double> &compositional_fields, /*composition*/
const Point<dim> &) const
{
- const double c = compositional_fields.size()>0?
- std::max(0.0, compositional_fields[0])
- :
- 0.0;
+ const double c1 = compositional_fields.size()>0?
+ std::max(0.0, compositional_fields[0])
+ :
+ 0.0;
+ const double c2 = compositional_fields.size()>1?
+ std::max(0.0, compositional_fields[1])
+ :
+ 0.0;
return reference_rho * (1 - thermal_alpha * (temperature - reference_T))
- + compositional_delta_rho * c;
+ + compositional_delta_rho_1 * c1 + compositional_delta_rho_2 * c2;
}
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
const std::vector<double> &, /*composition*/
@@ -150,7 +159,7 @@ namespace aspect
template <int dim>
double
- Simple<dim>::
+ CompositionReaction<dim>::
compressibility (const double,
const double,
const std::vector<double> &, /*composition*/
@@ -159,9 +168,37 @@ namespace aspect
return 0.0;
}
+
+ template <int dim>
+ double
+ CompositionReaction<dim>::
+ reaction_term (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position,
+ const unsigned int compositional_variable) const
+ {
+ const double depth = this->get_geometry_model().depth(position);
+ double delta_C = 0.0;
+ switch(compositional_variable)
+ {
+ case 0:
+ if(depth < reaction_depth) delta_C = -compositional_fields[0];
+ break;
+ case 1:
+ if(depth < reaction_depth) delta_C = compositional_fields[0];
+ break;
+ default:
+ delta_C = 0.0;
+ break;
+ }
+ return delta_C;
+ }
+
+
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
viscosity_depends_on (const NonlinearDependence::Dependence dependence) const
{
// compare this with the implementation of the viscosity() function
@@ -172,7 +209,7 @@ namespace aspect
return true;
else if (((dependence & NonlinearDependence::compositional_fields) != NonlinearDependence::none)
&&
- (composition_viscosity_prefactor != 1.0))
+ (composition_viscosity_prefactor_1 != 1.0 || composition_viscosity_prefactor_2 != 1.0))
return true;
else
return false;
@@ -181,7 +218,7 @@ namespace aspect
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
density_depends_on (const NonlinearDependence::Dependence dependence) const
{
// compare this with the implementation of the density() function
@@ -192,7 +229,7 @@ namespace aspect
return true;
else if (((dependence & NonlinearDependence::compositional_fields) != NonlinearDependence::none)
&&
- (compositional_delta_rho != 0))
+ (compositional_delta_rho_1 != 0 || compositional_delta_rho_2 != 0))
return true;
else
return false;
@@ -200,7 +237,7 @@ namespace aspect
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
compressibility_depends_on (const NonlinearDependence::Dependence) const
{
return false;
@@ -208,7 +245,7 @@ namespace aspect
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
specific_heat_depends_on (const NonlinearDependence::Dependence) const
{
return false;
@@ -216,7 +253,7 @@ namespace aspect
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
thermal_conductivity_depends_on (const NonlinearDependence::Dependence dependence) const
{
return false;
@@ -225,7 +262,7 @@ namespace aspect
template <int dim>
bool
- Simple<dim>::
+ CompositionReaction<dim>::
is_compressible () const
{
return false;
@@ -235,11 +272,11 @@ namespace aspect
template <int dim>
void
- Simple<dim>::declare_parameters (ParameterHandler &prm)
+ CompositionReaction<dim>::declare_parameters (ParameterHandler &prm)
{
prm.enter_subsection("Material model");
{
- prm.enter_subsection("Simple model");
+ prm.enter_subsection("Composition reaction model");
{
prm.declare_entry ("Reference density", "3300",
Patterns::Double (0),
@@ -250,11 +287,16 @@ namespace aspect
prm.declare_entry ("Viscosity", "5e24",
Patterns::Double (0),
"The value of the constant viscosity. Units: $kg/m/s$.");
- prm.declare_entry ("Composition viscosity prefactor", "1.0",
+ prm.declare_entry ("Composition viscosity prefactor 1", "1.0",
Patterns::Double (0),
"A linear dependency of viscosity on the first compositional field. "
"Dimensionless prefactor. With a value of 1.0 (the default) the "
"viscosity does not depend on the composition.");
+ prm.declare_entry ("Composition viscosity prefactor 2", "1.0",
+ Patterns::Double (0),
+ "A linear dependency of viscosity on the second compositional field. "
+ "Dimensionless prefactor. With a value of 1.0 (the default) the "
+ "viscosity does not depend on the composition.");
prm.declare_entry ("Thermal viscosity exponent", "0.0",
Patterns::Double (0),
"The temperature dependence of viscosity. Dimensionless exponent.");
@@ -277,10 +319,26 @@ namespace aspect
"model, we make the following assumptions: if no compositional fields "
"are used in the current simulation, then the density is simply the usual "
"one with its linear dependence on the temperature. If there are compositional "
- "fields, then the density only depends on the first one in such a way that "
+ "fields, then the density only depends on the first and the second one in such a way that "
+ "the density has an additional term of the kind $+\\Delta \\rho \\; c_1(\\mathbf x)$. "
+ "This parameter describes the value of $\\Delta \\rho$ for the first field. "
+ "Units: $kg/m^3/\\textrm{unit change in composition}$.");
+ prm.declare_entry ("Density differential for compositional field 2", "0",
+ Patterns::Double(),
+ "If compositional fields are used, then one would frequently want "
+ "to make the density depend on these fields. In this simple material "
+ "model, we make the following assumptions: if no compositional fields "
+ "are used in the current simulation, then the density is simply the usual "
+ "one with its linear dependence on the temperature. If there are compositional "
+ "fields, then the density only depends on the first and the second one in such a way that "
"the density has an additional term of the kind $+\\Delta \\rho \\; c_1(\\mathbf x)$. "
- "This parameter describes the value of $\\Delta \\rho$. Units: $kg/m^3/\\textrm{unit "
- "change in composition}$.");
+ "This parameter describes the value of $\\Delta \\rho$ for the second field. "
+ "Units: $kg/m^3/\\textrm{unit change in composition}$.");
+ prm.declare_entry ("Reaction depth", "0",
+ Patterns::Double (0),
+ "Above this depth the compositional fields react: "
+ "The first field gets converted to the second field. "
+ "Units: $m$.");
}
prm.leave_subsection();
}
@@ -291,24 +349,27 @@ namespace aspect
template <int dim>
void
- Simple<dim>::parse_parameters (ParameterHandler &prm)
+ CompositionReaction<dim>::parse_parameters (ParameterHandler &prm)
{
prm.enter_subsection("Material model");
{
- prm.enter_subsection("Simple model");
+ prm.enter_subsection("Composition reaction model");
{
reference_rho = prm.get_double ("Reference density");
reference_T = prm.get_double ("Reference temperature");
eta = prm.get_double ("Viscosity");
- composition_viscosity_prefactor = prm.get_double ("Composition viscosity prefactor");
+ composition_viscosity_prefactor_1 = prm.get_double ("Composition viscosity prefactor 1");
+ composition_viscosity_prefactor_2 = prm.get_double ("Composition viscosity prefactor 2");
thermal_viscosity_exponent = prm.get_double ("Thermal viscosity exponent");
k_value = prm.get_double ("Thermal conductivity");
reference_specific_heat = prm.get_double ("Reference specific heat");
thermal_alpha = prm.get_double ("Thermal expansion coefficient");
- compositional_delta_rho = prm.get_double ("Density differential for compositional field 1");
+ compositional_delta_rho_1 = prm.get_double ("Density differential for compositional field 1");
+ compositional_delta_rho_2 = prm.get_double ("Density differential for compositional field 2");
+ reaction_depth = prm.get_double ("Reaction depth");
if (thermal_viscosity_exponent!=0.0 && reference_T == 0.0)
- AssertThrow(false, ExcMessage("Error: Material model simple with Thermal viscosity exponent can not have reference_T=0."));
+ AssertThrow(false, ExcMessage("Error: Material model composition reaction with Thermal viscosity exponent can not have reference_T=0."));
}
prm.leave_subsection();
}
@@ -322,18 +383,12 @@ namespace aspect
{
namespace MaterialModel
{
- ASPECT_REGISTER_MATERIAL_MODEL(Simple,
- "simple",
- "A simple material model that has constant values "
- "for all coefficients but the density and viscosity. "
- "This model uses the formulation that assumes an incompressible"
- " medium despite the fact that the density follows the law "
- "$\\rho(T)=\\rho_0(1-\\beta(T-T_{\\text{ref}})$. "
- "The temperature dependency of viscosity is "
- " switched off by default and follows the formula "
- "$\\eta(T)=\\eta_0 e^{\\eta_T \\cdot \\Delta T / T_{\\text{ref}})}$."
- "The value for the components of this formula and additional "
- "parameters are read from the parameter file in subsection "
- "'Simple model'.")
+ ASPECT_REGISTER_MATERIAL_MODEL(CompositionReaction,
+ "composition reaction",
+ "A material model that behaves in the same way as "
+ "the simple material model, but includes two compositional "
+ "fields and a reaction between them. Above a depth given "
+ "in the input file, the first fields gets converted to the "
+ "second field. ")
}
}
More information about the CIG-COMMITS
mailing list