step-21.h

Go to the documentation of this file.

 
,
 *                            const unsigned int /*component*/ = 0) const override
 *       {
 *         return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Pressureboundaryvalues"></a> 
 * <h4>Pressure boundary values</h4>
 * 
 
 * 
 * The next are pressure boundary values. As mentioned in the introduction,
 * we choose a linear pressure field:
 * 
 * @code
 *     template <int dim>
 *     class PressureBoundaryValues : public Function<dim>
 *     {
 *     public:
 *       PressureBoundaryValues()
 *         : Function<dim>(1)
 *       {}
 *   
 *       virtual double value(const Point<dim> &p,
 *                            const unsigned int /*component*/ = 0) const override
 *       {
 *         return 1 - p[0];
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Saturationboundaryvalues"></a> 
 * <h4>Saturation boundary values</h4>
 * 
 
 * 
 * Then we also need boundary values on the inflow portions of the
 * boundary. The question whether something is an inflow part is decided
 * when assembling the right hand side, we only have to provide a functional
 * description of the boundary values. This is as explained in the
 * introduction:
 * 
 * @code
 *     template <int dim>
 *     class SaturationBoundaryValues : public Function<dim>
 *     {
 *     public:
 *       SaturationBoundaryValues()
 *         : Function<dim>(1)
 *       {}
 *   
 *       virtual double value(const Point<dim> &p,
 *                            const unsigned int /*component*/ = 0) const override
 *       {
 *         if (p[0] == 0)
 *           return 1;
 *         else
 *           return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Initialdata"></a> 
 * <h4>Initial data</h4>
 * 
 
 * 
 * Finally, we need initial data. In reality, we only need initial data for
 * the saturation, but we are lazy, so we will later, before the first time
 * step, simply interpolate the entire solution for the previous time step
 * from a function that contains all vector components.
 *   
 
 * 
 * We therefore simply create a function that returns zero in all
 * components. We do that by simply forward every function to the
 * Functions::ZeroFunction class. Why not use that right away in the places of
 * this program where we presently use the <code>InitialValues</code> class?
 * Because this way it is simpler to later go back and choose a different
 * function for initial values.
 * 
 * @code
 *     template <int dim>
 *     class InitialValues : public Function<dim>
 *     {
 *     public:
 *       InitialValues()
 *         : Function<dim>(dim + 2)
 *       {}
 *   
 *       virtual double value(const Point<dim>  &p,
 *                            const unsigned int component = 0) const override
 *       {
 *         return Functions::ZeroFunction<dim>(dim + 2).value(p, component);
 *       }
 *   
 *       virtual void vector_value(const Point<dim> &p,
 *                                 Vector<double>   &values) const override
 *       {
 *         Functions::ZeroFunction<dim>(dim + 2).vector_value(p, values);
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Theinversepermeabilitytensor"></a> 
 * <h3>The inverse permeability tensor</h3>
 * 
 
 * 
 * As announced in the introduction, we implement two different permeability
 * tensor fields. Each of them we put into a namespace of its own, so that
 * it will be easy later to replace use of one by the other in the code.
 * 
 
 * 
 * 
 * <a name="step_21-Singlecurvingcrackpermeability"></a> 
 * <h4>Single curving crack permeability</h4>
 * 
 
 * 
 * The first function for the permeability was the one that models a single
 * curving crack. It was already used at the end of @ref step_20 "step-20", and its
 * functional form is given in the introduction of the present tutorial
 * program. As in some previous programs, we have to declare a (seemingly
 * unnecessary) default constructor of the KInverse class to avoid warnings
 * from some compilers:
 * 
 * @code
 *     namespace SingleCurvingCrack
 *     {
 *       template <int dim>
 *       class KInverse : public TensorFunction<2, dim>
 *       {
 *       public:
 *         KInverse()
 *           : TensorFunction<2, dim>()
 *         {}
 *   
 *         virtual void
 *         value_list(const std::vector<Point<dim>> &points,
 *                    std::vector<Tensor<2, dim>>   &values) const override
 *         {
 *           AssertDimension(points.size(), values.size());
 *   
 *           for (unsigned int p = 0; p < points.size(); ++p)
 *             {
 *               values[p].clear();
 *   
 *               const double distance_to_flowline =
 *                 std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
 *   
 *               const double permeability =
 *                 std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
 *                                   (0.1 * 0.1)),
 *                          0.01);
 *   
 *               for (unsigned int d = 0; d < dim; ++d)
 *                 values[p][d][d] = 1. / permeability;
 *             }
 *         }
 *       };
 *     } // namespace SingleCurvingCrack
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Randommediumpermeability"></a> 
 * <h4>Random medium permeability</h4>
 * 
 
 * 
 * This function does as announced in the introduction, i.e. it creates an
 * overlay of exponentials at random places. There is one thing worth
 * considering for this class. The issue centers around the problem that the
 * class creates the centers of the exponentials using a random function. If
 * we therefore created the centers each time we create an object of the
 * present type, we would get a different list of centers each time. That's
 * not what we expect from classes of this type: they should reliably
 * represent the same function.
 *   
 
 * 
 * The solution to this problem is to make the list of centers a static
 * member variable of this class, i.e. there exists exactly one such
 * variable for the entire program, rather than for each object of this
 * type. That's exactly what we are going to do.
 *   
 
 * 
 * The next problem, however, is that we need a way to initialize this
 * variable. Since this variable is initialized at the beginning of the
 * program, we can't use a regular member function for that since there may
 * not be an object of this type around at the time. The C++ standard
 * therefore says that only non-member and static member functions can be
 * used to initialize a static variable. We use the latter possibility by
 * defining a function <code>get_centers</code> that computes the list of
 * center points when called.
 *   
 
 * 
 * Note that this class works just fine in both 2d and 3d, with the only
 * difference being that we use more points in 3d: by experimenting we find
 * that we need more exponentials in 3d than in 2d (we have more ground to
 * cover, after all, if we want to keep the distance between centers roughly
 * equal), so we choose 40 in 2d and 100 in 3d. For any other dimension, the
 * function does presently not know what to do so simply throws an exception
 * indicating exactly this.
 * 
 * @code
 *     namespace RandomMedium
 *     {
 *       template <int dim>
 *       class KInverse : public TensorFunction<2, dim>
 *       {
 *       public:
 *         KInverse()
 *           : TensorFunction<2, dim>()
 *         {}
 *   
 *         virtual void
 *         value_list(const std::vector<Point<dim>> &points,
 *                    std::vector<Tensor<2, dim>>   &values) const override
 *         {
 *           AssertDimension(points.size(), values.size());
 *   
 *           for (unsigned int p = 0; p < points.size(); ++p)
 *             {
 *               values[p].clear();
 *   
 *               double permeability = 0;
 *               for (unsigned int i = 0; i < centers.size(); ++i)
 *                 permeability += std::exp(-(points[p] - centers[i]).norm_square() /
 *                                          (0.05 * 0.05));
 *   
 *               const double normalized_permeability =
 *                 std::min(std::max(permeability, 0.01), 4.);
 *   
 *               for (unsigned int d = 0; d < dim; ++d)
 *                 values[p][d][d] = 1. / normalized_permeability;
 *             }
 *         }
 *   
 *       private:
 *         static std::vector<Point<dim>> centers;
 *   
 *         static std::vector<Point<dim>> get_centers()
 *         {
 *           const unsigned int N =
 *             (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
 *   
 *           std::vector<Point<dim>> centers_list(N);
 *           for (unsigned int i = 0; i < N; ++i)
 *             for (unsigned int d = 0; d < dim; ++d)
 *               centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
 *   
 *           return centers_list;
 *         }
 *       };
 *   
 *   
 *   
 *       template <int dim>
 *       std::vector<Point<dim>> KInverse<dim>::centers =
 *         KInverse<dim>::get_centers();
 *     } // namespace RandomMedium
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Theinversemobilityandsaturationfunctions"></a> 
 * <h3>The inverse mobility and saturation functions</h3>
 * 
 
 * 
 * There are two more pieces of data that we need to describe, namely the
 * inverse mobility function and the saturation curve. Their form is also
 * given in the introduction:
 * 
 * @code
 *     double mobility_inverse(const double S, const double viscosity)
 *     {
 *       return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
 *     }
 *   
 *     double fractional_flow(const double S, const double viscosity)
 *     {
 *       return S * S / (S * S + viscosity * (1 - S) * (1 - S));
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Linearsolversandpreconditioners"></a> 
 * <h3>Linear solvers and preconditioners</h3>
 * 
 
 * 
 * The linear solvers we use are also completely analogous to the ones used
 * in @ref step_20 "step-20". The following classes are therefore copied verbatim from
 * there. Note that the classes here are not only copied from
 * @ref step_20 "step-20", but also duplicate classes in deal.II. In a future version of this
 * example, they should be replaced by an efficient method, though. There is a
 * single change: if the size of a linear system is small, i.e. when the mesh
 * is very coarse, then it is sometimes not sufficient to set a maximum of
 * <code>src.size()</code> CG iterations before the solver in the
 * <code>vmult()</code> function converges. (This is, of course, a result of
 * numerical round-off, since we know that on paper, the CG method converges
 * in at most <code>src.size()</code> steps.) As a consequence, we set the
 * maximum number of iterations equal to the maximum of the size of the linear
 * system and 200.
 * 
 * @code
 *     template <class MatrixType>
 *     class InverseMatrix : public Subscriptor
 *     {
 *     public:
 *       InverseMatrix(const MatrixType &m)
 *         : matrix(&m)
 *       {}
 *   
 *       void vmult(Vector<double> &dst, const Vector<double> &src) const
 *       {
 *         SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
 *                                      1e-8 * src.l2_norm());
 *         SolverCG<Vector<double>> cg(solver_control);
 *   
 *         dst = 0;
 *   
 *         cg.solve(*matrix, dst, src, PreconditionIdentity());
 *       }
 *   
 *     private:
 *       const SmartPointer<const MatrixType> matrix;
 *     };
 *   
 *   
 *   
 *     class SchurComplement : public Subscriptor
 *     {
 *     public:
 *       SchurComplement(const BlockSparseMatrix<double>           &A,
 *                       const InverseMatrix<SparseMatrix<double>> &Minv)
 *         : system_matrix(&A)
 *         , m_inverse(&Minv)
 *         , tmp1(A.block(0, 0).m())
 *         , tmp2(A.block(0, 0).m())
 *       {}
 *   
 *       void vmult(Vector<double> &dst, const Vector<double> &src) const
 *       {
 *         system_matrix->block(0, 1).vmult(tmp1, src);
 *         m_inverse->vmult(tmp2, tmp1);
 *         system_matrix->block(1, 0).vmult(dst, tmp2);
 *       }
 *   
 *     private:
 *       const SmartPointer<const BlockSparseMatrix<double>>           system_matrix;
 *       const SmartPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
 *   
 *       mutable Vector<double> tmp1, tmp2;
 *     };
 *   
 *   
 *   
 *     class ApproximateSchurComplement : public Subscriptor
 *     {
 *     public:
 *       ApproximateSchurComplement(const BlockSparseMatrix<double> &A)
 *         : system_matrix(&A)
 *         , tmp1(A.block(0, 0).m())
 *         , tmp2(A.block(0, 0).m())
 *       {}
 *   
 *       void vmult(Vector<double> &dst, const Vector<double> &src) const
 *       {
 *         system_matrix->block(0, 1).vmult(tmp1, src);
 *         system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
 *         system_matrix->block(1, 0).vmult(dst, tmp2);
 *       }
 *   
 *     private:
 *       const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
 *   
 *       mutable Vector<double> tmp1, tmp2;
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-codeTwoPhaseFlowProblemcodeclassimplementation"></a> 
 * <h3><code>TwoPhaseFlowProblem</code> class implementation</h3>
 * 
 
 * 
 * Here now the implementation of the main class. Much of it is actually
 * copied from @ref step_20 "step-20", so we won't comment on it in much detail. You should
 * try to get familiar with that program first, then most of what is
 * happening here should be mostly clear.
 * 
 
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemTwoPhaseFlowProblem"></a> 
 * <h4>TwoPhaseFlowProblem::TwoPhaseFlowProblem</h4>
 * 
 
 * 
 * First for the constructor. We use @f$RT_k \times DQ_k \times DQ_k@f$
 * spaces. For initializing the DiscreteTime object, we don't set the time
 * step size in the constructor because we don't have its value yet.
 * The time step size is initially set to zero, but it will be computed
 * before it is needed to increment time, as described in a subsection of
 * the introduction. The time object internally prevents itself from being
 * incremented when @f$dt = 0@f$, forcing us to set a non-zero desired size for
 * @f$dt@f$ before advancing time.
 * 
 * @code
 *     template <int dim>
 *     TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem(const unsigned int degree)
 *       : degree(degree)
 *       , fe(FE_RaviartThomas<dim>(degree),
 *            FE_DGQ<dim>(degree),
 *            FE_DGQ<dim>(degree))
 *       , dof_handler(triangulation)
 *       , n_refinement_steps(5)
 *       , time(/*start time*/ 0., /*end time*/ 1.)
 *       , viscosity(0.2)
 *     {}
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemmake_grid_and_dofs"></a> 
 * <h4>TwoPhaseFlowProblem::make_grid_and_dofs</h4>
 * 
 
 * 
 * This next function starts out with well-known functions calls that create
 * and refine a mesh, and then associate degrees of freedom with it. It does
 * all the same things as in @ref step_20 "step-20", just now for three components instead
 * of two.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::make_grid_and_dofs()
 *     {
 *       GridGenerator::hyper_cube(triangulation, 0, 1);
 *       triangulation.refine_global(n_refinement_steps);
 *   
 *       dof_handler.distribute_dofs(fe);
 *       DoFRenumbering::component_wise(dof_handler);
 *   
 *       const std::vector<types::global_dof_index> dofs_per_component =
 *         DoFTools::count_dofs_per_fe_component(dof_handler);
 *       const unsigned int n_u = dofs_per_component[0],
 *                          n_p = dofs_per_component[dim],
 *                          n_s = dofs_per_component[dim + 1];
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl
 *                 << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << " (" << n_u << '+' << n_p << '+' << n_s << ')' << std::endl
 *                 << std::endl;
 *   
 *       const std::vector<types::global_dof_index> block_sizes = {n_u, n_p, n_s};
 *       BlockDynamicSparsityPattern                dsp(block_sizes, block_sizes);
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *   
 *       sparsity_pattern.copy_from(dsp);
 *       system_matrix.reinit(sparsity_pattern);
 *   
 *       solution.reinit(block_sizes);
 *       old_solution.reinit(block_sizes);
 *       system_rhs.reinit(block_sizes);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemassemble_system"></a> 
 * <h4>TwoPhaseFlowProblem::assemble_system</h4>
 * 
 
 * 
 * This is the function that assembles the linear system, or at least
 * everything except the (1,3) block that depends on the still-unknown
 * velocity computed during this time step (we deal with this in
 * <code>assemble_rhs_S</code>). Much of it is again as in @ref step_20 "step-20", but we
 * have to deal with some nonlinearity this time.  However, the top of the
 * function is pretty much as usual (note that we set matrix and right hand
 * side to zero at the beginning &mdash; something we didn't have to do for
 * stationary problems since there we use each matrix object only once and
 * it is empty at the beginning anyway).
 *   
 
 * 
 * Note that in its present form, the function uses the permeability
 * implemented in the RandomMedium::KInverse class. Switching to the single
 * curved crack permeability function is as simple as just changing the
 * namespace name.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::assemble_system()
 *     {
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       const QGauss<dim>     quadrature_formula(degree + 2);
 *       const QGauss<dim - 1> face_quadrature_formula(degree + 2);
 *   
 *       FEValues<dim>     fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       FEFaceValues<dim> fe_face_values(fe,
 *                                        face_quadrature_formula,
 *                                        update_values | update_normal_vectors |
 *                                          update_quadrature_points |
 *                                          update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *   
 *       const unsigned int n_q_points      = quadrature_formula.size();
 *       const unsigned int n_face_q_points = face_quadrature_formula.size();
 *   
 *       FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     local_rhs(dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       const PressureRightHandSide<dim>  pressure_right_hand_side;
 *       const PressureBoundaryValues<dim> pressure_boundary_values;
 *       const RandomMedium::KInverse<dim> k_inverse;
 *   
 *       std::vector<double>         pressure_rhs_values(n_q_points);
 *       std::vector<double>         boundary_values(n_face_q_points);
 *       std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
 *   
 *       std::vector<Vector<double>> old_solution_values(n_q_points,
 *                                                       Vector<double>(dim + 2));
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *       const FEValuesExtractors::Scalar pressure(dim);
 *       const FEValuesExtractors::Scalar saturation(dim + 1);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           local_matrix = 0;
 *           local_rhs    = 0;
 *   
 * @endcode
 * 
 * Here's the first significant difference: We have to get the values
 * of the saturation function of the previous time step at the
 * quadrature points. To this end, we can use the
 * FEValues::get_function_values (previously already used in @ref step_9 "step-9",
 * @ref step_14 "step-14" and @ref step_15 "step-15"), a function that takes a solution vector and
 * returns a list of function values at the quadrature points of the
 * present cell. In fact, it returns the complete vector-valued
 * solution at each quadrature point, i.e. not only the saturation but
 * also the velocities and pressure:
 * 
 * @code
 *           fe_values.get_function_values(old_solution, old_solution_values);
 *   
 * @endcode
 * 
 * Then we also have to get the values of the pressure right hand side
 * and of the inverse permeability tensor at the quadrature points:
 * 
 * @code
 *           pressure_right_hand_side.value_list(fe_values.get_quadrature_points(),
 *                                               pressure_rhs_values);
 *           k_inverse.value_list(fe_values.get_quadrature_points(),
 *                                k_inverse_values);
 *   
 * @endcode
 * 
 * With all this, we can now loop over all the quadrature points and
 * shape functions on this cell and assemble those parts of the matrix
 * and right hand side that we deal with in this function. The
 * individual terms in the contributions should be self-explanatory
 * given the explicit form of the bilinear form stated in the
 * introduction:
 * 
 * @code
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *                 const double old_s = old_solution_values[q](dim + 1);
 *   
 *                 const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
 *                 const double div_phi_i_u = fe_values[velocities].divergence(i, q);
 *                 const double phi_i_p     = fe_values[pressure].value(i, q);
 *                 const double phi_i_s     = fe_values[saturation].value(i, q);
 *   
 *                 for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                   {
 *                     const Tensor<1, dim> phi_j_u =
 *                       fe_values[velocities].value(j, q);
 *                     const double div_phi_j_u =
 *                       fe_values[velocities].divergence(j, q);
 *                     const double phi_j_p = fe_values[pressure].value(j, q);
 *                     const double phi_j_s = fe_values[saturation].value(j, q);
 *   
 *                     local_matrix(i, j) +=
 *                       (phi_i_u * k_inverse_values[q] *
 *                          mobility_inverse(old_s, viscosity) * phi_j_u -
 *                        div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u +
 *                        phi_i_s * phi_j_s) *
 *                       fe_values.JxW(q);
 *                   }
 *   
 *                 local_rhs(i) +=
 *                   (-phi_i_p * pressure_rhs_values[q]) * fe_values.JxW(q);
 *               }
 *   
 *   
 * @endcode
 * 
 * Next, we also have to deal with the pressure boundary values. This,
 * again is as in @ref step_20 "step-20":
 * 
 * @code
 *           for (const auto &face : cell->face_iterators())
 *             if (face->at_boundary())
 *               {
 *                 fe_face_values.reinit(cell, face);
 *   
 *                 pressure_boundary_values.value_list(
 *                   fe_face_values.get_quadrature_points(), boundary_values);
 *   
 *                 for (unsigned int q = 0; q < n_face_q_points; ++q)
 *                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                     {
 *                       const Tensor<1, dim> phi_i_u =
 *                         fe_face_values[velocities].value(i, q);
 *   
 *                       local_rhs(i) +=
 *                         -(phi_i_u * fe_face_values.normal_vector(q) *
 *                           boundary_values[q] * fe_face_values.JxW(q));
 *                     }
 *               }
 *   
 * @endcode
 * 
 * The final step in the loop over all cells is to transfer local
 * contributions into the global matrix and right hand side vector:
 * 
 * @code
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *               system_matrix.add(local_dof_indices[i],
 *                                 local_dof_indices[j],
 *                                 local_matrix(i, j));
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             system_rhs(local_dof_indices[i]) += local_rhs(i);
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * So much for assembly of matrix and right hand side. Note that we do not
 * have to interpolate and apply boundary values since they have all been
 * taken care of in the weak form already.
 * 
 
 * 
 * 
 
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemassemble_rhs_S"></a> 
 * <h4>TwoPhaseFlowProblem::assemble_rhs_S</h4>
 * 
 
 * 
 * As explained in the introduction, we can only evaluate the right hand
 * side of the saturation equation once the velocity has been computed. We
 * therefore have this separate function to this end.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::assemble_rhs_S()
 *     {
 *       const QGauss<dim>     quadrature_formula(degree + 2);
 *       const QGauss<dim - 1> face_quadrature_formula(degree + 2);
 *       FEValues<dim>         fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       FEFaceValues<dim>     fe_face_values(fe,
 *                                        face_quadrature_formula,
 *                                        update_values | update_normal_vectors |
 *                                          update_quadrature_points |
 *                                          update_JxW_values);
 *       FEFaceValues<dim>     fe_face_values_neighbor(fe,
 *                                                 face_quadrature_formula,
 *                                                 update_values);
 *   
 *       const unsigned int dofs_per_cell   = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points      = quadrature_formula.size();
 *       const unsigned int n_face_q_points = face_quadrature_formula.size();
 *   
 *       Vector<double> local_rhs(dofs_per_cell);
 *   
 *       std::vector<Vector<double>> old_solution_values(n_q_points,
 *                                                       Vector<double>(dim + 2));
 *       std::vector<Vector<double>> old_solution_values_face(n_face_q_points,
 *                                                            Vector<double>(dim +
 *                                                                           2));
 *       std::vector<Vector<double>> old_solution_values_face_neighbor(
 *         n_face_q_points, Vector<double>(dim + 2));
 *       std::vector<Vector<double>> present_solution_values(n_q_points,
 *                                                           Vector<double>(dim +
 *                                                                          2));
 *       std::vector<Vector<double>> present_solution_values_face(
 *         n_face_q_points, Vector<double>(dim + 2));
 *   
 *       std::vector<double>                  neighbor_saturation(n_face_q_points);
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       SaturationBoundaryValues<dim> saturation_boundary_values;
 *   
 *       const FEValuesExtractors::Scalar saturation(dim + 1);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           local_rhs = 0;
 *           fe_values.reinit(cell);
 *   
 *           fe_values.get_function_values(old_solution, old_solution_values);
 *           fe_values.get_function_values(solution, present_solution_values);
 *   
 * @endcode
 * 
 * First for the cell terms. These are, following the formulas in the
 * introduction, @f$(S^n,\sigma)-(F(S^n) \mathbf{v}^{n+1},\nabla
 * \sigma)@f$, where @f$\sigma@f$ is the saturation component of the test
 * function:
 * 
 * @code
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *                 const double   old_s = old_solution_values[q](dim + 1);
 *                 Tensor<1, dim> present_u;
 *                 for (unsigned int d = 0; d < dim; ++d)
 *                   present_u[d] = present_solution_values[q](d);
 *   
 *                 const double         phi_i_s = fe_values[saturation].value(i, q);
 *                 const Tensor<1, dim> grad_phi_i_s =
 *                   fe_values[saturation].gradient(i, q);
 *   
 *                 local_rhs(i) +=
 *                   (time.get_next_step_size() * fractional_flow(old_s, viscosity) *
 *                      present_u * grad_phi_i_s +
 *                    old_s * phi_i_s) *
 *                   fe_values.JxW(q);
 *               }
 *   
 * @endcode
 * 
 * Secondly, we have to deal with the flux parts on the face
 * boundaries. This was a bit more involved because we first have to
 * determine which are the influx and outflux parts of the cell
 * boundary. If we have an influx boundary, we need to evaluate the
 * saturation on the other side of the face (or the boundary values,
 * if we are at the boundary of the domain).
 *         
 
 * 
 * All this is a bit tricky, but has been explained in some detail
 * already in @ref step_9 "step-9". Take a look there how this is supposed to work!
 * 
 * @code
 *           for (const auto face_no : cell->face_indices())
 *             {
 *               fe_face_values.reinit(cell, face_no);
 *   
 *               fe_face_values.get_function_values(old_solution,
 *                                                  old_solution_values_face);
 *               fe_face_values.get_function_values(solution,
 *                                                  present_solution_values_face);
 *   
 *               if (cell->at_boundary(face_no))
 *                 saturation_boundary_values.value_list(
 *                   fe_face_values.get_quadrature_points(), neighbor_saturation);
 *               else
 *                 {
 *                   const auto         neighbor = cell->neighbor(face_no);
 *                   const unsigned int neighbor_face =
 *                     cell->neighbor_of_neighbor(face_no);
 *   
 *                   fe_face_values_neighbor.reinit(neighbor, neighbor_face);
 *   
 *                   fe_face_values_neighbor.get_function_values(
 *                     old_solution, old_solution_values_face_neighbor);
 *   
 *                   for (unsigned int q = 0; q < n_face_q_points; ++q)
 *                     neighbor_saturation[q] =
 *                       old_solution_values_face_neighbor[q](dim + 1);
 *                 }
 *   
 *   
 *               for (unsigned int q = 0; q < n_face_q_points; ++q)
 *                 {
 *                   Tensor<1, dim> present_u_face;
 *                   for (unsigned int d = 0; d < dim; ++d)
 *                     present_u_face[d] = present_solution_values_face[q](d);
 *   
 *                   const double normal_flux =
 *                     present_u_face * fe_face_values.normal_vector(q);
 *   
 *                   const bool is_outflow_q_point = (normal_flux >= 0);
 *   
 *                   for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                     local_rhs(i) -=
 *                       time.get_next_step_size() * normal_flux *
 *                       fractional_flow((is_outflow_q_point == true ?
 *                                          old_solution_values_face[q](dim + 1) :
 *                                          neighbor_saturation[q]),
 *                                       viscosity) *
 *                       fe_face_values[saturation].value(i, q) *
 *                       fe_face_values.JxW(q);
 *                 }
 *             }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             system_rhs(local_dof_indices[i]) += local_rhs(i);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemsolve"></a> 
 * <h4>TwoPhaseFlowProblem::solve</h4>
 * 
 
 * 
 * After all these preparations, we finally solve the linear system for
 * velocity and pressure in the same way as in @ref step_20 "step-20". After that, we have
 * to deal with the saturation equation (see below):
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::solve()
 *     {
 *       const InverseMatrix<SparseMatrix<double>> m_inverse(
 *         system_matrix.block(0, 0));
 *       Vector<double> tmp(solution.block(0).size());
 *       Vector<double> schur_rhs(solution.block(1).size());
 *       Vector<double> tmp2(solution.block(2).size());
 *   
 *   
 * @endcode
 * 
 * First the pressure, using the pressure Schur complement of the first
 * two equations:
 * 
 * @code
 *       {
 *         m_inverse.vmult(tmp, system_rhs.block(0));
 *         system_matrix.block(1, 0).vmult(schur_rhs, tmp);
 *         schur_rhs -= system_rhs.block(1);
 *   
 *   
 *         SchurComplement schur_complement(system_matrix, m_inverse);
 *   
 *         ApproximateSchurComplement approximate_schur_complement(system_matrix);
 *   
 *         InverseMatrix<ApproximateSchurComplement> preconditioner(
 *           approximate_schur_complement);
 *   
 *   
 *         SolverControl            solver_control(solution.block(1).size(),
 *                                      1e-12 * schur_rhs.l2_norm());
 *         SolverCG<Vector<double>> cg(solver_control);
 *   
 *         cg.solve(schur_complement, solution.block(1), schur_rhs, preconditioner);
 *   
 *         std::cout << "   " << solver_control.last_step()
 *                   << " CG Schur complement iterations for pressure." << std::endl;
 *       }
 *   
 * @endcode
 * 
 * Now the velocity:
 * 
 * @code
 *       {
 *         system_matrix.block(0, 1).vmult(tmp, solution.block(1));
 *         tmp *= -1;
 *         tmp += system_rhs.block(0);
 *   
 *         m_inverse.vmult(solution.block(0), tmp);
 *       }
 *   
 * @endcode
 * 
 * Finally, we have to take care of the saturation equation. The first
 * business we have here is to determine the time step using the formula
 * in the introduction. Knowing the shape of our domain and that we
 * created the mesh by regular subdivision of cells, we can compute the
 * diameter of each of our cells quite easily (in fact we use the linear
 * extensions in coordinate directions of the cells, not the
 * diameter). Note that we will learn a more general way to do this in
 * @ref step_24 "step-24", where we use the GridTools::minimal_cell_diameter function.
 *     
 
 * 
 * The maximal velocity we compute using a helper function to compute the
 * maximal velocity defined below, and with all this we can evaluate our
 * new time step length. We use the method
 * DiscreteTime::set_desired_next_time_step() to suggest the new
 * calculated value of the time step to the DiscreteTime object. In most
 * cases, the time object uses the exact provided value to increment time.
 * It some case, the step size may be modified further by the time object.
 * For example, if the calculated time increment overshoots the end time,
 * it is truncated accordingly.
 * 
 * @code
 *       time.set_desired_next_step_size(std::pow(0.5, double(n_refinement_steps)) /
 *                                       get_maximal_velocity());
 *   
 * @endcode
 * 
 * The next step is to assemble the right hand side, and then to pass
 * everything on for solution. At the end, we project back saturations
 * onto the physically reasonable range:
 * 
 * @code
 *       assemble_rhs_S();
 *       {
 *         SolverControl            solver_control(system_matrix.block(2, 2).m(),
 *                                      1e-8 * system_rhs.block(2).l2_norm());
 *         SolverCG<Vector<double>> cg(solver_control);
 *         cg.solve(system_matrix.block(2, 2),
 *                  solution.block(2),
 *                  system_rhs.block(2),
 *                  PreconditionIdentity());
 *   
 *         project_back_saturation();
 *   
 *         std::cout << "   " << solver_control.last_step()
 *                   << " CG iterations for saturation." << std::endl;
 *       }
 *   
 *   
 *       old_solution = solution;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemoutput_results"></a> 
 * <h4>TwoPhaseFlowProblem::output_results</h4>
 * 
 
 * 
 * There is nothing surprising here. Since the program will do a lot of time
 * steps, we create an output file only every fifth time step and skip all
 * other time steps at the top of the file already.
 *   
 
 * 
 * When creating file names for output close to the bottom of the function,
 * we convert the number of the time step to a string representation that
 * is padded by leading zeros to four digits. We do this because this way
 * all output file names have the same length, and consequently sort well
 * when creating a directory listing.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::output_results() const
 *     {
 *       if (time.get_step_number() % 5 != 0)
 *         return;
 *   
 *       std::vector<std::string> solution_names;
 *       switch (dim)
 *         {
 *           case 2:
 *             solution_names = {"u", "v", "p", "S"};
 *             break;
 *   
 *           case 3:
 *             solution_names = {"u", "v", "w", "p", "S"};
 *             break;
 *   
 *           default:
 *             DEAL_II_NOT_IMPLEMENTED();
 *         }
 *   
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, solution_names);
 *   
 *       data_out.build_patches(degree + 1);
 *   
 *       std::ofstream output("solution-" +
 *                            Utilities::int_to_string(time.get_step_number(), 4) +
 *                            ".vtk");
 *       data_out.write_vtk(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemproject_back_saturation"></a> 
 * <h4>TwoPhaseFlowProblem::project_back_saturation</h4>
 * 
 
 * 
 * In this function, we simply run over all saturation degrees of freedom
 * and make sure that if they should have left the physically reasonable
 * range, that they be reset to the interval @f$[0,1]@f$. To do this, we only
 * have to loop over all saturation components of the solution vector; these
 * are stored in the block 2 (block 0 are the velocities, block 1 are the
 * pressures).
 *   
 
 * 
 * It may be instructive to note that this function almost never triggers
 * when the time step is chosen as mentioned in the introduction. However,
 * if we choose the timestep only slightly larger, we get plenty of values
 * outside the proper range. Strictly speaking, the function is therefore
 * unnecessary if we choose the time step small enough. In a sense, the
 * function is therefore only a safety device to avoid situations where our
 * entire solution becomes unphysical because individual degrees of freedom
 * have become unphysical a few time steps earlier.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::project_back_saturation()
 *     {
 *       for (unsigned int i = 0; i < solution.block(2).size(); ++i)
 *         if (solution.block(2)(i) < 0)
 *           solution.block(2)(i) = 0;
 *         else if (solution.block(2)(i) > 1)
 *           solution.block(2)(i) = 1;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemget_maximal_velocity"></a> 
 * <h4>TwoPhaseFlowProblem::get_maximal_velocity</h4>
 * 
 
 * 
 * The following function is used in determining the maximal allowable time
 * step. What it does is to loop over all quadrature points in the domain
 * and find what the maximal magnitude of the velocity is.
 * 
 * @code
 *     template <int dim>
 *     double TwoPhaseFlowProblem<dim>::get_maximal_velocity() const
 *     {
 *       const QGauss<dim>  quadrature_formula(degree + 2);
 *       const unsigned int n_q_points = quadrature_formula.size();
 *   
 *       FEValues<dim> fe_values(fe, quadrature_formula, update_values);
 *       std::vector<Vector<double>> solution_values(n_q_points,
 *                                                   Vector<double>(dim + 2));
 *       double                      max_velocity = 0;
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           fe_values.get_function_values(solution, solution_values);
 *   
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               Tensor<1, dim> velocity;
 *               for (unsigned int i = 0; i < dim; ++i)
 *                 velocity[i] = solution_values[q](i);
 *   
 *               max_velocity = std::max(max_velocity, velocity.norm());
 *             }
 *         }
 *   
 *       return max_velocity;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-TwoPhaseFlowProblemrun"></a> 
 * <h4>TwoPhaseFlowProblem::run</h4>
 * 
 
 * 
 * This is the final function of our main class. Its brevity speaks for
 * itself. There are only two points worth noting: First, the function
 * projects the initial values onto the finite element space at the
 * beginning; the VectorTools::project function doing this requires an
 * argument indicating the hanging node constraints. We have none in this
 * program (we compute on a uniformly refined mesh), but the function
 * requires the argument anyway, of course. So we have to create a
 * constraint object. In its original state, constraint objects are
 * unsorted, and have to be sorted (using the AffineConstraints::close
 * function) before they can be used. This is what we do here, and which is
 * why we can't simply call the VectorTools::project function with an
 * anonymous temporary object <code>AffineConstraints<double>()</code> as the
 * second argument.
 *   
 
 * 
 * The second point worth mentioning is that we only compute the length of
 * the present time step in the middle of solving the linear system
 * corresponding to each time step. We can therefore output the present
 * time of a time step only at the end of the time step.
 * We increment time by calling the method DiscreteTime::advance_time()
 * inside the loop. Since we are reporting the time and dt after we
 * increment it, we have to call the method
 * DiscreteTime::get_previous_step_size() instead of
 * DiscreteTime::get_next_step_size(). After many steps, when the simulation
 * reaches the end time, the last dt is chosen by the DiscreteTime class in
 * such a way that the last step finishes exactly at the end time.
 * 
 * @code
 *     template <int dim>
 *     void TwoPhaseFlowProblem<dim>::run()
 *     {
 *       make_grid_and_dofs();
 *   
 *       {
 *         AffineConstraints<double> constraints;
 *         constraints.close();
 *   
 *         VectorTools::project(dof_handler,
 *                              constraints,
 *                              QGauss<dim>(degree + 2),
 *                              InitialValues<dim>(),
 *                              old_solution);
 *       }
 *   
 *       do
 *         {
 *           std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
 *   
 *           assemble_system();
 *   
 *           solve();
 *   
 *           output_results();
 *   
 *           time.advance_time();
 *           std::cout << "   Now at t=" << time.get_current_time()
 *                     << ", dt=" << time.get_previous_step_size() << '.'
 *                     << std::endl
 *                     << std::endl;
 *         }
 *       while (time.is_at_end() == false);
 *     }
 *   } // namespace Step21
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_21-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * That's it. In the main function, we pass the degree of the finite element
 * space to the constructor of the TwoPhaseFlowProblem object.  Here, we use
 * zero-th degree elements, i.e. @f$RT_0\times DQ_0 \times DQ_0@f$. The rest is as
 * in all the other programs.
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step21;
 *   
 *         TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
 *         two_phase_flow_problem.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_21-Results"></a><h1>Results</h1>
 
 
The code as presented here does not actually compute the results
found on the web page. The reason is, that even on a decent
computer it runs more than a day. If you want to reproduce these
results, modify the end time of the DiscreteTime object to `250` within the
constructor of TwoPhaseFlowProblem.
 
If we run the program, we get the following kind of output:
@code
Number of active cells: 1024
Number of degrees of freedom: 4160 (2112+1024+1024)
 
Timestep 1
   22 CG Schur complement iterations for pressure.
   1 CG iterations for saturation.
   Now at t=0.0326742, dt=0.0326742.
 
Timestep 2
   17 CG Schur complement iterations for pressure.
   1 CG iterations for saturation.
   Now at t=0.0653816, dt=0.0327074.
 
Timestep 3
   17 CG Schur complement iterations for pressure.
   1 CG iterations for saturation.
   Now at t=0.0980651, dt=0.0326836.
 
...
@endcode
As we can see, the time step is pretty much constant right from the start,
which indicates that the velocities in the domain are not strongly dependent
on changes in saturation, although they certainly are through the factor
@f$\lambda(S)@f$ in the pressure equation.
 
Our second observation is that the number of CG iterations needed to solve the
pressure Schur complement equation drops from 22 to 17 between the first and
the second time step (in fact, it remains around 17 for the rest of the
computations). The reason is actually simple: Before we solve for the pressure
during a time step, we don't reset the <code>solution</code> variable to
zero. The pressure (and the other variables) therefore have the previous time
step's values at the time we get into the CG solver. Since the velocities and
pressures don't change very much as computations progress, the previous time
step's pressure is actually a good initial guess for this time step's
pressure. Consequently, the number of iterations we need once we have computed
the pressure once is significantly reduced.
 
The final observation concerns the number of iterations needed to solve for
the saturation, i.e. one. This shouldn't surprise us too much: the matrix we
have to solve with is the mass matrix. However, this is the mass matrix for
the @f$DGQ_0@f$ element of piecewise constants where no element couples with the
degrees of freedom on neighboring cells. The matrix is therefore a diagonal
one, and it is clear that we should be able to invert this matrix in a single
CG iteration.
 
 
With all this, here are a few movies that show how the saturation progresses
over time. First, this is for the single crack model, as implemented in the
<code>SingleCurvingCrack::KInverse</code> class:
 
<img src="https://www.dealii.org/images/steps/developer/step-21.centerline.gif" alt="">
 
As can be seen, the water rich fluid snakes its way mostly along the
high-permeability zone in the middle of the domain, whereas the rest of the
domain is mostly impermeable. This and the next movie are generated using
<code>n_refinement_steps=7</code>, leading to a @f$128\times 128@f$ mesh with some
16,000 cells and about 66,000 unknowns in total.
 
 
The second movie shows the saturation for the random medium model of class
<code>RandomMedium::KInverse</code>, where we have randomly distributed
centers of high permeability and fluid hops from one of these zones to
the next:
 
<img src="https://www.dealii.org/images/steps/developer/step-21.random2d.gif" alt="">
 
 
Finally, here is the same situation in three space dimensions, on a mesh with
<code>n_refinement_steps=5</code>, which produces a mesh of some 32,000 cells
and 167,000 degrees of freedom:
 
<img src="https://www.dealii.org/images/steps/developer/step-21.random3d.gif" alt="">
 
To repeat these computations, all you have to do is to change the line
@code
      TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
@endcode
in the main function to
@code
      TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
@endcode
The visualization uses a cloud technique, where the saturation is indicated by
colored but transparent clouds for each cell. This way, one can also see
somewhat what happens deep inside the domain. A different way of visualizing
would have been to show isosurfaces of the saturation evolving over
time. There are techniques to plot isosurfaces transparently, so that one can
see several of them at the same time like the layers of an onion.
 
So why don't we show such isosurfaces? The problem lies in the way isosurfaces
are computed: they require that the field to be visualized is continuous, so
that the isosurfaces can be generated by following contours at least across a
single cell. However, our saturation field is piecewise constant and
discontinuous. If we wanted to plot an isosurface for a saturation @f$S=0.5@f$,
chances would be that there is no single point in the domain where that
saturation is actually attained. If we had to define isosurfaces in that
context at all, we would have to take the interfaces between cells, where one
of the two adjacent cells has a saturation greater than and the other cell a
saturation less than 0.5. However, it appears that most visualization programs
are not equipped to do this kind of transformation.
 
 
<a name="step-21-extensions"></a>
<a name="step_21-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
There are a number of areas where this program can be improved. Three of them
are listed below. All of them are, in fact, addressed in a tutorial program
that forms the continuation of the current one: @ref step_43 "step-43".
 
 
<a name="step_21-Solvers"></a><h4>Solvers</h4>
 
 
At present, the program is not particularly fast: the 2d random medium
computation took about a day for the 1,000 or so time steps. The corresponding
3d computation took almost two days for 800 time steps. The reason why it
isn't faster than this is twofold. First, we rebuild the entire matrix in
every time step, although some parts such as the @f$B@f$, @f$B^T@f$, and @f$M^S@f$ blocks
never change.
 
Second, we could do a lot better with the solver and
preconditioners. Presently, we solve the Schur complement @f$B^TM^u(S)^{-1}B@f$
with a CG method, using @f$[B^T (\textrm{diag}(M^u(S)))^{-1} B]^{-1}@f$ as a
preconditioner. Applying this preconditioner is expensive, since it involves
solving a linear system each time. This may have been appropriate for @ref
step_20 "step-20", where we have to solve the entire problem only
once. However, here we have to solve it hundreds of times, and in such cases
it is worth considering a preconditioner that is more expensive to set up the
first time, but cheaper to apply later on.
 
One possibility would be to realize that the matrix we use as preconditioner,
@f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ is still sparse, and symmetric on top of
that. If one looks at the flow field evolve over time, we also see that while
@f$S@f$ changes significantly over time, the pressure hardly does and consequently
@f$B^T (\textrm{diag}(M^u(S)))^{-1} B \approx B^T (\textrm{diag}(M^u(S^0)))^{-1}
B@f$. In other words, the matrix for the first time step should be a good
preconditioner also for all later time steps.  With a bit of
back-and-forthing, it isn't hard to actually get a representation of it as a
SparseMatrix object. We could then hand it off to the SparseMIC class to form
a sparse incomplete Cholesky decomposition. To form this decomposition is
expensive, but we have to do it only once in the first time step, and can then
use it as a cheap preconditioner in the future. We could do better even by
using the SparseDirectUMFPACK class that produces not only an incomplete, but
a complete decomposition of the matrix, which should yield an even better
preconditioner.
 
Finally, why use the approximation @f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ to
precondition @f$B^T M^u(S)^{-1} B@f$? The latter matrix, after all, is the mixed
form of the Laplace operator on the pressure space, for which we use linear
elements. We could therefore build a separate matrix @f$A^p@f$ on the side that
directly corresponds to the non-mixed formulation of the Laplacian, for
example using the bilinear form @f$(\mathbf{K}\lambda(S^n) \nabla
\varphi_i,\nabla\varphi_j)@f$. We could then form an incomplete or complete
decomposition of this non-mixed matrix and use it as a preconditioner of the
mixed form.
 
Using such techniques, it can reasonably be expected that the solution process
will be faster by at least an order of magnitude.
 
 
<a name="step_21-Timestepping"></a><h4>Time stepping</h4>
 
 
In the introduction we have identified the time step restriction
@f[
  \triangle t_{n+1} \le \frac h{|\mathbf{u}^{n+1}(\mathbf{x})|}
@f]
that has to hold globally, i.e. for all @f$\mathbf x@f$. After discretization, we
satisfy it by choosing
@f[
  \triangle t_{n+1} = \frac {\min_K h_K}{\max_{\mathbf{x}}|\mathbf{u}^{n+1}(\mathbf{x})|}.
@f]
 
This restriction on the time step is somewhat annoying: the finer we make the
mesh the smaller the time step; in other words, we get punished twice: each
time step is more expensive to solve and we have to do more time steps.
 
This is particularly annoying since the majority of the additional work is
spent solving the implicit part of the equations, i.e. the pressure-velocity
system, whereas it is the hyperbolic transport equation for the saturation
that imposes the time step restriction.
 
To avoid this bottleneck, people have invented a number of approaches. For
example, they may only re-compute the pressure-velocity field every few time
steps (or, if you want, use different time step sizes for the
pressure/velocity and saturation equations). This keeps the time step
restriction on the cheap explicit part while it makes the solution of the
implicit part less frequent. Experiments in this direction are
certainly worthwhile; one starting point for such an approach is the paper by
Zhangxin Chen, Guanren Huan and Baoyan Li: <i>An improved IMPES method for
two-phase flow in porous media</i>, Transport in Porous Media, 54 (2004),
pp. 361&mdash;376. There are certainly many other papers on this topic as well, but
this one happened to land on our desk a while back.
 
 
 
<a name="step_21-Adaptivity"></a><h4>Adaptivity</h4>
 
 
Adaptivity would also clearly help. Looking at the movies, one clearly sees
that most of the action is confined to a relatively small part of the domain
(this particularly obvious for the saturation, but also holds for the
velocities and pressures). Adaptivity can therefore be expected to keep the
necessary number of degrees of freedom low, or alternatively increase the
accuracy.
 
On the other hand, adaptivity for time dependent problems is not a trivial
thing: we would have to change the mesh every few time steps, and we would
have to transport our present solution to the next mesh every time we change
it (something that the SolutionTransfer class can help with). These are not
insurmountable obstacles, but they do require some additional coding and more
than we felt comfortable was worth packing into this tutorial program.
 *
 *
<a name="step_21-PlainProg"></a>
<h1> The plain program</h1>
@include "step-21.cc"
*/