implement "hydrostatic compression" mass term

doc/manual/manual.tex

@@ -1882,6 +1882,17 @@ \subsubsection{Mass conservation approximation}
 is defined in the material model. This is the explicit compressible mass equation where 
 the velocity $\textbf{u}$ on the right-hand side is an extrapolated velocity from the last timesteps.
 
+\item
+``hydrostatic compression'':
+\[
+ \nabla \cdot \textbf{u}
+= - \left( \frac{1}{\rho} \frac{\partial \rho}{\partial p} \rho \textbf{g} + \frac{1}{\rho} \frac{\partial \rho}{\partial T} \nabla T \right) \cdot \textbf{u}
+= - \left( \kappa \rho \textbf{g} - \alpha \nabla T \right) \cdot \textbf{u}
+\]
+where $\frac{1}{\rho} \frac{\partial \rho}{\partial p} = \kappa$ is the compressibility,
+$- \frac{1}{\rho}\frac{\partial \rho}{\partial T} = \alpha$ is the thermal expansion coefficient, 
+and both are defined in the material model. 
+
 \item
 ``reference density profile'':
 \[

doc/modules/changes/20171122_tjhei

@@ -0,0 +1,3 @@
+New: Implement new mass conservation formulation "hydrostatic compression".
+<br>
+(Juliane Dannberg, Rene Gassmoeller, Timo Heister, 2017/11/22)

include/aspect/material_model/interface.h

@@ -372,7 +372,8 @@ namespace aspect
       std::vector<double> densities;
 
       /**
-       * Thermal expansion coefficients at the given positions.
+       * Thermal expansion coefficients at the given positions. It is defined
+       * as $\alpha = - \frac{1}{\rho} \frac{\partial\rho}{\partial T}$
        */
       std::vector<double> thermal_expansion_coefficients;
 
@@ -387,8 +388,8 @@ namespace aspect
       std::vector<double> thermal_conductivities;
 
       /**
-       * Compressibility at the given positions. The compressibility is given
-       * as $\frac 1\rho \frac{\partial\rho}{\partial p}$.
+       * Compressibility at the given positions. The compressibility is defined
+       * as $\kappa = \frac{1}{\rho} \frac{\partial\rho}{\partial p}$.
        */
       std::vector<double> compressibilities;
 

include/aspect/parameters.h

@@ -167,6 +167,7 @@ namespace aspect
         enum Kind
         {
           isothermal_compression,
+          hydrostatic_compression,
           reference_density_profile,
           implicit_reference_density_profile,
           incompressible,
@@ -182,6 +183,8 @@ namespace aspect
         {
           if (input == "isothermal compression")
             return Formulation::MassConservation::isothermal_compression;
+          else if (input == "hydrostatic compression")
+            return Formulation::MassConservation::hydrostatic_compression;
           else if (input == "reference density profile")
             return Formulation::MassConservation::reference_density_profile;
           else if (input == "implicit reference density profile")

include/aspect/simulator/assemblers/interface.h

@@ -163,6 +163,7 @@ namespace aspect
           std::vector<Tensor<1,dim> >          phi_u;
           std::vector<Tensor<1,dim> >          velocity_values;
           std::vector<double>                  velocity_divergence;
+          std::vector<Tensor<1,dim> >          temperature_gradients;
 
           /**
            * Material model inputs and outputs computed at the current

include/aspect/simulator/assemblers/stokes.h

@@ -120,7 +120,7 @@ namespace aspect
      * includes this term implicitly in the matrix,
      * which is therefore not longer symmetric.
      * This class approximates this term as
-     * $ - \nabla \mathbf{u} - \frac{1}{\rho^{\ast}} * \frac{\partial rho{^\ast}}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$
+     * $ - \nabla \cdot \mathbf{u} - \frac{1}{\rho^{\ast}} * \frac{\partial rho{^\ast}}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$
      */
     template <int dim>
     class StokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
@@ -137,7 +137,7 @@ namespace aspect
      * This class assembles the right-hand-side term of the Stokes equation
      * that is caused by the compressibility in the mass conservation equation.
      * This class approximates this term as
-     * $ - \nabla \mathbf{u} = \kappa \rho \mathbf{g} \cdot \mathbf{u}$
+     * $ - \nabla \cdot \mathbf{u} = \kappa \rho \mathbf{g} \cdot \mathbf{u}$
      * where $\kappa$ is the compressibility provided by the material model,
      * which is frequently computed as
      * $\kappa = \frac{1}{\rho} * \frac{\partial rho}{\partial p}$.
@@ -153,6 +153,28 @@ namespace aspect
                 internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
     };
 
+
+    /**
+     * This class assembles the right-hand-side term of the Stokes equation
+     * that is caused by the compressibility in the mass conservation equation.
+     * This class approximates this term as
+     * $ -\nabla \cdot \mathbf{u} =  - \left( \kappa \rho \textbf{g} - \alpha \nabla T \right) \cdot \textbf{u}$
+     *
+     * where $\frac{1}{\rho} \frac{\partial \rho}{\partial p} = \kappa$ is the compressibility,
+     * $- \frac{1}{\rho}\frac{\partial \rho}{\partial T} = \alpha$ is the thermal expansion coefficient,
+     * and both are defined in the material model.
+     */
+    template <int dim>
+    class StokesHydrostaticCompressionTerm : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        virtual
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
+    };
+
     /**
      * This class assembles the right-hand-side terms that are used to weakly
      * prescribe the boundary tractions.

source/simulator/assemblers/interface.cc

@@ -125,6 +125,7 @@ namespace aspect
           phi_u (stokes_dofs_per_cell, numbers::signaling_nan<Tensor<1,dim> >()),
           velocity_values (quadrature.size(), numbers::signaling_nan<Tensor<1,dim> >()),
           velocity_divergence(quadrature.size(), numbers::signaling_nan<double>()),
+          temperature_gradients (quadrature.size(), numbers::signaling_nan<Tensor<1,dim> >()),
           face_material_model_inputs(face_quadrature.size(), n_compositional_fields),
           face_material_model_outputs(face_quadrature.size(), n_compositional_fields),
           reference_densities(use_reference_density_profile ? quadrature.size() : 0, numbers::signaling_nan<double>()),
@@ -148,6 +149,7 @@ namespace aspect
           phi_u (scratch.phi_u),
           velocity_values (scratch.velocity_values),
           velocity_divergence (scratch.velocity_divergence),
+          temperature_gradients (scratch.temperature_gradients),
           face_material_model_inputs(scratch.face_material_model_inputs),
           face_material_model_outputs(scratch.face_material_model_outputs),
           reference_densities(scratch.reference_densities),

source/simulator/assemblers/stokes.cc

@@ -447,7 +447,11 @@ namespace aspect
       internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>& > (data_base);
 
       // assemble RHS of:
-      //  - div u = 1/rho * drho/dp rho * g * u
+      //  - div \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}
+      // to the manual, this term seems to have the wrong sign, but this
+      // is because we negate the entire equation to make sure we get
+      // -div(u) as the adjoint operator of grad(p)
+
       Assert(this->get_parameters().formulation_mass_conservation ==
              Parameters<dim>::Formulation::MassConservation::isothermal_compression,
              ExcInternalError());
@@ -481,10 +485,6 @@ namespace aspect
 
           for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
             data.local_rhs(i) += (
-                                   // add the term that results from the compressibility. compared
-                                   // to the manual, this term seems to have the wrong sign, but this
-                                   // is because we negate the entire equation to make sure we get
-                                   // -div(u) as the adjoint operator of grad(p)
                                    (pressure_scaling *
                                     compressibility * density *
                                     (scratch.velocity_values[q] * gravity) *
@@ -495,6 +495,76 @@ namespace aspect
     }
 
 
+    template <int dim>
+    void
+    StokesHydrostaticCompressionTerm<dim>::
+    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
+             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
+    {
+      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>& > (scratch_base);
+      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>& > (data_base);
+
+      // assemble RHS of:
+      //  - div \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}
+      // + \frac{1}{\rho} * \frac{\partial rho}{\partial T} \nabla T \cdot \mathbf{u}
+
+      // Compared
+      // to the manual, this term seems to have the wrong sign, but this
+      // is because we negate the entire equation to make sure we get
+      // -div(u) as the adjoint operator of grad(p)
+
+      Assert(this->get_parameters().formulation_mass_conservation ==
+             Parameters<dim>::Formulation::MassConservation::hydrostatic_compression,
+             ExcInternalError());
+
+      const Introspection<dim> &introspection = this->introspection();
+      const FiniteElement<dim> &fe = this->get_fe();
+      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
+      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
+      const double pressure_scaling = this->get_pressure_scaling();
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
+            {
+              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
+                {
+                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
+                  ++i_stokes;
+                }
+              ++i;
+            }
+
+          const Tensor<1,dim>
+          gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
+
+          const double compressibility
+            = scratch.material_model_outputs.compressibilities[q];
+
+          const double thermal_alpha
+            = scratch.material_model_outputs.thermal_expansion_coefficients[q];
+
+          const double density = scratch.material_model_outputs.densities[q];
+          const double JxW = scratch.finite_element_values.JxW(q);
+
+          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
+            data.local_rhs(i) += (
+                                   (pressure_scaling *
+                                    (
+                                      // pressure part:
+                                      compressibility * density *
+                                      (scratch.velocity_values[q] * gravity)
+                                      // temperature part:
+                                      - thermal_alpha *
+                                      (scratch.velocity_values[q] * scratch.temperature_gradients[q])
+                                    ) * scratch.phi_p[i])
+                                 )
+                                 * JxW;
+
+        }
+    }
+
+
     template <int dim>
     void
     StokesPressureRHSCompatibilityModification<dim>::execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
@@ -588,6 +658,7 @@ namespace aspect
   template class StokesReferenceDensityCompressibilityTerm<dim>; \
   template class StokesImplicitReferenceDensityCompressibilityTerm<dim>; \
   template class StokesIsothermalCompressionTerm<dim>; \
+  template class StokesHydrostaticCompressionTerm<dim>; \
   template class StokesPressureRHSCompatibilityModification<dim>; \
   template class StokesBoundaryTraction<dim>;
 

source/simulator/assembly.cc

@@ -248,15 +248,25 @@ namespace aspect
           {
             // do nothing, because we assembled div u =0 above already
           }
-        else
+        else if (parameters.formulation_mass_conservation ==
+                 Parameters<dim>::Formulation::MassConservation::isothermal_compression)
           {
             aspect::Assemblers::StokesIsothermalCompressionTerm<dim> *compressible_terms
               = new aspect::Assemblers::StokesIsothermalCompressionTerm<dim>();
 
             assemblers->stokes_system.push_back(
               std_cxx11::unique_ptr<aspect::Assemblers::StokesIsothermalCompressionTerm<dim> > (compressible_terms));
           }
-
+        else if (parameters.formulation_mass_conservation ==
+                 Parameters<dim>::Formulation::MassConservation::hydrostatic_compression)
+          {
+            aspect::Assemblers::StokesHydrostaticCompressionTerm<dim> *compressible_terms
+              = new aspect::Assemblers::StokesHydrostaticCompressionTerm<dim>();
+            assemblers->stokes_system.push_back(
+              std_cxx11::unique_ptr<aspect::Assemblers::StokesHydrostaticCompressionTerm<dim> > (compressible_terms));
+          }
+        else
+          Assert(false, ExcNotImplemented());
       }
 
     // add the boundary integral for melt migration
@@ -617,6 +627,10 @@ namespace aspect
         scratch.velocity_values);
     if (assemble_newton_stokes_system)
       scratch.finite_element_values[introspection.extractors.velocities].get_function_divergences(current_linearization_point,scratch.velocity_divergence);
+    if (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::hydrostatic_compression)
+      scratch.finite_element_values[introspection.extractors.temperature].get_function_gradients(current_linearization_point,
+          scratch.temperature_gradients);
+
 
     const bool use_reference_density_profile = (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::reference_density_profile)
                                                || (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile);

source/simulator/parameters.cc

@@ -457,7 +457,7 @@ namespace aspect
                          "density that depends on temperature and depth and not on the pressure.}");
 
       prm.declare_entry ("Mass conservation", "ask material model",
-                         Patterns::Selection ("incompressible|isothermal compression|"
+                         Patterns::Selection ("incompressible|isothermal compression|hydrostatic compression|"
                                               "reference density profile|implicit reference density profile|"
                                               "ask material model"),
                          "Possible approximations for the density derivatives in the mass "