add multicomponent equation of state

include/aspect/material_model/dynamic_friction.h

@@ -23,7 +23,7 @@
 
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
-#include <deal.II/base/parameter_handler.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 
 namespace aspect
 {
@@ -137,10 +137,7 @@ namespace aspect
 
         MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;
 
-        /**
-         * Vector for field densities, read from parameter file .
-         */
-        std::vector<double> densities;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
 
         /**
          * The dynamic coefficient of friction
@@ -174,21 +171,10 @@ namespace aspect
         double reference_strain_rate;
         double minimum_strain_rate;
 
-        /**
-         * Vector for field thermal expnsivities, read from parameter file.
-         */
-        std::vector<double> thermal_expansivities;
-
         /**
          * Vector for field thermal conductivities, read from parameter file.
          */
         std::vector<double> thermal_conductivities;
-
-        /**
-         * Vector for field specific heats, read from parameter file.
-         */
-        std::vector<double> specific_heats;
-
     };
 
   }

include/aspect/material_model/equation_of_state/multicomponent_incompressible.h

@@ -0,0 +1,118 @@
+/*
+  Copyright (C) 2011 - 2019 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 LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _aspect_material_model_equation_of_state_multicomponent_incompressible_h
+#define _aspect_material_model_equation_of_state_multicomponent_incompressible_h
+
+#include <aspect/material_model/interface.h>
+#include <aspect/simulator_access.h>
+#include <aspect/material_model/equation_of_state/interface.h>
+
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    namespace EquationOfState
+    {
+      using namespace dealii;
+
+      /**
+       * An incompressible equation of state that is intended for use with multiple compositional
+       * fields. For each material property, the user supplies a comma delimited list of
+       * length N+1, where N is the number of compositional fields used in the computation.
+       * The first entry corresponds to the "background" (which is also why there are N+1 entries).
+       *
+       * If a single value is given, then all the compositional fields are given
+       * that value. Other lengths of lists are not allowed. For a given
+       * compositional field the material parameters are treated as constant,
+       * except density, which varies linearly with temperature according to the equation:
+       *
+       * $\rho(p,T,\mathfrak c) = \left(1-\alpha_i (T-T_0)\right) \rho_0(\mathfrak c_i).$
+       *
+       * There is no pressure-dependence of the density.
+       */
+      template <int dim>
+      class MulticomponentIncompressible :  public ::aspect::SimulatorAccess<dim>
+      {
+        public:
+          /**
+           * A function that computes the output of the equation of state @p out
+           * for all compositions, given the inputs in @p in and an index q that
+           * determines which entry of the vector of inputs is used.
+           */
+          void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+                        const unsigned int q,
+                        MaterialModel::EquationOfStateOutputs<dim> &out) const;
+
+          /**
+           * Return whether the model is compressible or not. Incompressibility
+           * does not necessarily imply that the density is constant; rather, it
+           * may still depend on temperature or pressure. In the current
+           * context, compressibility means whether we should solve the continuity
+           * equation as $\nabla \cdot (\rho \mathbf u)=0$ (compressible Stokes)
+           * or as $\nabla \cdot \mathbf{u}=0$ (incompressible Stokes).
+           */
+          bool is_compressible () const;
+
+          /**
+           * Declare the parameters this class takes through input files.
+           * The optional parameter @p n_compositions determines the maximum
+           * number of compositions the equation of state is set up with,
+           * in other words, how many compositional fields influence the
+           * density.
+           */
+          static
+          void
+          declare_parameters (ParameterHandler &prm,
+                              const double default_thermal_expansion = 3.5e-5);
+
+          /**
+           * Read the parameters this class declares from the parameter file.
+           * The optional parameter @p n_compositions determines the maximum
+           * number of compositions the equation of state is set up with,
+           * and should have the same value as the parameter with the same
+           * name in the declare_parameters() function.
+           */
+          void
+          parse_parameters (ParameterHandler &prm);
+
+
+        private:
+          // Vector of reference densities $\rho_0$ with one entry per composition,
+          // used in the computation of the density.
+          std::vector<double> densities;
+
+          // The reference temperature $T_0$ used in the computation of the density.
+          // All component use the same reference temperature.
+          double reference_T;
+
+          // Vector of constant thermal expansivities $\alpha$ with one entry per composition,
+          // used in the computation of the density.
+          std::vector<double> thermal_expansivities;
+
+          // Vector of specific heat capacities with one entry per composition.
+          std::vector<double> specific_heats;
+      };
+    }
+  }
+}
+
+#endif

include/aspect/material_model/multicomponent.h

@@ -23,7 +23,8 @@
 
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
-#include <deal.II/base/parameter_handler.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
+
 
 namespace aspect
 {
@@ -128,30 +129,17 @@ namespace aspect
          */
         MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;
 
-        /**
-         * Vector for field densities, read from parameter file .
-         */
-        std::vector<double> densities;
-
         /**
          * Vector for field viscosities, read from parameter file.
          */
         std::vector<double> viscosities;
 
-        /**
-         * Vector for field thermal expansivities, read from parameter file.
-         */
-        std::vector<double> thermal_expansivities;
-
         /**
          * Vector for field thermal conductivities, read from parameter file.
          */
         std::vector<double> thermal_conductivities;
 
-        /**
-         * Vector for field specific heats, read from parameter file.
-         */
-        std::vector<double> specific_heats;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
     };
 
   }

include/aspect/material_model/visco_plastic.h

@@ -23,6 +23,7 @@
 
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 
 #include<deal.II/fe/component_mask.h>
 
@@ -178,10 +179,9 @@ namespace aspect
 
         double grain_size;
 
-        std::vector<double> densities;
-        std::vector<double> thermal_expansivities;
         std::vector<double> thermal_diffusivities;
-        std::vector<double> heat_capacities;
+
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
 
         /**
          * Enumeration for selecting which viscosity averaging scheme to use.

include/aspect/material_model/viscoelastic.h

@@ -23,6 +23,7 @@
 
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 
 namespace aspect
 {
@@ -249,31 +250,18 @@ namespace aspect
                                                          const double dte) const;
 
 
-        /**
-         * Vector for field densities, read from parameter file.
-         */
-        std::vector<double> densities;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
 
         /**
          * Vector for field viscosities, read from parameter file.
          */
         std::vector<double> viscosities;
 
-        /**
-         * Vector for field thermal expnsivities, read from parameter file.
-         */
-        std::vector<double> thermal_expansivities;
-
         /**
          * Vector for field thermal conductivities, read from parameter file.
          */
         std::vector<double> thermal_conductivities;
 
-        /**
-         * Vector for field specific heats, read from parameter file.
-         */
-        std::vector<double> specific_heats;
-
         /**
          * Vector for field elastic shear moduli, read from parameter file.
          */

include/aspect/material_model/viscoelastic_plastic.h

@@ -24,6 +24,7 @@
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/viscoelastic.h>
+#include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
 
 namespace aspect
 {
@@ -268,26 +269,13 @@ namespace aspect
          */
         MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;
 
-        /**
-         * Vector for field densities, read from parameter file.
-         */
-        std::vector<double> densities;
-
-        /**
-         * Vector for field thermal expansivities, read from parameter file.
-         */
-        std::vector<double> thermal_expansivities;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
 
         /**
          * Vector for field thermal conductivities, read from parameter file.
          */
         std::vector<double> thermal_conductivities;
 
-        /**
-         * Vector for field specific heats, read from parameter file.
-         */
-        std::vector<double> specific_heats;
-
         /**
          * Vector for field linear (fixed) viscosities, read from parameter file.
          */