[cig-commits] commit 2053 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Mon Nov 25 06:37:01 PST 2013
Revision 2053
The last piece of Juliane's r1679 from her branch: Implement latent heat and boundary values for compositional fields.
_U trunk/aspect/
U trunk/aspect/doc/modules/changes.h
A trunk/aspect/include/aspect/boundary_composition/
D trunk/aspect/include/aspect/boundary_composition/box.h
A trunk/aspect/include/aspect/boundary_composition/box.h
D trunk/aspect/include/aspect/boundary_composition/interface.h
A trunk/aspect/include/aspect/boundary_composition/interface.h
U trunk/aspect/include/aspect/material_model/interface.h
U trunk/aspect/include/aspect/simulator.h
A trunk/aspect/source/boundary_composition/
D trunk/aspect/source/boundary_composition/box.cc
A trunk/aspect/source/boundary_composition/box.cc
D trunk/aspect/source/boundary_composition/interface.cc
A trunk/aspect/source/boundary_composition/interface.cc
U trunk/aspect/source/material_model/interface.cc
U trunk/aspect/source/simulator/assembly.cc
U trunk/aspect/source/simulator/core.cc
U trunk/aspect/source/simulator/parameters.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2053&peg=2053
Diff:
Index: trunk/aspect
===================================================================
--- trunk/aspect 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect 2013-11-25 14:36:21 UTC (rev 2053)
Property changes on: trunk/aspect
___________________________________________________________________
Modified: svn:mergeinfo
## -1,3 +1,4 ##
/branches/active_compositions:1257-1296
/branches/compositional:1141-1251
/branches/fully-nonlinear:542-728
+/branches/j-dannberg:1679
\ No newline at end of property
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/doc/modules/changes.h 2013-11-25 14:36:21 UTC (rev 2053)
@@ -8,6 +8,20 @@
</p>
<ol>
+ <li>New: One can now select in the input file that the model should
+ include latent heat. The generation of latent heat then needs to be
+ described in the material model.
+ <br>
+ (Juliane Dannberg 2013/11/24)
+
+ <li>New: It is now possible to prescribe boundary values for
+ compositional fields in cases where there is inflow through
+ a segment of the boundary. This is implemented through a set
+ of plugins for boundary values in the same way as is done for
+ temperature boundary values.
+ <br>
+ (Juliane Dannberg 2013/11/24)
+
<li>New: There is now a refinement criterion "topography" that makes sure
the mesh is always refined at the surface of the domain.
<br>
Copied: trunk/aspect/include/aspect/boundary_composition/box.h (from rev 1679, branches/j-dannberg/include/aspect/boundary_composition/box.h)
===================================================================
--- trunk/aspect/include/aspect/boundary_composition/box.h (rev 0)
+++ trunk/aspect/include/aspect/boundary_composition/box.h 2013-11-25 14:36:21 UTC (rev 2053)
@@ -0,0 +1,84 @@
+/*
+ Copyright (C) 2013 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id: box.h 1538 2013-01-06 03:12:23Z bangerth $ */
+
+
+#ifndef __aspect__boundary_composition_box_h
+#define __aspect__boundary_composition_box_h
+
+#include <aspect/boundary_composition/interface.h>
+
+
+namespace aspect
+{
+ namespace BoundaryComposition
+ {
+ /**
+ * A class that implements a composition boundary condition for a box
+ * geometry.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ template <int dim>
+ class Box : public Interface<dim>
+ {
+ public:
+ /**
+ * Return the composition that is to hold at a particular location on the
+ * boundary of the domain. This function returns constant compositions
+ * at the left and right boundaries.
+ *
+ * @param geometry_model The geometry model that describes the domain. This may
+ * be used to determine whether the boundary composition model is implemented
+ * for this geometry.
+ * @param boundary_indicator The boundary indicator of the part of the boundary
+ * of the domain on which the point is located at which we are requesting the
+ * composition.
+ * @param location The location of the point at which we ask for the composition.
+ **/
+ virtual
+ double composition (const GeometryModel::Interface<dim> &geometry_model,
+ const unsigned int boundary_indicator,
+ const Point<dim> &location) const;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ * This class declares the inner and outer boundary compositions.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ private:
+ double composition_[2*dim];
+ };
+ }
+}
+
+
+#endif
Copied: trunk/aspect/include/aspect/boundary_composition/interface.h (from rev 1679, branches/j-dannberg/include/aspect/boundary_composition/interface.h)
===================================================================
--- trunk/aspect/include/aspect/boundary_composition/interface.h (rev 0)
+++ trunk/aspect/include/aspect/boundary_composition/interface.h 2013-11-25 14:36:21 UTC (rev 2053)
@@ -0,0 +1,163 @@
+/*
+ Copyright (C) 2013 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id: interface.h 1538 2013-01-06 03:12:23Z bangerth $ */
+
+
+#ifndef __aspect__boundary_composition_interface_h
+#define __aspect__boundary_composition_interface_h
+
+#include <aspect/plugins.h>
+#include <aspect/geometry_model/interface.h>
+
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/distributed/tria.h>
+
+
+namespace aspect
+{
+ /**
+ * A namespace for the definition of things that have to do with describing
+ * the boundary values for the composition.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ namespace BoundaryComposition
+ {
+ using namespace dealii;
+
+ /**
+ * Base class for classes that describe composition boundary values.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ template <int dim>
+ class Interface
+ {
+ public:
+ /**
+ * Destructor. Made virtual to enforce that derived classes also have
+ * virtual destructors.
+ */
+ virtual ~Interface();
+
+ /**
+ * Return the composition that is to hold at a particular location on the
+ * boundary of the domain.
+ *
+ * @param geometry_model The geometry model that describes the domain. This may
+ * be used to determine whether the boundary composition model is implemented
+ * for this geometry.
+ * @param boundary_indicator The boundary indicator of the part of the boundary
+ * of the domain on which the point is located at which we are requesting the
+ * composition.
+ * @param location The location of the point at which we ask for the composition.
+ **/
+ virtual
+ double composition (const GeometryModel::Interface<dim> &geometry_model,
+ const unsigned int boundary_indicator,
+ const Point<dim> &location) const = 0;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ * The default implementation of this function does not describe
+ * any parameters. Consequently, derived classes do not have to
+ * overload this function if they do not take any runtime parameters.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file. The default implementation of this function does not read
+ * any parameters. Consequently, derived classes do not have to
+ * overload this function if they do not take any runtime parameters.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+ };
+
+
+ /**
+ * Register a boundary composition model so that it can be selected from the parameter file.
+ *
+ * @param name A string that identifies the boundary composition model
+ * @param description A text description of what this model
+ * does and that will be listed in the documentation of
+ * the parameter file.
+ * @param declare_parameters_function A pointer to a function that can be used to
+ * declare the parameters that this geometry model wants to read from input files.
+ * @param factory_function A pointer to a function that can create an object of
+ * this boundary composition model.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ template <int dim>
+ void
+ register_boundary_composition (const std::string &name,
+ const std::string &description,
+ void (*declare_parameters_function) (ParameterHandler &),
+ Interface<dim> *(*factory_function) ());
+
+ /**
+ * A function that given the name of a model returns a pointer to an object
+ * that describes it. Ownership of the pointer is transferred to the caller.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ template <int dim>
+ Interface<dim> *
+ create_boundary_composition (ParameterHandler &prm);
+
+
+ /**
+ * Declare the runtime parameters of the registered boundary composition models.
+ *
+ * @ingroup BoundaryCompositions
+ */
+ template <int dim>
+ void
+ declare_parameters (ParameterHandler &prm);
+
+
+ /**
+ * Given a class name, a name, and a description for the parameter file for a boundary composition model, register it with
+ * the functions that can declare their parameters and create these objects.
+ *
+ * @ingroup BoundaryCompositions
+ */
+#define ASPECT_REGISTER_BOUNDARY_COMPOSITION_MODEL(classname, name, description) \
+ template class classname<2>; \
+ template class classname<3>; \
+ namespace ASPECT_REGISTER_BOUNDARY_COMPOSITION_MODEL_ ## classname \
+ { \
+ aspect::internal::Plugins::RegisterHelper<Interface<2>,classname<2> > \
+ dummy_ ## classname ## _2d (&aspect::BoundaryComposition::register_boundary_composition<2>, \
+ name, description); \
+ aspect::internal::Plugins::RegisterHelper<Interface<3>,classname<3> > \
+ dummy_ ## classname ## _3d (&aspect::BoundaryComposition::register_boundary_composition<3>, \
+ name, description); \
+ }
+ }
+}
+
+
+#endif
Modified: trunk/aspect/include/aspect/material_model/interface.h
===================================================================
--- trunk/aspect/include/aspect/material_model/interface.h 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/include/aspect/material_model/interface.h 2013-11-25 14:36:21 UTC (rev 2053)
@@ -564,6 +564,14 @@
* $ rac 1
ho rac{\partial
ho}{\partial p}$.
*/
std::vector<double> compressibilities;
+ /**
+ * Pressure derivative of entropy at the given positions.
+ */
+ std::vector<double> entropy_derivative_pressure;
+ /**
+ * Temperature derivative of entropy at the given positions.
+ */
+ std::vector<double> entropy_derivative_temperature;
};
/**
@@ -681,7 +689,36 @@
const std::vector<double> &compositional_fields,
const Point<dim> &position) const=0;
- /**
+ /**
+ * Return the product of the change in entropy across phase
+ * transitions, the pressure derivative of the phase function
+ * (if this is the pressure derivative) or the product of the
+ * former two and the Clapeyron slope (if this is the
+ * temperature derivative).
+ * The entropy change across a phase transition can be calculated
+ * as $ rac{\gamma \Delta
ho}{
ho_ ext{light}
ho_ ext{heavy}}$.
+ * $\gamma$ is the Clapeyron slope of the phase transition,
+ * $\Delta
ho$ is the density jump across the phase transition,
+ * $
ho_ ext{light}$ is the density of the light material
+ * (above the phase transition) and $
ho_ ext{heavy}$ the
+ * density of the heavy material (below the phase transition).
+ * The phase function hat values ranging from 0 to 1 indicating
+ * which percentage of the material has already undergone the phase
+ * transition. Its argument is usually the excess pressure
+ * $\pi = p - p_0 - \gamma T$, where $p_0$ is the zero-degree
+ * transition pressure.
+ *
+ * This function has a default implementation that computes sets
+ * the entropy gradient to zero (assuming no phase changes).
+ */
+ virtual double entropy_derivative (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position,
+ const NonlinearDependence::Dependence dependence) const;
+
+
+ /**
* Return the thermal conductivity $k$ of the model as a function of temperature,
* pressure and position. The units of $k$ are $ extrm{W} / extrm{m} / extrm{K}$
* in 3d, and $ extrm{W} / extrm{K}$ in 2d. This is easily see by considering that
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/include/aspect/simulator.h 2013-11-25 14:36:21 UTC (rev 2053)
@@ -45,6 +45,7 @@
#include <aspect/geometry_model/interface.h>
#include <aspect/gravity_model/interface.h>
#include <aspect/boundary_temperature/interface.h>
+#include <aspect/boundary_composition/interface.h>
#include <aspect/initial_conditions/interface.h>
#include <aspect/compositional_initial_conditions/interface.h>
#include <aspect/velocity_boundary_conditions/interface.h>
@@ -183,8 +184,10 @@
*/
bool include_shear_heating;
bool include_adiabatic_heating;
+ bool include_latent_heat;
double radiogenic_heating_rate;
std::set<types::boundary_id> fixed_temperature_boundary_indicators;
+ std::set<types::boundary_id> fixed_composition_boundary_indicators;
std::set<types::boundary_id> zero_velocity_boundary_indicators;
std::set<types::boundary_id> tangential_velocity_boundary_indicators;
/**
@@ -1063,6 +1066,7 @@
const std::auto_ptr<MaterialModel::Interface<dim> > material_model;
const std::auto_ptr<GravityModel::Interface<dim> > gravity_model;
const std::auto_ptr<BoundaryTemperature::Interface<dim> > boundary_temperature;
+ const std::auto_ptr<BoundaryComposition::Interface<dim> > boundary_composition;
std::auto_ptr<CompositionalInitialConditions::Interface<dim> > compositional_initial_conditions;
std::auto_ptr<const AdiabaticConditions<dim> > adiabatic_conditions;
std::auto_ptr<InitialConditions::Interface<dim> > initial_conditions;
Copied: trunk/aspect/source/boundary_composition/box.cc (from rev 1679, branches/j-dannberg/source/boundary_composition/box.cc)
===================================================================
--- trunk/aspect/source/boundary_composition/box.cc (rev 0)
+++ trunk/aspect/source/boundary_composition/box.cc 2013-11-25 14:36:21 UTC (rev 2053)
@@ -0,0 +1,140 @@
+/*
+ Copyright (C) 2013 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id: box.cc 1545 2013-01-21 16:14:41Z heister $ */
+
+
+#include <aspect/boundary_composition/box.h>
+#include <aspect/geometry_model/box.h>
+
+#include <utility>
+#include <limits>
+
+
+namespace aspect
+{
+ namespace BoundaryComposition
+ {
+// ------------------------------ Box -------------------
+
+ template <int dim>
+ double
+ Box<dim>::
+ composition (const GeometryModel::Interface<dim> &geometry_model,
+ const unsigned int boundary_indicator,
+ const Point<dim> &location) const
+ {
+ // verify that the geometry is in fact a box since only
+ // for this geometry do we know for sure what boundary indicators it
+ // uses and what they mean
+ Assert (dynamic_cast<const GeometryModel::Box<dim>*>(&geometry_model)
+ != 0,
+ ExcMessage ("This boundary model is only implemented if the geometry is "
+ "in fact a box."));
+
+ Assert (boundary_indicator<2*dim, ExcMessage ("Unknown boundary indicator."));
+ return composition_[boundary_indicator];
+ }
+
+ template <int dim>
+ void
+ Box<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Boundary composition model");
+ {
+ prm.enter_subsection("Box");
+ {
+ prm.declare_entry ("Left composition", "0",
+ Patterns::Double (),
+ "Composition at the left boundary (at minimal x-value). Units: K.");
+ prm.declare_entry ("Right composition", "0",
+ Patterns::Double (),
+ "Composition at the right boundary (at maximal x-value). Units: K.");
+ prm.declare_entry ("Bottom composition", "1",
+ Patterns::Double (),
+ "Composition at the bottom boundary (at minimal z-value). Units: K.");
+ prm.declare_entry ("Top composition", "0",
+ Patterns::Double (),
+ "Composition at the top boundary (at maximal x-value). Units: K.");
+ if (dim==3)
+ {
+ prm.declare_entry ("Front composition", "0",
+ Patterns::Double (),
+ "Composition at the front boundary (at minimal y-value). Units: K.");
+ prm.declare_entry ("Back composition", "0",
+ Patterns::Double (),
+ "Composition at the back boundary (at maximal y-value). Units: K.");
+ }
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+
+
+ template <int dim>
+ void
+ Box<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Boundary composition model");
+ {
+ prm.enter_subsection("Box");
+ {
+ switch (dim)
+ {
+ case 2:
+ composition_[0] = prm.get_double ("Left composition");
+ composition_[1] = prm.get_double ("Right composition");
+ composition_[2] = prm.get_double ("Bottom composition");
+ composition_[3] = prm.get_double ("Top composition");
+ break;
+
+ case 3:
+ composition_[0] = prm.get_double ("Left composition");
+ composition_[1] = prm.get_double ("Right composition");
+ composition_[2] = prm.get_double ("Front composition");
+ composition_[3] = prm.get_double ("Back composition");
+ composition_[4] = prm.get_double ("Bottom composition");
+ composition_[5] = prm.get_double ("Top composition");
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+
+
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace BoundaryComposition
+ {
+ ASPECT_REGISTER_BOUNDARY_COMPOSITION_MODEL(Box,
+ "box",
+ "A model in which the composition is chosen constant on "
+ "all the sides of a box.")
+ }
+}
Copied: trunk/aspect/source/boundary_composition/interface.cc (from rev 1679, branches/j-dannberg/source/boundary_composition/interface.cc)
===================================================================
--- trunk/aspect/source/boundary_composition/interface.cc (rev 0)
+++ trunk/aspect/source/boundary_composition/interface.cc 2013-11-25 14:36:21 UTC (rev 2053)
@@ -0,0 +1,162 @@
+/*
+ Copyright (C) 2013 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id: interface.cc 885 2012-04-03 13:35:03Z bangerth $ */
+
+
+#include <aspect/global.h>
+#include <aspect/boundary_composition/interface.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx1x/tuple.h>
+
+#include <list>
+
+
+namespace aspect
+{
+ namespace BoundaryComposition
+ {
+ template <int dim>
+ Interface<dim>::~Interface ()
+ {}
+
+
+ template <int dim>
+
+ void
+ Interface<dim>::
+ declare_parameters (dealii::ParameterHandler &prm)
+ {}
+
+
+ template <int dim>
+ void
+ Interface<dim>::parse_parameters (dealii::ParameterHandler &prm)
+ {}
+
+
+// -------------------------------- Deal with registering models and automating
+// -------------------------------- their setup and selection at run time
+
+ namespace
+ {
+ std_cxx1x::tuple
+ <void *,
+ void *,
+ internal::Plugins::PluginList<Interface<2> >,
+ internal::Plugins::PluginList<Interface<3> > > registered_plugins;
+ }
+
+
+
+ template <int dim>
+ void
+ register_boundary_composition (const std::string &name,
+ const std::string &description,
+ void (*declare_parameters_function) (ParameterHandler &),
+ Interface<dim> *(*factory_function) ())
+ {
+ std_cxx1x::get<dim>(registered_plugins).register_plugin (name,
+ description,
+ declare_parameters_function,
+ factory_function);
+ }
+
+
+ template <int dim>
+ Interface<dim> *
+ create_boundary_composition (ParameterHandler &prm)
+ {
+ std::string model_name;
+ prm.enter_subsection ("Boundary composition model");
+ {
+ model_name = prm.get ("Model name");
+ }
+ prm.leave_subsection ();
+
+ return std_cxx1x::get<dim>(registered_plugins).create_plugin (model_name,
+ "Boundary composition model::Model name",
+ prm);
+ }
+
+
+
+ template <int dim>
+ void
+ declare_parameters (ParameterHandler &prm)
+ {
+ // declare the entry in the parameter file
+ prm.enter_subsection ("Boundary composition model");
+ {
+ const std::string pattern_of_names
+ = std_cxx1x::get<dim>(registered_plugins).get_pattern_of_names ();
+ prm.declare_entry ("Model name", "box",
+ Patterns::Selection (pattern_of_names),
+ "Select one of the following models:
"
+ +
+ std_cxx1x::get<dim>(registered_plugins).get_description_string());
+ }
+ prm.leave_subsection ();
+
+ std_cxx1x::get<dim>(registered_plugins).declare_parameters (prm);
+ }
+
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace internal
+ {
+ namespace Plugins
+ {
+ template <>
+ std::list<internal::Plugins::PluginList<BoundaryComposition::Interface<2> >::PluginInfo> *
+ internal::Plugins::PluginList<BoundaryComposition::Interface<2> >::plugins = 0;
+ template <>
+ std::list<internal::Plugins::PluginList<BoundaryComposition::Interface<3> >::PluginInfo> *
+ internal::Plugins::PluginList<BoundaryComposition::Interface<3> >::plugins = 0;
+ }
+ }
+
+ namespace BoundaryComposition
+ {
+#define INSTANTIATE(dim) \
+ template class Interface<dim>; \
+ \
+ template \
+ void \
+ register_boundary_composition<dim> (const std::string &, \
+ const std::string &, \
+ void ( *) (ParameterHandler &), \
+ Interface<dim> *( *) ()); \
+ \
+ template \
+ void \
+ declare_parameters<dim> (ParameterHandler &); \
+ \
+ template \
+ Interface<dim> * \
+ create_boundary_composition<dim> (ParameterHandler &prm);
+
+ ASPECT_INSTANTIATE(INSTANTIATE)
+ }
+}
Modified: trunk/aspect/source/material_model/interface.cc
===================================================================
--- trunk/aspect/source/material_model/interface.cc 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/source/material_model/interface.cc 2013-11-25 14:36:21 UTC (rev 2053)
@@ -296,6 +296,19 @@
template <int dim>
+ double
+ InterfaceCompatibility<dim>::
+ entropy_derivative (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position,
+ const NonlinearDependence::Dependence dependence) const
+ {
+ return 0.0;
+ }
+
+
+ template <int dim>
void
declare_parameters (ParameterHandler &prm)
{
@@ -344,6 +357,8 @@
specific_heat.resize(n_points);
thermal_conductivities.resize(n_points);
compressibilities.resize(n_points);
+ entropy_derivative_pressure.resize(n_points);
+ entropy_derivative_temperature.resize(n_points);
}
@@ -375,6 +390,8 @@
out.specific_heat[i] = specific_heat (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
out.thermal_conductivities[i] = thermal_conductivity (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
out.compressibilities[i] = compressibility (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.entropy_derivative_pressure[i] = entropy_derivative (in.temperature[i], in.pressure[i], in.composition[i], in.position[i], NonlinearDependence::pressure);
+ out.entropy_derivative_temperature[i] = entropy_derivative (in.temperature[i], in.pressure[i], in.composition[i], in.position[i], NonlinearDependence::temperature);
}
}
}
Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/source/simulator/assembly.cc 2013-11-25 14:36:21 UTC (rev 2053)
@@ -207,6 +207,8 @@
std::vector<double> old_pressure;
std::vector<double> old_old_pressure;
+ std::vector<Tensor<1,dim> > old_pressure_gradients;
+ std::vector<Tensor<1,dim> > old_old_pressure_gradients;
std::vector<SymmetricTensor<2,dim> > old_strain_rates;
std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
@@ -228,6 +230,7 @@
std::vector<Tensor<1,dim> > current_velocity_values;
std::vector<SymmetricTensor<2,dim> > current_strain_rates;
std::vector<double> current_pressure_values;
+ std::vector<Tensor<1,dim> > current_pressure_gradients;
std::vector<std::vector<double> > current_composition_values;
typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
@@ -263,6 +266,8 @@
old_old_velocity_values (quadrature.size()),
old_pressure (quadrature.size()),
old_old_pressure (quadrature.size()),
+ old_pressure_gradients (quadrature.size()),
+ old_old_pressure_gradients (quadrature.size()),
old_strain_rates (quadrature.size()),
old_old_strain_rates (quadrature.size()),
old_temperature_values (quadrature.size()),
@@ -279,6 +284,7 @@
current_velocity_values(quadrature.size()),
current_strain_rates(quadrature.size()),
current_pressure_values(quadrature.size()),
+ current_pressure_gradients(quadrature.size()),
current_composition_values(n_compositional_fields,
std::vector<double>(quadrature.size())),
material_model_inputs(quadrature.size(), n_compositional_fields),
@@ -306,6 +312,8 @@
old_old_velocity_values (scratch.old_old_velocity_values),
old_pressure (scratch.old_pressure),
old_old_pressure (scratch.old_old_pressure),
+ old_pressure_gradients (scratch.old_pressure_gradients),
+ old_old_pressure_gradients (scratch.old_old_pressure_gradients),
old_strain_rates (scratch.old_strain_rates),
old_old_strain_rates (scratch.old_old_strain_rates),
old_temperature_values (scratch.old_temperature_values),
@@ -320,6 +328,7 @@
current_velocity_values(scratch.current_velocity_values),
current_strain_rates(scratch.current_strain_rates),
current_pressure_values(scratch.current_pressure_values),
+ current_pressure_gradients(scratch.current_pressure_gradients),
current_composition_values(scratch.current_composition_values),
material_model_inputs(scratch.material_model_inputs),
material_model_outputs(scratch.material_model_outputs),
@@ -1284,6 +1293,7 @@
const double current_T = material_model_inputs.temperature[q];
const SymmetricTensor<2,dim> current_strain_rate = material_model_inputs.strain_rate[q];
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
+ const Tensor<1,dim> current_grad_p = scratch.current_pressure_gradients[q];
const double alpha = material_model_outputs.thermal_expansion_coefficients[q];
const double density = material_model_outputs.densities[q];
@@ -1294,6 +1304,7 @@
material_model_outputs.compressibilities[q]
:
std::numeric_limits<double>::quiet_NaN() );
+ const double entropy_gradient = material_model_outputs.entropy_derivative_pressure[q];
const Tensor<1,dim>
gravity = gravity_model->gravity_vector (scratch.finite_element_values.quadrature_point(q));
@@ -1329,15 +1340,27 @@
:
0)
+
- // finally add the term from adiabatic compression heating
+ // add the term from adiabatic compression heating
// - drho/dT T (u . g)
// where we use the definition of
// alpha = - 1/rho drho/dT
(parameters.include_adiabatic_heating
?
- alpha * density * (current_u*gravity) * current_T
+ (current_u * current_grad_p) * alpha * current_T
:
0)
+ +
+ // finally add the right-hand side term from latent heating
+ // DeltaS dLambda/dpi T rho (v . grad p)
+ // DeltaS: change of entropy across phase transition
+ // dLambda/dpi: derivative of the phase function
+ // pi: excess pressure (argument of the phase function)
+ // formulation modified after Christensen & Yuen, 1985
+ (parameters.include_latent_heat
+ ?
+ current_T * density * entropy_gradient * (current_u * current_grad_p)
+ :
+ 0)
);
return gamma;
@@ -1410,6 +1433,13 @@
scratch.finite_element_values[introspection.extractors.pressure].get_function_values (old_old_solution,
scratch.old_old_pressure);
+ scratch.finite_element_values[introspection.extractors.pressure].get_function_gradients (old_solution,
+ scratch.old_pressure_gradients);
+ scratch.finite_element_values[introspection.extractors.pressure].get_function_gradients (old_old_solution,
+ scratch.old_old_pressure_gradients);
+ scratch.finite_element_values[introspection.extractors.pressure].get_function_gradients (current_linearization_point,
+ scratch.current_pressure_gradients);
+
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
scratch.finite_element_values[introspection.extractors.compositional_fields[c]].get_function_values(old_solution,
@@ -1504,6 +1534,14 @@
scratch.material_model_outputs.thermal_conductivities[q]
:
0.0);
+ const double latent_heat_LHS =
+ (parameters.include_latent_heat && temperature_or_composition.is_temperature())
+ ?
+ scratch.material_model_outputs.densities[q] *
+ scratch.material_model_inputs.temperature[q] *
+ scratch.material_model_outputs.entropy_derivative_temperature[q]
+ :
+ 0.0;
const double gamma = compute_heating_term(scratch,
scratch.material_model_inputs,
scratch.material_model_outputs,
@@ -1521,7 +1559,7 @@
:
(*scratch.old_field_values)[q])
*
- density_c_P;
+ (density_c_P + latent_heat_LHS);
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
@@ -1548,7 +1586,7 @@
* (scratch.grad_phi_field[i] * scratch.grad_phi_field[j]))
+ ((time_step * (current_u * scratch.grad_phi_field[j] * scratch.phi_field[i]))
+ (factor * scratch.phi_field[i] * scratch.phi_field[j])) *
- density_c_P
+ (density_c_P + latent_heat_LHS)
)
* scratch.finite_element_values.JxW(q);
}
@@ -1701,8 +1739,8 @@
template void Simulator<dim>::copy_local_to_global_advection_system ( \
const internal::Assembly::CopyData::AdvectionSystem<dim> &data); \
template void Simulator<dim>::assemble_advection_system (const TemperatureOrComposition &temperature_or_composition); \
-
+
ASPECT_INSTANTIATE(INSTANTIATE)
}
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2013-11-25 14:30:03 UTC (rev 2052)
+++ trunk/aspect/source/simulator/core.cc 2013-11-25 14:36:21 UTC (rev 2053)
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+ Copyright (C) 2011, 2012, 2013 by the authors of the ASPECT code.
This file is part of ASPECT.
@@ -99,6 +99,7 @@
material_model (MaterialModel::create_material_model<dim>(prm)),
gravity_model (GravityModel::create_gravity_model<dim>(prm)),
boundary_temperature (BoundaryTemperature::create_boundary_temperature<dim>(prm)),
+ boundary_composition (BoundaryComposition::create_boundary_composition<dim>(prm)),
compositional_initial_conditions (CompositionalInitialConditions::create_initial_conditions (prm,
*geometry_model)),
adiabatic_conditions(),
@@ -183,14 +184,22 @@
ExcMessage ("One of the boundary indicators listed in the input file "
"is not used by the geometry model."));
- // now do the same for the fixed temperature indicators
+ // now do the same for the fixed temperature indicators and the
+ // compositional indicators
for (typename std::set<types::boundary_id>::const_iterator
p = parameters.fixed_temperature_boundary_indicators.begin();
p != parameters.fixed_temperature_boundary_indicators.end(); ++p)
AssertThrow (all_boundary_indicators.find (*p)
!= all_boundary_indicators.end(),
- ExcMessage ("One of the fixed boundary indicators listed in the input file "
+ ExcMessage ("One of the fixed boundary temperature indicators listed in the input file "
"is not used by the geometry model."));
+ for (typename std::set<types::boundary_id>::const_iterator
+ p = parameters.fixed_composition_boundary_indicators.begin();
+ p != parameters.fixed_composition_boundary_indicators.end(); ++p)
+ AssertThrow (all_boundary_indicators.find (*p)
+ != all_boundary_indicators.end(),
+ ExcMessage ("One of the fixed boundary composition indicators listed in the input file "
+ "is not used by the geometry model."));
}
// continue with initializing members that can't be initialized for one reason
@@ -205,6 +214,8 @@
sim->initialize (*this);
if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*boundary_temperature))
sim->initialize (*this);
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*boundary_composition))
+ sim->initialize (*this);
adiabatic_conditions.reset(new AdiabaticConditions<dim> (*geometry_model,
*gravity_model,
@@ -676,10 +687,10 @@
constraints.reinit(introspection.index_sets.system_relevant_set);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
-
+
//Now set up the constraints for periodic boundary conditions
{
- typedef std::set< std::pair< std::pair< types::boundary_id, types::boundary_id>, unsigned int> >
+ typedef std::set< std::pair< std::pair< types::boundary_id, types::boundary_id>, unsigned int> >
periodic_boundary_set;
periodic_boundary_set pbs = geometry_model->get_periodic_boundary_pairs();
@@ -695,19 +706,19 @@
parameters.prescribed_velocity_boundary_indicators.find( (*p).first.first)
== parameters.prescribed_velocity_boundary_indicators.end() &&
parameters.prescribed_velocity_boundary_indicators.find( (*p).first.second)
- == parameters.prescribed_velocity_boundary_indicators.end(),
+ == parameters.prescribed_velocity_boundary_indicators.end(),
ExcInternalError());
#if (DEAL_II_MAJOR*100 + DEAL_II_MINOR) >= 801
- DoFTools::make_periodicity_constraints(dof_handler,
+ DoFTools::make_periodicity_constraints(dof_handler,
(*p).first.first, //first boundary id
(*p).first.second, //second boundary id
(*p).second, //cartesian direction for translational symmetry
constraints);
#endif
}
-
+
}
// then compute constraints for the velocity. the constraints we compute
// here are the ones that are the same for all following time steps. in
@@ -762,8 +773,29 @@
}
- // we do nothing with the compositional fields: homogeneous Neumann boundary conditions
-
More information about the CIG-COMMITS
mailing list