Add AdditionalMaterialOutputsStokesRHS

doc/modules/changes/20170509_tjhei

@@ -0,0 +1,5 @@
+New: Allow additional RHS force terms in the Stokes system by enabling the
+parameter "Enable additional Stokes RHS" and filling
+AdditionalMaterialOutputsStokesRHS in the material model.
+<br>
+(Timo Heister, 2017/05/09)

include/aspect/assembly.h

@@ -711,6 +711,9 @@ namespace aspect
       public SimulatorAccess<dim>
     {
       public:
+
+        virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const;
+
         /**
          * This function assembles the terms of the Stokes preconditioner matrix for the current cell.
          */

include/aspect/material_model/interface.h

@@ -589,6 +589,53 @@ namespace aspect
         std::vector<double> vp;
     };
 
+
+    /**
+     * A class for additional output fields to be added to the RHS of the
+     * Stokes system, which can be attached to the
+     * MaterialModel::MaterialModelOutputs structure and filled in the
+     * MaterialModel::Interface::evaluate() function.
+     */
+    template<int dim>
+    class AdditionalMaterialOutputsStokesRHS: public AdditionalMaterialOutputs<dim>
+    {
+      public:
+        AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
+          : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
+        {}
+
+        virtual ~AdditionalMaterialOutputsStokesRHS()
+        {}
+
+        virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
+                              const FullMatrix<double>  &/*projection_matrix*/,
+                              const FullMatrix<double>  &/*expansion_matrix*/)
+        {
+          // TODO: not implemented
+        }
+
+        /**
+         * Force tensor on the right-hand side for the conservation of
+         * momentum equation (first part of the Stokes equation) in each
+         * quadrature point.
+         */
+        std::vector<Tensor<1,dim> > rhs_u;
+
+        /**
+         * Force value for the conservation of mass equation (second Stokes
+        * equation) in each quadrature point.
+        */
+        std::vector<double> rhs_p;
+
+        /**
+        * Force for the compaction pressure equation (when using melt
+        * transport) in each quadrature point.
+        */
+        std::vector<double> rhs_melt_pc;
+    };
+
+
+
     /**
      * A base class for parameterizations of material models. Classes derived
      * from this class will need to implement functions that provide material

include/aspect/parameters.h

@@ -367,6 +367,7 @@ namespace aspect
      * @{
      */
     bool                           include_melt_transport;
+    bool                           enable_additional_stokes_rhs;
 
     std::set<types::boundary_id> fixed_temperature_boundary_indicators;
     std::set<types::boundary_id> fixed_composition_boundary_indicators;

source/simulator/assemblers/stokes.cc

@@ -22,11 +22,32 @@
 #include <aspect/utilities.h>
 #include <aspect/assembly.h>
 #include <aspect/simulator_access.h>
+#include <deal.II/base/signaling_nan.h>
 
 namespace aspect
 {
   namespace Assemblers
   {
+    template <int dim>
+    void
+    StokesAssembler<dim>::
+    create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
+    {
+      const unsigned int n_points = outputs.viscosities.size();
+
+      if (this->get_parameters().enable_additional_stokes_rhs
+          && outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >() == NULL)
+        {
+          outputs.additional_outputs.push_back(
+            std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >
+            (new MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> (n_points)));
+        }
+      Assert(!this->get_parameters().enable_additional_stokes_rhs
+             ||
+             outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()->rhs_u.size()
+             == n_points, ExcInternalError());
+    }
+
     template <int dim>
     void
     StokesAssembler<dim>::
@@ -113,6 +134,9 @@ namespace aspect
       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 MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>
+      *force = scratch.material_model_outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >();
+
       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*/)
@@ -147,14 +171,18 @@ namespace aspect
           gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(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) += (density * gravity * scratch.phi_u[i])
                                    * JxW;
 
+              if (force != NULL)
+                data.local_rhs(i) += (force->rhs_u[q] * scratch.phi_u[i]
+                                      + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i])
+                                     * JxW;
+
               if (rebuild_stokes_matrix)
                 for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
                   {

source/simulator/assembly.cc

@@ -500,6 +500,36 @@ namespace aspect
 
   namespace Assemblers
   {
+    virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
+    {
+      if (this->get_parameters().enable_additional_stokes_rhs
+          && !outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >())
+        {
+          const unsigned int n_points = outputs.viscosities.size();
+          const unsigned int n_stokes_comp = dim + 2;
+          outputs.additional_outputs.push_back(
+            std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >
+            (new MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> (n_points, n_stokes_comp)));
+        }
+    }
+
+    MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>
+    *force = scratch.material_model_outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >();
+
+    // fill values from Stokes RHS additional outputs
+    Tensor<1,dim> force_u = numbers::signaling_nan<Tensor<1,dim> >();
+    if (force != NULL)
+      {
+        for (unsigned int d=0; d<dim; ++d)
+          force_u[d] = force->rhs[q][d];
+      }
+    const double force_p = (force!=NULL) ?
+                           force->rhs[q][dim]
+                           : std::numeric_limits<double>::quiet_NaN();
+    if (force != NULL)
+      data.local_rhs(i) += (force_u *scratch.phi_u[i]
+                            + pressure_scaling *force_p *scratch.phi_p[i])
+                           *JxW;
     /**
      * A namespace for the definition of functions that implement various
      * other terms that need to occasionally or always be assembled.

source/simulator/melt.cc

@@ -21,6 +21,7 @@
 
 #include <aspect/melt.h>
 #include <aspect/simulator.h>
+#include <deal.II/base/signaling_nan.h>
 
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/lac/sparsity_tools.h>
@@ -88,6 +89,21 @@ namespace aspect
     MeltEquations<dim>::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
     {
       MeltHandler<dim>::create_material_model_outputs(outputs);
+
+      const unsigned int n_points = outputs.viscosities.size();
+
+      if (this->get_parameters().enable_additional_stokes_rhs
+          && outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >() == NULL)
+        {
+          outputs.additional_outputs.push_back(
+            std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >
+            (new MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> (n_points)));
+        }
+
+      Assert(!this->get_parameters().enable_additional_stokes_rhs
+             ||
+             outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()->rhs_u.size()
+             == n_points, ExcInternalError());
     }
 
 
@@ -215,6 +231,8 @@ namespace aspect
       const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index;
 
       MaterialModel::MeltOutputs<dim> *melt_outputs = scratch.material_model_outputs.template get_additional_output<MaterialModel::MeltOutputs<dim> >();
+      const MaterialModel::AdditionalMaterialOutputsStokesRHS<dim>
+      *force = scratch.material_model_outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >();
 
       for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
         {
@@ -300,6 +318,13 @@ namespace aspect
               data.local_rhs(i) += (
                                      (bulk_density * gravity * scratch.phi_u[i])
                                      +
+                                     // add a force to the RHS if present
+                                     (force!=NULL ?
+                                      (force->rhs_u[q] * scratch.phi_u[i]
+                                       + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]
+                                       + pressure_scaling * force->rhs_melt_pc[q] * scratch.phi_p_c[i])
+                                      : 0.0)
+                                     +
                                      // 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

source/simulator/parameters.cc

@@ -546,6 +546,13 @@ namespace aspect
                          "\n\n"
                          "Note that while more than one operation can be selected it only makes sense to "
                          "pick one rotational and one translational operation.");
+      prm.declare_entry ("Enable additional Stokes RHS", "false",
+                         Patterns::Bool (),
+                         "Whether to ask the material model for additional terms for the right-hand side "
+                         "of the Stokes equation. This can be used to implement force vectors for "
+                         "manufactured solution problems. The assembler for the Stokes system will "
+                         "create an additional output of type AdditionalMaterialOutputsStokesRHS.");
+
     }
     prm.leave_subsection ();
 
@@ -1065,6 +1072,7 @@ namespace aspect
     prm.enter_subsection ("Model settings");
     {
       include_melt_transport = prm.get_bool ("Include melt transport");
+      enable_additional_stokes_rhs = prm.get_bool ("Enable additional Stokes RHS");
 
       {
         nullspace_removal = NullspaceRemoval::none;