Merge pull request geodynamics#1765 from jdannberg/melt_free_surface

Make melt models work with free surface

Author: gassmoeller

doc/modules/changes/20180130_jdannberg

@@ -0,0 +1,4 @@
+New: Aspect now supports using the free surface in models with 
+melt migration, and with all nonlinear solver schemes.
+<br>
+(Juliane Dannberg, 2018/01/30)

include/aspect/melt.h

@@ -326,6 +326,15 @@ namespace aspect
        */
       bool is_porosity (const typename Simulator<dim>::AdvectionField &advection_field) const;
 
+      /**
+       * Apply free surface stabilization to a cell of the system matrix when melt
+       * transport is used in the computation. Called during assembly of the system matrix.
+       */
+      void apply_free_surface_stabilization_with_melt (const double free_surface_theta,
+                                                       const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                                       internal::Assembly::Scratch::StokesSystem<dim>       &scratch,
+                                                       internal::Assembly::CopyData::StokesSystem<dim>      &data) const;
+
       /**
        * The porosity limit for melt migration. For smaller porosities, the equations
        * reduce to the Stokes equations and neglect melt transport. In practice, this

source/simulator/assembly.cc

@@ -237,8 +237,9 @@ namespace aspect
       {
         melt_handler->set_assemblers(*assemblers);
 
-        AssertThrow (parameters.free_surface_enabled == false,
-                     ExcMessage("The melt implementation does not support free surface computations."));
+        // Let the free surface add its assembler:
+        if (parameters.free_surface_enabled)
+          free_surface->set_assemblers();
       }
     else if (!parameters.include_melt_transport
              && assemble_newton_stokes_system)
@@ -437,9 +438,9 @@ namespace aspect
 
     /*  The stabilization term for the free surface (Kaus et. al., 2010)
      *  makes changes to the system matrix which are of the same form as
-     *  boundary stresses.  If these stresses are not also added to the
-     *  system_preconditioner_matrix, then  if fails to be very good as a
-     *  preconditioner.  Instead, we just pass the system_matrix to the
+     *  boundary stresses. If these stresses are not also added to the
+     *  system_preconditioner_matrix, then it fails to be very good as a
+     *  preconditioner. Instead, we just pass the system_matrix to the
      *  AMG precondition initialization so that it builds the preconditioner
      *  directly from that. However, we still need the mass matrix for the
      *  pressure block which is assembled in the preconditioner matrix.

source/simulator/core.cc

@@ -548,8 +548,6 @@ namespace aspect
         AssertThrow( !parameters.use_discontinuous_temperature_discretization &&
                      !parameters.use_discontinuous_composition_discretization,
                      ExcMessage("Melt transport can not be used with discontinuous elements.") );
-        AssertThrow( !parameters.free_surface_enabled,
-                     ExcMessage("Melt transport together with a free surface has not been tested.") );
         melt_handler->initialize_simulator (*this);
       }
 

source/simulator/free_surface.cc

@@ -22,6 +22,8 @@
 #include <aspect/simulator.h>
 #include <aspect/free_surface.h>
 #include <aspect/global.h>
+#include <aspect/simulator/assemblers/interface.h>
+#include <aspect/melt.h>
 
 #include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_accessor.h>
@@ -52,6 +54,15 @@ namespace aspect
       if (!this->get_parameters().free_surface_enabled)
         return;
 
+      if (this->get_parameters().include_melt_transport)
+        {
+          this->get_melt_handler().apply_free_surface_stabilization_with_melt (this->get_free_surface_handler().get_stabilization_term(),
+                                                                               scratch.cell,
+                                                                               scratch,
+                                                                               data);
+          return;
+        }
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
 

source/simulator/melt.cc

@@ -475,6 +475,9 @@ namespace aspect
       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);
 
+      const types::boundary_id boundary_indicator
+        = scratch.cell->face(scratch.face_number)->boundary_id();
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
 
@@ -1396,6 +1399,88 @@ namespace aspect
       std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputs<dim> >
       (new MaterialModel::MeltOutputs<dim> (n_points, n_comp)));
   }
+
+
+  template <int dim>
+  void
+  MeltHandler<dim>::
+  apply_free_surface_stabilization_with_melt (const double free_surface_theta,
+                                              const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                              internal::Assembly::Scratch::StokesSystem<dim>       &scratch,
+                                              internal::Assembly::CopyData::StokesSystem<dim>      &data) const
+  {
+    if (!this->get_parameters().free_surface_enabled)
+      return;
+
+    const unsigned int n_face_q_points = scratch.face_finite_element_values.n_quadrature_points;
+    const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
+
+    // only apply on free surface faces
+    if (cell->at_boundary() && cell->is_locally_owned())
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+        if (cell->face(face_no)->at_boundary())
+          {
+            const types::boundary_id boundary_indicator
+              = cell->face(face_no)->boundary_id();
+
+            if (this->get_parameters().free_surface_boundary_indicators.find(boundary_indicator)
+                == this->get_parameters().free_surface_boundary_indicators.end())
+              continue;
+
+            scratch.face_finite_element_values.reinit(cell, face_no);
+
+            this->compute_material_model_input_values (this->get_solution(),
+                                                       scratch.face_finite_element_values,
+                                                       cell,
+                                                       true,
+                                                       scratch.face_material_model_inputs);
+
+            this->get_material_model().evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs);
+
+            const unsigned int p_f_component_index = this->introspection().variable("fluid pressure").first_component_index;
+            const unsigned int p_c_component_index = this->introspection().variable("compaction pressure").first_component_index;
+
+            for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
+              {
+                for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
+                  {
+                    const unsigned int component_index_i = this->get_fe().system_to_component_index(i).first;
+
+                    if (Assemblers::is_velocity_or_pressures(this->introspection(),p_c_component_index,p_f_component_index,component_index_i))
+                      {
+                        scratch.phi_u[i_stokes] = scratch.face_finite_element_values[this->introspection().extractors.velocities].value(i, q_point);
+                        ++i_stokes;
+                      }
+                    ++i;
+                  }
+
+                const Tensor<1,dim>
+                gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point));
+                const double g_norm = gravity.norm();
+
+                // construct the relevant vectors
+                const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point);
+                const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm);
+
+                const double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] *
+                                                     this->get_timestep() * free_surface_theta * g_norm;
+
+                // see Kaus et al 2010 for details of the stabilization term
+                for (unsigned int i=0; i< stokes_dofs_per_cell; ++i)
+                  for (unsigned int j=0; j< stokes_dofs_per_cell; ++j)
+                    {
+                      // The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n)
+                      const double stress_value = - pressure_perturbation
+                                                  * (scratch.phi_u[i] * g_hat) * (scratch.phi_u[j] * n_hat)
+                                                  * scratch.face_finite_element_values.JxW(q_point);
+
+                      data.local_matrix(i,j) += stress_value;
+                    }
+              }
+          }
+
+
+  }
 }