deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-43.h
Go to the documentation of this file.
1) const
741 *   {
742 *   return 1 - p[0];
743 *   }
744 *  
745 *  
746 *   template <int dim>
747 *   class SaturationBoundaryValues : public Function<dim>
748 *   {
749 *   public:
750 *   SaturationBoundaryValues()
751 *   : Function<dim>(1)
752 *   {}
753 *  
754 *   virtual double value(const Point<dim> &p,
755 *   const unsigned int component = 0) const override;
756 *   };
757 *  
758 *  
759 *  
760 *   template <int dim>
761 *   double
762 *   SaturationBoundaryValues<dim>::value(const Point<dim> &p,
763 *   const unsigned int /*component*/) const
764 *   {
765 *   if (p[0] == 0)
766 *   return 1;
767 *   else
768 *   return 0;
769 *   }
770 *  
771 *  
772 *   template <int dim>
773 *   class SaturationInitialValues : public Function<dim>
774 *   {
775 *   public:
776 *   SaturationInitialValues()
777 *   : Function<dim>(1)
778 *   {}
779 *  
780 *   virtual double value(const Point<dim> &p,
781 *   const unsigned int component = 0) const override;
782 *  
783 *   virtual void vector_value(const Point<dim> &p,
784 *   Vector<double> &value) const override;
785 *   };
786 *  
787 *  
788 *   template <int dim>
789 *   double
790 *   SaturationInitialValues<dim>::value(const Point<dim> & /*p*/,
791 *   const unsigned int /*component*/) const
792 *   {
793 *   return 0.2;
794 *   }
795 *  
796 *  
797 *   template <int dim>
798 *   void SaturationInitialValues<dim>::vector_value(const Point<dim> &p,
799 *   Vector<double> &values) const
800 *   {
801 *   for (unsigned int c = 0; c < this->n_components; ++c)
802 *   values(c) = SaturationInitialValues<dim>::value(p, c);
803 *   }
804 *  
805 *  
806 * @endcode
807 *
808 *
809 * <a name="step_43-Permeabilitymodels"></a>
810 * <h3>Permeability models</h3>
811 *
812
813 *
814 * In this tutorial, we still use the two permeability models previously
815 * used in @ref step_21 "step-21" so we again refrain from commenting in detail about them.
816 *
817 * @code
818 *   namespace SingleCurvingCrack
819 *   {
820 *   template <int dim>
821 *   class KInverse : public TensorFunction<2, dim>
822 *   {
823 *   public:
824 *   KInverse()
825 *   : TensorFunction<2, dim>()
826 *   {}
827 *  
828 *   virtual void
829 *   value_list(const std::vector<Point<dim>> &points,
830 *   std::vector<Tensor<2, dim>> &values) const override;
831 *   };
832 *  
833 *  
834 *   template <int dim>
835 *   void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
836 *   std::vector<Tensor<2, dim>> &values) const
837 *   {
838 *   AssertDimension(points.size(), values.size());
839 *  
840 *   for (unsigned int p = 0; p < points.size(); ++p)
841 *   {
842 *   values[p].clear();
843 *  
844 *   const double distance_to_flowline =
845 *   std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
846 *  
847 *   const double permeability =
848 *   std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
849 *   (0.1 * 0.1)),
850 *   0.01);
851 *  
852 *   for (unsigned int d = 0; d < dim; ++d)
853 *   values[p][d][d] = 1. / permeability;
854 *   }
855 *   }
856 *   } // namespace SingleCurvingCrack
857 *  
858 *  
859 *   namespace RandomMedium
860 *   {
861 *   template <int dim>
862 *   class KInverse : public TensorFunction<2, dim>
863 *   {
864 *   public:
865 *   KInverse()
866 *   : TensorFunction<2, dim>()
867 *   {}
868 *  
869 *   virtual void
870 *   value_list(const std::vector<Point<dim>> &points,
871 *   std::vector<Tensor<2, dim>> &values) const override;
872 *  
873 *   private:
874 *   static std::vector<Point<dim>> centers;
875 *   };
876 *  
877 *  
878 *  
879 *   template <int dim>
880 *   std::vector<Point<dim>> KInverse<dim>::centers = []() {
881 *   const unsigned int N =
882 *   (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
883 *  
884 *   std::vector<Point<dim>> centers_list(N);
885 *   for (unsigned int i = 0; i < N; ++i)
886 *   for (unsigned int d = 0; d < dim; ++d)
887 *   centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
888 *  
889 *   return centers_list;
890 *   }();
891 *  
892 *  
893 *  
894 *   template <int dim>
895 *   void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
896 *   std::vector<Tensor<2, dim>> &values) const
897 *   {
898 *   AssertDimension(points.size(), values.size());
899 *  
900 *   for (unsigned int p = 0; p < points.size(); ++p)
901 *   {
902 *   values[p].clear();
903 *  
904 *   double permeability = 0;
905 *   for (unsigned int i = 0; i < centers.size(); ++i)
906 *   permeability +=
907 *   std::exp(-(points[p] - centers[i]).norm_square() / (0.05 * 0.05));
908 *  
909 *   const double normalized_permeability =
910 *   std::clamp(permeability, 0.01, 4.);
911 *  
912 *   for (unsigned int d = 0; d < dim; ++d)
913 *   values[p][d][d] = 1. / normalized_permeability;
914 *   }
915 *   }
916 *   } // namespace RandomMedium
917 *  
918 *  
919 * @endcode
920 *
921 *
922 * <a name="step_43-Physicalquantities"></a>
923 * <h3>Physical quantities</h3>
924 *
925
926 *
927 * The implementations of all the physical quantities such as total mobility
928 * @f$\lambda_t@f$ and fractional flow of water @f$F@f$ are taken from @ref step_21 "step-21" so
929 * again we don't have do any comment about them. Compared to @ref step_21 "step-21" we
930 * have added checks that the saturation passed to these functions is in
931 * fact within the physically valid range. Furthermore, given that the
932 * wetting phase moves at speed @f$\mathbf u F'(S)@f$ it is clear that @f$F'(S)@f$
933 * must be greater or equal to zero, so we assert that as well to make sure
934 * that our calculations to get at the formula for the derivative made
935 * sense.
936 *
937 * @code
938 *   double mobility_inverse(const double S, const double viscosity)
939 *   {
940 *   return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
941 *   }
942 *  
943 *  
944 *   double fractional_flow(const double S, const double viscosity)
945 *   {
946 *   Assert((S >= 0) && (S <= 1),
947 *   ExcMessage("Saturation is outside its physically valid range."));
948 *  
949 *   return S * S / (S * S + viscosity * (1 - S) * (1 - S));
950 *   }
951 *  
952 *  
953 *   double fractional_flow_derivative(const double S, const double viscosity)
954 *   {
955 *   Assert((S >= 0) && (S <= 1),
956 *   ExcMessage("Saturation is outside its physically valid range."));
957 *  
958 *   const double temp = (S * S + viscosity * (1 - S) * (1 - S));
959 *  
960 *   const double numerator =
961 *   2.0 * S * temp - S * S * (2.0 * S - 2.0 * viscosity * (1 - S));
962 *   const double denominator = Utilities::fixed_power<2>(temp);
963 *  
964 *   const double F_prime = numerator / denominator;
965 *  
966 *   Assert(F_prime >= 0, ExcInternalError());
967 *  
968 *   return F_prime;
969 *   }
970 *  
971 *  
972 * @endcode
973 *
974 *
975 * <a name="step_43-Helperclassesforsolversandpreconditioners"></a>
976 * <h3>Helper classes for solvers and preconditioners</h3>
977 *
978
979 *
980 * In this first part we define a number of classes that we need in the
981 * construction of linear solvers and preconditioners. This part is
982 * essentially the same as that used in @ref step_31 "step-31". The only difference is that
983 * the original variable name stokes_matrix is replaced by another name
984 * darcy_matrix to match our problem.
985 *
986 * @code
987 *   namespace LinearSolvers
988 *   {
989 *   template <class MatrixType, class PreconditionerType>
990 *   class InverseMatrix : public EnableObserverPointer
991 *   {
992 *   public:
993 *   InverseMatrix(const MatrixType &m,
994 *   const PreconditionerType &preconditioner);
995 *  
996 *  
997 *   template <typename VectorType>
998 *   void vmult(VectorType &dst, const VectorType &src) const;
999 *  
1000 *   private:
1001 *   const ObserverPointer<const MatrixType> matrix;
1002 *   const PreconditionerType &preconditioner;
1003 *   };
1004 *  
1005 *  
1006 *   template <class MatrixType, class PreconditionerType>
1007 *   InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
1008 *   const MatrixType &m,
1009 *   const PreconditionerType &preconditioner)
1010 *   : matrix(&m)
1011 *   , preconditioner(preconditioner)
1012 *   {}
1013 *  
1014 *  
1015 *  
1016 *   template <class MatrixType, class PreconditionerType>
1017 *   template <typename VectorType>
1018 *   void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1019 *   VectorType &dst,
1020 *   const VectorType &src) const
1021 *   {
1022 *   SolverControl solver_control(src.size(), 1e-7 * src.l2_norm());
1023 *   SolverCG<VectorType> cg(solver_control);
1024 *  
1025 *   dst = 0;
1026 *  
1027 *   try
1028 *   {
1029 *   cg.solve(*matrix, dst, src, preconditioner);
1030 *   }
1031 *   catch (std::exception &e)
1032 *   {
1033 *   Assert(false, ExcMessage(e.what()));
1034 *   }
1035 *   }
1036 *  
1037 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1038 *   class BlockSchurPreconditioner : public EnableObserverPointer
1039 *   {
1040 *   public:
1041 *   BlockSchurPreconditioner(
1042 *   const TrilinosWrappers::BlockSparseMatrix &S,
1043 *   const InverseMatrix<TrilinosWrappers::SparseMatrix,
1044 *   PreconditionerTypeMp> &Mpinv,
1045 *   const PreconditionerTypeA &Apreconditioner);
1046 *  
1047 *   void vmult(TrilinosWrappers::MPI::BlockVector &dst,
1048 *   const TrilinosWrappers::MPI::BlockVector &src) const;
1049 *  
1050 *   private:
1051 *   const ObserverPointer<const TrilinosWrappers::BlockSparseMatrix>
1052 *   darcy_matrix;
1053 *   const ObserverPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
1054 *   PreconditionerTypeMp>>
1055 *   m_inverse;
1056 *   const PreconditionerTypeA &a_preconditioner;
1057 *  
1058 *   mutable TrilinosWrappers::MPI::Vector tmp;
1059 *   };
1060 *  
1061 *  
1062 *  
1063 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1064 *   BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::
1065 *   BlockSchurPreconditioner(
1066 *   const TrilinosWrappers::BlockSparseMatrix &S,
1067 *   const InverseMatrix<TrilinosWrappers::SparseMatrix,
1068 *   PreconditionerTypeMp> &Mpinv,
1069 *   const PreconditionerTypeA &Apreconditioner)
1070 *   : darcy_matrix(&S)
1071 *   , m_inverse(&Mpinv)
1072 *   , a_preconditioner(Apreconditioner)
1073 *   , tmp(complete_index_set(darcy_matrix->block(1, 1).m()))
1074 *   {}
1075 *  
1076 *  
1077 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1078 *   void
1079 *   BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::vmult(
1080 *   TrilinosWrappers::MPI::BlockVector &dst,
1081 *   const TrilinosWrappers::MPI::BlockVector &src) const
1082 *   {
1083 *   a_preconditioner.vmult(dst.block(0), src.block(0));
1084 *   darcy_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
1085 *   tmp *= -1;
1086 *   m_inverse->vmult(dst.block(1), tmp);
1087 *   }
1088 *   } // namespace LinearSolvers
1089 *  
1090 *  
1091 * @endcode
1092 *
1093 *
1094 * <a name="step_43-TheTwoPhaseFlowProblemclass"></a>
1095 * <h3>The TwoPhaseFlowProblem class</h3>
1096 *
1097
1098 *
1099 * The definition of the class that defines the top-level logic of solving
1100 * the time-dependent advection-dominated two-phase flow problem (or
1101 * Buckley-Leverett problem @cite Buckley1942) is mainly based on tutorial
1102 * programs @ref step_21 "step-21" and @ref step_33 "step-33", and in particular on @ref step_31 "step-31" where we have
1103 * used basically the same general structure as done here. As in @ref step_31 "step-31",
1104 * the key routines to look for in the implementation below are the
1105 * <code>run()</code> and <code>solve()</code> functions.
1106 *
1107
1108 *
1109 * The main difference to @ref step_31 "step-31" is that, since adaptive operator splitting
1110 * is considered, we need a couple more member variables to hold the last
1111 * two computed Darcy (velocity/pressure) solutions in addition to the
1112 * current one (which is either computed directly, or extrapolated from the
1113 * previous two), and we need to remember the last two times we computed the
1114 * Darcy solution. We also need a helper function that figures out whether
1115 * we do indeed need to recompute the Darcy solution.
1116 *
1117
1118 *
1119 * Unlike @ref step_31 "step-31", this step uses one more AffineConstraints object called
1120 * darcy_preconditioner_constraints. This constraint object is used only for
1121 * assembling the matrix for the Darcy preconditioner and includes hanging
1122 * node constraints as well as Dirichlet boundary value constraints for the
1123 * pressure variable. We need this because we are building a Laplace matrix
1124 * for the pressure as an approximation of the Schur complement) which is
1125 * only positive definite if boundary conditions are applied.
1126 *
1127
1128 *
1129 * The collection of member functions and variables thus declared in this
1130 * class is then rather similar to those in @ref step_31 "step-31":
1131 *
1132 * @code
1133 *   template <int dim>
1134 *   class TwoPhaseFlowProblem
1135 *   {
1136 *   public:
1137 *   TwoPhaseFlowProblem(const unsigned int degree);
1138 *   void run();
1139 *  
1140 *   private:
1141 *   void setup_dofs();
1142 *   void assemble_darcy_preconditioner();
1143 *   void build_darcy_preconditioner();
1144 *   void assemble_darcy_system();
1145 *   void assemble_saturation_system();
1146 *   void assemble_saturation_matrix();
1147 *   void assemble_saturation_rhs();
1148 *   void assemble_saturation_rhs_cell_term(
1149 *   const FEValues<dim> &saturation_fe_values,
1150 *   const FEValues<dim> &darcy_fe_values,
1151 *   const double global_max_u_F_prime,
1152 *   const double global_S_variation,
1153 *   const std::vector<types::global_dof_index> &local_dof_indices);
1154 *   void assemble_saturation_rhs_boundary_term(
1155 *   const FEFaceValues<dim> &saturation_fe_face_values,
1156 *   const FEFaceValues<dim> &darcy_fe_face_values,
1157 *   const std::vector<types::global_dof_index> &local_dof_indices);
1158 *   void solve();
1159 *   void refine_mesh(const unsigned int min_grid_level,
1160 *   const unsigned int max_grid_level);
1161 *   void output_results() const;
1162 *  
1163 * @endcode
1164 *
1165 * We follow with a number of helper functions that are used in a variety
1166 * of places throughout the program:
1167 *
1168 * @code
1169 *   double get_max_u_F_prime() const;
1170 *   std::pair<double, double> get_extrapolated_saturation_range() const;
1171 *   bool determine_whether_to_solve_for_pressure_and_velocity() const;
1172 *   void project_back_saturation();
1173 *   double compute_viscosity(
1174 *   const std::vector<double> &old_saturation,
1175 *   const std::vector<double> &old_old_saturation,
1176 *   const std::vector<Tensor<1, dim>> &old_saturation_grads,
1177 *   const std::vector<Tensor<1, dim>> &old_old_saturation_grads,
1178 *   const std::vector<Vector<double>> &present_darcy_values,
1179 *   const double global_max_u_F_prime,
1180 *   const double global_S_variation,
1181 *   const double cell_diameter) const;
1182 *  
1183 *  
1184 * @endcode
1185 *
1186 * This all is followed by the member variables, most of which are similar
1187 * to the ones in @ref step_31 "step-31", with the exception of the ones that pertain to
1188 * the macro time stepping for the velocity/pressure system:
1189 *
1190 * @code
1191 *   Triangulation<dim> triangulation;
1192 *   double global_Omega_diameter;
1193 *  
1194 *   const unsigned int degree;
1195 *  
1196 *   const unsigned int darcy_degree;
1197 *   const FESystem<dim> darcy_fe;
1198 *   DoFHandler<dim> darcy_dof_handler;
1199 *   AffineConstraints<double> darcy_constraints;
1200 *  
1201 *   AffineConstraints<double> darcy_preconditioner_constraints;
1202 *  
1203 *   TrilinosWrappers::BlockSparseMatrix darcy_matrix;
1204 *   TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix;
1205 *  
1206 *   TrilinosWrappers::MPI::BlockVector darcy_solution;
1207 *   TrilinosWrappers::MPI::BlockVector darcy_rhs;
1208 *  
1209 *   TrilinosWrappers::MPI::BlockVector last_computed_darcy_solution;
1210 *   TrilinosWrappers::MPI::BlockVector second_last_computed_darcy_solution;
1211 *  
1212 *  
1213 *   const unsigned int saturation_degree;
1214 *   const FE_Q<dim> saturation_fe;
1215 *   DoFHandler<dim> saturation_dof_handler;
1216 *   AffineConstraints<double> saturation_constraints;
1217 *  
1218 *   TrilinosWrappers::SparseMatrix saturation_matrix;
1219 *  
1220 *  
1221 *   TrilinosWrappers::MPI::Vector saturation_solution;
1222 *   TrilinosWrappers::MPI::Vector old_saturation_solution;
1223 *   TrilinosWrappers::MPI::Vector old_old_saturation_solution;
1224 *   TrilinosWrappers::MPI::Vector saturation_rhs;
1225 *  
1226 *   TrilinosWrappers::MPI::Vector
1227 *   saturation_matching_last_computed_darcy_solution;
1228 *  
1229 *   const double saturation_refinement_threshold;
1230 *  
1231 *   double time;
1232 *   const double end_time;
1233 *  
1234 *   double current_macro_time_step;
1235 *   double old_macro_time_step;
1236 *  
1237 *   double time_step;
1238 *   double old_time_step;
1239 *   unsigned int timestep_number;
1240 *  
1241 *   const double viscosity;
1242 *   const double porosity;
1243 *   const double AOS_threshold;
1244 *  
1245 *   #ifdef DEAL_II_TRILINOS_WITH_EPETRA
1246 *   using PreconditionType = TrilinosWrappers::PreconditionIC;
1247 *   #else
1248 * @endcode
1249 *
1250 * For Tpetra, IC is only available through Ifpack.
1251 *
1252 * @code
1253 *   using PreconditionType =
1254 *   LinearAlgebra::TpetraWrappers::PreconditionIfpack<double>;
1255 *   #endif
1256 *   std::shared_ptr<PreconditionType> top_left_preconditioner;
1257 *   std::shared_ptr<PreconditionType> bottom_right_preconditioner;
1258 *  
1259 *   bool rebuild_saturation_matrix;
1260 *  
1261 * @endcode
1262 *
1263 * At the very end we declare a variable that denotes the material
1264 * model. Compared to @ref step_21 "step-21", we do this here as a member variable since
1265 * we will want to use it in a variety of places and so having a central
1266 * place where such a variable is declared will make it simpler to replace
1267 * one class by another (e.g. replace RandomMedium::KInverse by
1268 * SingleCurvingCrack::KInverse).
1269 *
1270 * @code
1271 *   const RandomMedium::KInverse<dim> k_inverse;
1272 *   };
1273 *  
1274 *  
1275 * @endcode
1276 *
1277 *
1278 * <a name="step_43-TwoPhaseFlowProblemdimTwoPhaseFlowProblem"></a>
1279 * <h3>TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem</h3>
1280 *
1281
1282 *
1283 * The constructor of this class is an extension of the constructors in
1284 * @ref step_21 "step-21" and @ref step_31 "step-31". We need to add the various variables that concern
1285 * the saturation. As discussed in the introduction, we are going to use
1286 * @f$Q_2 \times Q_1@f$ (Taylor-Hood) elements again for the Darcy system, an
1287 * element combination that fulfills the Ladyzhenskaya-Babuska-Brezzi (LBB)
1288 * conditions [Brezzi and Fortin 1991, Chen 2005], and @f$Q_1@f$ elements for
1289 * the saturation. However, by using variables that store the polynomial
1290 * degree of the Darcy and temperature finite elements, it is easy to
1291 * consistently modify the degree of the elements as well as all quadrature
1292 * formulas used on them downstream. Moreover, we initialize the time
1293 * stepping variables related to operator splitting as well as the option
1294 * for matrix assembly and preconditioning:
1295 *
1296 * @code
1297 *   template <int dim>
1298 *   TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem(const unsigned int degree)
1299 *   : triangulation(Triangulation<dim>::maximum_smoothing)
1300 *   , global_Omega_diameter(std::numeric_limits<double>::quiet_NaN())
1301 *   , degree(degree)
1302 *   , darcy_degree(degree)
1303 *   , darcy_fe(FE_Q<dim>(darcy_degree + 1) ^ dim, FE_Q<dim>(darcy_degree))
1304 *   , darcy_dof_handler(triangulation)
1305 *   ,
1306 *  
1307 *   saturation_degree(degree + 1)
1308 *   , saturation_fe(saturation_degree)
1309 *   , saturation_dof_handler(triangulation)
1310 *   ,
1311 *  
1312 *   saturation_refinement_threshold(0.5)
1313 *   ,
1314 *  
1315 *   time(0)
1316 *   , end_time(10)
1317 *   ,
1318 *  
1319 *   current_macro_time_step(0)
1320 *   , old_macro_time_step(0)
1321 *   ,
1322 *  
1323 *   time_step(0)
1324 *   , old_time_step(0)
1325 *   , timestep_number(0)
1326 *   , viscosity(0.2)
1327 *   , porosity(1.0)
1328 *   , AOS_threshold(3.0)
1329 *   ,
1330 *  
1331 *   rebuild_saturation_matrix(true)
1332 *   {}
1333 *  
1334 *  
1335 * @endcode
1336 *
1337 *
1338 * <a name="step_43-TwoPhaseFlowProblemdimsetup_dofs"></a>
1339 * <h3>TwoPhaseFlowProblem<dim>::setup_dofs</h3>
1340 *
1341
1342 *
1343 * This is the function that sets up the DoFHandler objects we have here
1344 * (one for the Darcy part and one for the saturation part) as well as set
1345 * to the right sizes the various objects required for the linear algebra in
1346 * this program. Its basic operations are similar to what @ref step_31 "step-31" did.
1347 *
1348
1349 *
1350 * The body of the function first enumerates all degrees of freedom for the
1351 * Darcy and saturation systems. For the Darcy part, degrees of freedom are
1352 * then sorted to ensure that velocities precede pressure DoFs so that we
1353 * can partition the Darcy matrix into a @f$2 \times 2@f$ matrix.
1354 *
1355
1356 *
1357 * Then, we need to incorporate hanging node constraints and Dirichlet
1358 * boundary value constraints into darcy_preconditioner_constraints. The
1359 * boundary condition constraints are only set on the pressure component
1360 * since the Schur complement preconditioner that corresponds to the porous
1361 * media flow operator in non-mixed form, @f$-\nabla \cdot [\mathbf K
1362 * \lambda_t(S)]\nabla@f$, acts only on the pressure variable. Therefore, we
1363 * use a component_mask that filters out the velocity component, so that the
1364 * condensation is performed on pressure degrees of freedom only.
1365 *
1366
1367 *
1368 * After having done so, we count the number of degrees of freedom in the
1369 * various blocks. This information is then used to create the sparsity
1370 * pattern for the Darcy and saturation system matrices as well as the
1371 * preconditioner matrix from which we build the Darcy preconditioner. As in
1372 * @ref step_31 "step-31", we choose to create the pattern using the blocked version of
1373 * DynamicSparsityPattern. So, for this, we follow the same way as @ref step_31 "step-31"
1374 * did and we don't have to repeat descriptions again for the rest of the
1375 * member function.
1376 *
1377 * @code
1378 *   template <int dim>
1379 *   void TwoPhaseFlowProblem<dim>::setup_dofs()
1380 *   {
1381 *   std::vector<unsigned int> darcy_block_component(dim + 1, 0);
1382 *   darcy_block_component[dim] = 1;
1383 *   {
1384 *   darcy_dof_handler.distribute_dofs(darcy_fe);
1385 *   DoFRenumbering::Cuthill_McKee(darcy_dof_handler);
1386 *   DoFRenumbering::component_wise(darcy_dof_handler, darcy_block_component);
1387 *  
1388 *   darcy_constraints.clear();
1389 *   DoFTools::make_hanging_node_constraints(darcy_dof_handler,
1390 *   darcy_constraints);
1391 *   darcy_constraints.close();
1392 *   }
1393 *   {
1394 *   saturation_dof_handler.distribute_dofs(saturation_fe);
1395 *  
1396 *   saturation_constraints.clear();
1397 *   DoFTools::make_hanging_node_constraints(saturation_dof_handler,
1398 *   saturation_constraints);
1399 *   saturation_constraints.close();
1400 *   }
1401 *   {
1402 *   darcy_preconditioner_constraints.clear();
1403 *  
1404 *   const FEValuesExtractors::Scalar pressure(dim);
1405 *  
1406 *   DoFTools::make_hanging_node_constraints(darcy_dof_handler,
1407 *   darcy_preconditioner_constraints);
1408 *   DoFTools::make_zero_boundary_constraints(darcy_dof_handler,
1409 *   darcy_preconditioner_constraints,
1410 *   darcy_fe.component_mask(
1411 *   pressure));
1412 *  
1413 *   darcy_preconditioner_constraints.close();
1414 *   }
1415 *  
1416 *  
1417 *   const std::vector<types::global_dof_index> darcy_dofs_per_block =
1418 *   DoFTools::count_dofs_per_fe_block(darcy_dof_handler,
1419 *   darcy_block_component);
1420 *   const types::global_dof_index n_u = darcy_dofs_per_block[0],
1421 *   n_p = darcy_dofs_per_block[1],
1422 *   n_s = saturation_dof_handler.n_dofs();
1423 *  
1424 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
1425 *   << " (on " << triangulation.n_levels() << " levels)" << std::endl
1426 *   << "Number of degrees of freedom: " << n_u + n_p + n_s << " ("
1427 *   << n_u << '+' << n_p << '+' << n_s << ')' << std::endl
1428 *   << std::endl;
1429 *  
1430 *   {
1431 *   darcy_matrix.clear();
1432 *  
1433 *   BlockDynamicSparsityPattern dsp(darcy_dofs_per_block,
1434 *   darcy_dofs_per_block);
1435 *  
1436 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
1437 *   for (unsigned int c = 0; c < dim + 1; ++c)
1438 *   for (unsigned int d = 0; d < dim + 1; ++d)
1439 *   if (!((c == dim) && (d == dim)))
1440 *   coupling[c][d] = DoFTools::always;
1441 *   else
1442 *   coupling[c][d] = DoFTools::none;
1443 *  
1444 *  
1445 *   DoFTools::make_sparsity_pattern(
1446 *   darcy_dof_handler, coupling, dsp, darcy_constraints, false);
1447 *  
1448 *   darcy_matrix.reinit(dsp);
1449 *   }
1450 *  
1451 *   {
1452 *   top_left_preconditioner.reset();
1453 *   bottom_right_preconditioner.reset();
1454 *   darcy_preconditioner_matrix.clear();
1455 *  
1456 *   BlockDynamicSparsityPattern dsp(darcy_dofs_per_block,
1457 *   darcy_dofs_per_block);
1458 *  
1459 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
1460 *   for (unsigned int c = 0; c < dim + 1; ++c)
1461 *   for (unsigned int d = 0; d < dim + 1; ++d)
1462 *   if (c == d)
1463 *   coupling[c][d] = DoFTools::always;
1464 *   else
1465 *   coupling[c][d] = DoFTools::none;
1466 *  
1467 *   DoFTools::make_sparsity_pattern(
1468 *   darcy_dof_handler, coupling, dsp, darcy_constraints, false);
1469 *  
1470 *   darcy_preconditioner_matrix.reinit(dsp);
1471 *   }
1472 *  
1473 *  
1474 *   {
1475 *   saturation_matrix.clear();
1476 *  
1477 *   DynamicSparsityPattern dsp(n_s, n_s);
1478 *  
1479 *   DoFTools::make_sparsity_pattern(saturation_dof_handler,
1480 *   dsp,
1481 *   saturation_constraints,
1482 *   false);
1483 *  
1484 *  
1485 *   saturation_matrix.reinit(dsp);
1486 *   }
1487 *  
1488 *   const std::vector<IndexSet> darcy_partitioning = {complete_index_set(n_u),
1489 *   complete_index_set(n_p)};
1490 *  
1491 *   darcy_solution.reinit(darcy_partitioning, MPI_COMM_WORLD);
1492 *  
1493 *   last_computed_darcy_solution.reinit(darcy_partitioning, MPI_COMM_WORLD);
1494 *  
1495 *   second_last_computed_darcy_solution.reinit(darcy_partitioning,
1496 *   MPI_COMM_WORLD);
1497 *  
1498 *   darcy_rhs.reinit(darcy_partitioning, MPI_COMM_WORLD);
1499 *  
1500 *   const IndexSet saturation_partitioning = complete_index_set(n_s);
1501 *   saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1502 *   old_saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1503 *   old_old_saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1504 *  
1505 *   saturation_matching_last_computed_darcy_solution.reinit(
1506 *   saturation_partitioning, MPI_COMM_WORLD);
1507 *  
1508 *   saturation_rhs.reinit(saturation_partitioning, MPI_COMM_WORLD);
1509 *   }
1510 *  
1511 *  
1512 * @endcode
1513 *
1514 *
1515 * <a name="step_43-Assemblingmatricesandpreconditioners"></a>
1516 * <h3>Assembling matrices and preconditioners</h3>
1517 *
1518
1519 *
1520 * The next few functions are devoted to setting up the various system and
1521 * preconditioner matrices and right hand sides that we have to deal with in
1522 * this program.
1523 *
1524
1525 *
1526 *
1527 * <a name="step_43-TwoPhaseFlowProblemdimassemble_darcy_preconditioner"></a>
1528 * <h4>TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner</h4>
1529 *
1530
1531 *
1532 * This function assembles the matrix we use for preconditioning the Darcy
1533 * system. What we need are a vector @ref GlossMassMatrix "mass matrix" weighted by
1534 * @f$\left(\mathbf{K} \lambda_t\right)^{-1}@f$ on the velocity components and a
1535 * mass matrix weighted by @f$\left(\mathbf{K} \lambda_t\right)@f$ on the
1536 * pressure component. We start by generating a quadrature object of
1537 * appropriate order, the FEValues object that can give values and gradients
1538 * at the quadrature points (together with quadrature weights). Next we
1539 * create data structures for the cell matrix and the relation between local
1540 * and global DoFs. The vectors phi_u and grad_phi_p are going to hold the
1541 * values of the basis functions in order to faster build up the local
1542 * matrices, as was already done in @ref step_22 "step-22". Before we start the loop over
1543 * all active cells, we have to specify which components are pressure and
1544 * which are velocity.
1545 *
1546
1547 *
1548 * The creation of the local matrix is rather simple. There are only a term
1549 * weighted by @f$\left(\mathbf{K} \lambda_t\right)^{-1}@f$ (on the velocity)
1550 * and a Laplace matrix weighted by @f$\left(\mathbf{K} \lambda_t\right)@f$ to
1551 * be generated, so the creation of the local matrix is done in essentially
1552 * two lines. Since the material model functions at the top of this file
1553 * only provide the inverses of the permeability and mobility, we have to
1554 * compute @f$\mathbf K@f$ and @f$\lambda_t@f$ by hand from the given values, once
1555 * per quadrature point.
1556 *
1557
1558 *
1559 * Once the local matrix is ready (loop over rows and columns in the local
1560 * matrix on each quadrature point), we get the local DoF indices and write
1561 * the local information into the global matrix. We do this by directly
1562 * applying the constraints (i.e. darcy_preconditioner_constraints) that
1563 * takes care of hanging node and zero Dirichlet boundary condition
1564 * constraints. By doing so, we don't have to do that afterwards, and we
1565 * later don't have to use AffineConstraints::condense and
1566 * MatrixTools::apply_boundary_values, both functions that would need to
1567 * modify matrix and vector entries and so are difficult to write for the
1568 * Trilinos classes where we don't immediately have access to individual
1569 * memory locations.
1570 *
1571 * @code
1572 *   template <int dim>
1573 *   void TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner()
1574 *   {
1575 *   std::cout << " Rebuilding darcy preconditioner..." << std::endl;
1576 *  
1577 *   darcy_preconditioner_matrix = 0;
1578 *  
1579 *   const QGauss<dim> quadrature_formula(darcy_degree + 2);
1580 *   FEValues<dim> darcy_fe_values(darcy_fe,
1581 *   quadrature_formula,
1582 *   update_JxW_values | update_values |
1583 *   update_gradients |
1584 *   update_quadrature_points);
1585 *   FEValues<dim> saturation_fe_values(saturation_fe,
1586 *   quadrature_formula,
1587 *   update_values);
1588 *  
1589 *   const unsigned int dofs_per_cell = darcy_fe.n_dofs_per_cell();
1590 *   const unsigned int n_q_points = quadrature_formula.size();
1591 *  
1592 *   std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1593 *  
1594 *   std::vector<double> old_saturation_values(n_q_points);
1595 *  
1596 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1597 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1598 *  
1599 *   std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
1600 *   std::vector<Tensor<1, dim>> grad_phi_p(dofs_per_cell);
1601 *  
1602 *   const FEValuesExtractors::Vector velocities(0);
1603 *   const FEValuesExtractors::Scalar pressure(dim);
1604 *  
1605 *   auto cell = darcy_dof_handler.begin_active();
1606 *   const auto endc = darcy_dof_handler.end();
1607 *   auto saturation_cell = saturation_dof_handler.begin_active();
1608 *  
1609 *   for (; cell != endc; ++cell, ++saturation_cell)
1610 *   {
1611 *   darcy_fe_values.reinit(cell);
1612 *   saturation_fe_values.reinit(saturation_cell);
1613 *  
1614 *   local_matrix = 0;
1615 *  
1616 *   saturation_fe_values.get_function_values(old_saturation_solution,
1617 *   old_saturation_values);
1618 *  
1619 *   k_inverse.value_list(darcy_fe_values.get_quadrature_points(),
1620 *   k_inverse_values);
1621 *  
1622 *   for (unsigned int q = 0; q < n_q_points; ++q)
1623 *   {
1624 *   const double old_s = old_saturation_values[q];
1625 *  
1626 *   const double inverse_mobility = mobility_inverse(old_s, viscosity);
1627 *   const double mobility = 1.0 / inverse_mobility;
1628 *   const Tensor<2, dim> permeability = invert(k_inverse_values[q]);
1629 *  
1630 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1631 *   {
1632 *   phi_u[k] = darcy_fe_values[velocities].value(k, q);
1633 *   grad_phi_p[k] = darcy_fe_values[pressure].gradient(k, q);
1634 *   }
1635 *  
1636 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1637 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1638 *   {
1639 *   local_matrix(i, j) +=
1640 *   (k_inverse_values[q] * inverse_mobility * phi_u[i] *
1641 *   phi_u[j] +
1642 *   permeability * mobility * grad_phi_p[i] * grad_phi_p[j]) *
1643 *   darcy_fe_values.JxW(q);
1644 *   }
1645 *   }
1646 *  
1647 *   cell->get_dof_indices(local_dof_indices);
1648 *   darcy_preconditioner_constraints.distribute_local_to_global(
1649 *   local_matrix, local_dof_indices, darcy_preconditioner_matrix);
1650 *   }
1651 *  
1652 *   darcy_preconditioner_matrix.compress(VectorOperation::add);
1653 *   }
1654 *  
1655 *  
1656 * @endcode
1657 *
1658 *
1659 * <a name="step_43-TwoPhaseFlowProblemdimbuild_darcy_preconditioner"></a>
1660 * <h4>TwoPhaseFlowProblem<dim>::build_darcy_preconditioner</h4>
1661 *
1662
1663 *
1664 * After calling the above functions to assemble the preconditioner matrix,
1665 * this function generates the inner preconditioners that are going to be
1666 * used for the Schur complement block preconditioner. The preconditioners
1667 * need to be regenerated at every saturation time step since they depend on
1668 * the saturation @f$S@f$ that varies with time.
1669 *
1670
1671 *
1672 * In here, we set up the preconditioner for the velocity-velocity matrix
1673 * @f$\mathbf{M}^{\mathbf{u}}@f$ and the Schur complement @f$\mathbf{S}@f$. As
1674 * explained in the introduction, we are going to use an IC preconditioner
1675 * based on the vector matrix @f$\mathbf{M}^{\mathbf{u}}@f$ and another based on
1676 * the scalar Laplace matrix @f$\tilde{\mathbf{S}}^p@f$ (which is spectrally
1677 * close to the Schur complement of the Darcy matrix). Usually, the
1678 * TrilinosWrappers::PreconditionIC class can be seen as a good black-box
1679 * preconditioner which does not need any special knowledge of the matrix
1680 * structure and/or the operator that's behind it.
1681 *
1682 * @code
1683 *   template <int dim>
1684 *   void TwoPhaseFlowProblem<dim>::build_darcy_preconditioner()
1685 *   {
1686 *   assemble_darcy_preconditioner();
1687 *  
1688 *   #ifdef DEAL_II_TRILINOS_WITH_EPETRA
1689 *   top_left_preconditioner = std::make_shared<PreconditionType>();
1690 *   #else
1691 *   top_left_preconditioner = std::make_shared<PreconditionType>("FAST_IC");
1692 *   #endif
1693 *   top_left_preconditioner->initialize(
1694 *   darcy_preconditioner_matrix.block(0, 0));
1695 *  
1696 *   #ifdef DEAL_II_TRILINOS_WITH_EPETRA
1697 *   bottom_right_preconditioner = std::make_shared<PreconditionType>();
1698 *   #else
1699 *   bottom_right_preconditioner = std::make_shared<PreconditionType>("FAST_IC");
1700 *   #endif
1701 *   bottom_right_preconditioner->initialize(
1702 *   darcy_preconditioner_matrix.block(1, 1));
1703 *   }
1704 *  
1705 *  
1706 * @endcode
1707 *
1708 *
1709 * <a name="step_43-TwoPhaseFlowProblemdimassemble_darcy_system"></a>
1710 * <h4>TwoPhaseFlowProblem<dim>::assemble_darcy_system</h4>
1711 *
1712
1713 *
1714 * This is the function that assembles the linear system for the Darcy
1715 * system.
1716 *
1717
1718 *
1719 * Regarding the technical details of implementation, the procedures are
1720 * similar to those in @ref step_22 "step-22" and @ref step_31 "step-31". We reset matrix and vector,
1721 * create a quadrature formula on the cells, and then create the respective
1722 * FEValues object.
1723 *
1724
1725 *
1726 * There is one thing that needs to be commented: since we have a separate
1727 * finite element and DoFHandler for the saturation, we need to generate a
1728 * second FEValues object for the proper evaluation of the saturation
1729 * solution. This isn't too complicated to realize here: just use the
1730 * saturation structures and set an update flag for the basis function
1731 * values which we need for evaluation of the saturation solution. The only
1732 * important part to remember here is that the same quadrature formula is
1733 * used for both FEValues objects to ensure that we get matching information
1734 * when we loop over the quadrature points of the two objects.
1735 *
1736
1737 *
1738 * The declarations proceed with some shortcuts for array sizes, the
1739 * creation of the local matrix, right hand side as well as the vector for
1740 * the indices of the local dofs compared to the global system.
1741 *
1742 * @code
1743 *   template <int dim>
1744 *   void TwoPhaseFlowProblem<dim>::assemble_darcy_system()
1745 *   {
1746 *   darcy_matrix = 0;
1747 *   darcy_rhs = 0;
1748 *  
1749 *   const QGauss<dim> quadrature_formula(darcy_degree + 2);
1750 *   const QGauss<dim - 1> face_quadrature_formula(darcy_degree + 2);
1751 *  
1752 *   FEValues<dim> darcy_fe_values(darcy_fe,
1753 *   quadrature_formula,
1754 *   update_values | update_gradients |
1755 *   update_quadrature_points |
1756 *   update_JxW_values);
1757 *  
1758 *   FEValues<dim> saturation_fe_values(saturation_fe,
1759 *   quadrature_formula,
1760 *   update_values);
1761 *  
1762 *   FEFaceValues<dim> darcy_fe_face_values(darcy_fe,
1763 *   face_quadrature_formula,
1764 *   update_values |
1765 *   update_normal_vectors |
1766 *   update_quadrature_points |
1767 *   update_JxW_values);
1768 *  
1769 *   const unsigned int dofs_per_cell = darcy_fe.n_dofs_per_cell();
1770 *  
1771 *   const unsigned int n_q_points = quadrature_formula.size();
1772 *   const unsigned int n_face_q_points = face_quadrature_formula.size();
1773 *  
1774 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1775 *   Vector<double> local_rhs(dofs_per_cell);
1776 *  
1777 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1778 *  
1779 *   const Functions::ZeroFunction<dim> pressure_right_hand_side;
1780 *   const PressureBoundaryValues<dim> pressure_boundary_values;
1781 *  
1782 *   std::vector<double> pressure_rhs_values(n_q_points);
1783 *   std::vector<double> boundary_values(n_face_q_points);
1784 *   std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1785 *  
1786 * @endcode
1787 *
1788 * Next we need a vector that will contain the values of the saturation
1789 * solution at the previous time level at the quadrature points to
1790 * assemble the saturation dependent coefficients in the Darcy equations.
1791 *
1792
1793 *
1794 * The set of vectors we create next hold the evaluations of the basis
1795 * functions as well as their gradients that will be used for creating the
1796 * matrices. Putting these into their own arrays rather than asking the
1797 * FEValues object for this information each time it is needed is an
1798 * optimization to accelerate the assembly process, see @ref step_22 "step-22" for
1799 * details.
1800 *
1801
1802 *
1803 * The last two declarations are used to extract the individual blocks
1804 * (velocity, pressure, saturation) from the total FE system.
1805 *
1806 * @code
1807 *   std::vector<double> old_saturation_values(n_q_points);
1808 *  
1809 *   std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
1810 *   std::vector<double> div_phi_u(dofs_per_cell);
1811 *   std::vector<double> phi_p(dofs_per_cell);
1812 *  
1813 *   const FEValuesExtractors::Vector velocities(0);
1814 *   const FEValuesExtractors::Scalar pressure(dim);
1815 *  
1816 * @endcode
1817 *
1818 * Now start the loop over all cells in the problem. We are working on two
1819 * different DoFHandlers for this assembly routine, so we must have two
1820 * different cell iterators for the two objects in use. This might seem a
1821 * bit peculiar, but since both the Darcy system and the saturation system
1822 * use the same grid we can assume that the two iterators run in sync over
1823 * the cells of the two DoFHandler objects.
1824 *
1825
1826 *
1827 * The first statements within the loop are again all very familiar, doing
1828 * the update of the finite element data as specified by the update flags,
1829 * zeroing out the local arrays and getting the values of the old solution
1830 * at the quadrature points. At this point we also have to get the values
1831 * of the saturation function of the previous time step at the quadrature
1832 * points. To this end, we can use the FEValues::get_function_values
1833 * (previously already used in @ref step_9 "step-9", @ref step_14 "step-14" and @ref step_15 "step-15"), a function
1834 * that takes a solution vector and returns a list of function values at
1835 * the quadrature points of the present cell. In fact, it returns the
1836 * complete vector-valued solution at each quadrature point, i.e. not only
1837 * the saturation but also the velocities and pressure.
1838 *
1839
1840 *
1841 * Then we are ready to loop over the quadrature points on the cell to do
1842 * the integration. The formula for this follows in a straightforward way
1843 * from what has been discussed in the introduction.
1844 *
1845
1846 *
1847 * Once this is done, we start the loop over the rows and columns of the
1848 * local matrix and feed the matrix with the relevant products.
1849 *
1850
1851 *
1852 * The last step in the loop over all cells is to enter the local
1853 * contributions into the global matrix and vector structures to the
1854 * positions specified in local_dof_indices. Again, we let the
1855 * AffineConstraints class do the insertion of the cell matrix
1856 * elements to the global matrix, which already condenses the hanging node
1857 * constraints.
1858 *
1859 * @code
1860 *   auto cell = darcy_dof_handler.begin_active();
1861 *   const auto endc = darcy_dof_handler.end();
1862 *   auto saturation_cell = saturation_dof_handler.begin_active();
1863 *  
1864 *   for (; cell != endc; ++cell, ++saturation_cell)
1865 *   {
1866 *   darcy_fe_values.reinit(cell);
1867 *   saturation_fe_values.reinit(saturation_cell);
1868 *  
1869 *   local_matrix = 0;
1870 *   local_rhs = 0;
1871 *  
1872 *   saturation_fe_values.get_function_values(old_saturation_solution,
1873 *   old_saturation_values);
1874 *  
1875 *   pressure_right_hand_side.value_list(
1876 *   darcy_fe_values.get_quadrature_points(), pressure_rhs_values);
1877 *   k_inverse.value_list(darcy_fe_values.get_quadrature_points(),
1878 *   k_inverse_values);
1879 *  
1880 *   for (unsigned int q = 0; q < n_q_points; ++q)
1881 *   {
1882 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1883 *   {
1884 *   phi_u[k] = darcy_fe_values[velocities].value(k, q);
1885 *   div_phi_u[k] = darcy_fe_values[velocities].divergence(k, q);
1886 *   phi_p[k] = darcy_fe_values[pressure].value(k, q);
1887 *   }
1888 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1889 *   {
1890 *   const double old_s = old_saturation_values[q];
1891 *   for (unsigned int j = 0; j <= i; ++j)
1892 *   {
1893 *   local_matrix(i, j) +=
1894 *   (phi_u[i] * k_inverse_values[q] *
1895 *   mobility_inverse(old_s, viscosity) * phi_u[j] -
1896 *   div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
1897 *   darcy_fe_values.JxW(q);
1898 *   }
1899 *  
1900 *   local_rhs(i) +=
1901 *   (-phi_p[i] * pressure_rhs_values[q]) * darcy_fe_values.JxW(q);
1902 *   }
1903 *   }
1904 *  
1905 *   for (const auto &face : cell->face_iterators())
1906 *   if (face->at_boundary())
1907 *   {
1908 *   darcy_fe_face_values.reinit(cell, face);
1909 *  
1910 *   pressure_boundary_values.value_list(
1911 *   darcy_fe_face_values.get_quadrature_points(), boundary_values);
1912 *  
1913 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1914 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1915 *   {
1916 *   const Tensor<1, dim> phi_i_u =
1917 *   darcy_fe_face_values[velocities].value(i, q);
1918 *  
1919 *   local_rhs(i) +=
1920 *   -(phi_i_u * darcy_fe_face_values.normal_vector(q) *
1921 *   boundary_values[q] * darcy_fe_face_values.JxW(q));
1922 *   }
1923 *   }
1924 *  
1925 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1926 *   for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
1927 *   local_matrix(i, j) = local_matrix(j, i);
1928 *  
1929 *   cell->get_dof_indices(local_dof_indices);
1930 *  
1931 *   darcy_constraints.distribute_local_to_global(
1932 *   local_matrix, local_rhs, local_dof_indices, darcy_matrix, darcy_rhs);
1933 *   }
1934 *  
1935 *   darcy_matrix.compress(VectorOperation::add);
1936 *   darcy_rhs.compress(VectorOperation::add);
1937 *   }
1938 *  
1939 *  
1940 * @endcode
1941 *
1942 *
1943 * <a name="step_43-TwoPhaseFlowProblemdimassemble_saturation_system"></a>
1944 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_system</h4>
1945 *
1946
1947 *
1948 * This function is to assemble the linear system for the saturation
1949 * transport equation. It calls, if necessary, two other member functions:
1950 * assemble_saturation_matrix() and assemble_saturation_rhs(). The former
1951 * function then assembles the saturation matrix that only needs to be
1952 * changed occasionally. On the other hand, the latter function that
1953 * assembles the right hand side must be called at every saturation time
1954 * step.
1955 *
1956 * @code
1957 *   template <int dim>
1958 *   void TwoPhaseFlowProblem<dim>::assemble_saturation_system()
1959 *   {
1960 *   if (rebuild_saturation_matrix == true)
1961 *   {
1962 *   saturation_matrix = 0;
1963 *   assemble_saturation_matrix();
1964 *   }
1965 *  
1966 *   saturation_rhs = 0;
1967 *   assemble_saturation_rhs();
1968 *   }
1969 *  
1970 *  
1971 *  
1972 * @endcode
1973 *
1974 *
1975 * <a name="step_43-TwoPhaseFlowProblemdimassemble_saturation_matrix"></a>
1976 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_matrix</h4>
1977 *
1978
1979 *
1980 * This function is easily understood since it only forms a simple mass
1981 * matrix for the left hand side of the saturation linear system by basis
1982 * functions phi_i_s and phi_j_s only. Finally, as usual, we enter the local
1983 * contribution into the global matrix by specifying the position in
1984 * local_dof_indices. This is done by letting the AffineConstraints class do
1985 * the insertion of the cell matrix elements to the global matrix, which
1986 * already condenses the hanging node constraints.
1987 *
1988 * @code
1989 *   template <int dim>
1990 *   void TwoPhaseFlowProblem<dim>::assemble_saturation_matrix()
1991 *   {
1992 *   const QGauss<dim> quadrature_formula(saturation_degree + 2);
1993 *  
1994 *   FEValues<dim> saturation_fe_values(saturation_fe,
1995 *   quadrature_formula,
1996 *   update_values | update_JxW_values);
1997 *  
1998 *   const unsigned int dofs_per_cell = saturation_fe.n_dofs_per_cell();
1999 *  
2000 *   const unsigned int n_q_points = quadrature_formula.size();
2001 *  
2002 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2003 *   Vector<double> local_rhs(dofs_per_cell);
2004 *  
2005 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2006 *  
2007 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2008 *   {
2009 *   saturation_fe_values.reinit(cell);
2010 *   local_matrix = 0;
2011 *   local_rhs = 0;
2012 *  
2013 *   for (unsigned int q = 0; q < n_q_points; ++q)
2014 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2015 *   {
2016 *   const double phi_i_s = saturation_fe_values.shape_value(i, q);
2017 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
2018 *   {
2019 *   const double phi_j_s = saturation_fe_values.shape_value(j, q);
2020 *   local_matrix(i, j) +=
2021 *   porosity * phi_i_s * phi_j_s * saturation_fe_values.JxW(q);
2022 *   }
2023 *   }
2024 *   cell->get_dof_indices(local_dof_indices);
2025 *  
2026 *   saturation_constraints.distribute_local_to_global(local_matrix,
2027 *   local_dof_indices,
2028 *   saturation_matrix);
2029 *   }
2030 *  
2031 *   saturation_matrix.compress(VectorOperation::add);
2032 *   }
2033 *  
2034 *  
2035 *  
2036 * @endcode
2037 *
2038 *
2039 * <a name="step_43-TwoPhaseFlowProblemdimassemble_saturation_rhs"></a>
2040 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs</h4>
2041 *
2042
2043 *
2044 * This function is to assemble the right hand side of the saturation
2045 * transport equation. Before going about it, we have to create two FEValues
2046 * objects for the Darcy and saturation systems respectively and, in
2047 * addition, two FEFaceValues objects for the two systems because we have a
2048 * boundary integral term in the weak form of saturation equation. For the
2049 * FEFaceValues object of the saturation system, we also require normal
2050 * vectors, which we request using the update_normal_vectors flag.
2051 *
2052
2053 *
2054 * Next, before looping over all the cells, we have to compute some
2055 * parameters (e.g. global_u_infty, global_S_variation, and
2056 * global_Omega_diameter) that the artificial viscosity @f$\nu@f$ needs. This is
2057 * largely the same as was done in @ref step_31 "step-31", so you may see there for more
2058 * information.
2059 *
2060
2061 *
2062 * The real works starts with the loop over all the saturation and Darcy
2063 * cells to put the local contributions into the global vector. In this
2064 * loop, in order to simplify the implementation, we split some of the work
2065 * into two helper functions: assemble_saturation_rhs_cell_term and
2066 * assemble_saturation_rhs_boundary_term. We note that we insert cell or
2067 * boundary contributions into the global vector in the two functions rather
2068 * than in this present function.
2069 *
2070 * @code
2071 *   template <int dim>
2072 *   void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs()
2073 *   {
2074 *   const QGauss<dim> quadrature_formula(saturation_degree + 2);
2075 *   const QGauss<dim - 1> face_quadrature_formula(saturation_degree + 2);
2076 *  
2077 *   FEValues<dim> saturation_fe_values(saturation_fe,
2078 *   quadrature_formula,
2079 *   update_values | update_gradients |
2080 *   update_quadrature_points |
2081 *   update_JxW_values);
2082 *   FEValues<dim> darcy_fe_values(darcy_fe, quadrature_formula, update_values);
2083 *   FEFaceValues<dim> saturation_fe_face_values(saturation_fe,
2084 *   face_quadrature_formula,
2085 *   update_values |
2086 *   update_normal_vectors |
2087 *   update_quadrature_points |
2088 *   update_JxW_values);
2089 *   FEFaceValues<dim> darcy_fe_face_values(darcy_fe,
2090 *   face_quadrature_formula,
2091 *   update_values);
2092 *   FEFaceValues<dim> saturation_fe_face_values_neighbor(
2093 *   saturation_fe, face_quadrature_formula, update_values);
2094 *  
2095 *   const unsigned int dofs_per_cell =
2096 *   saturation_dof_handler.get_fe().n_dofs_per_cell();
2097 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2098 *  
2099 *   const double global_max_u_F_prime = get_max_u_F_prime();
2100 *   const std::pair<double, double> global_S_range =
2101 *   get_extrapolated_saturation_range();
2102 *   const double global_S_variation =
2103 *   global_S_range.second - global_S_range.first;
2104 *  
2105 *   auto cell = saturation_dof_handler.begin_active();
2106 *   const auto endc = saturation_dof_handler.end();
2107 *   auto darcy_cell = darcy_dof_handler.begin_active();
2108 *   for (; cell != endc; ++cell, ++darcy_cell)
2109 *   {
2110 *   saturation_fe_values.reinit(cell);
2111 *   darcy_fe_values.reinit(darcy_cell);
2112 *  
2113 *   cell->get_dof_indices(local_dof_indices);
2114 *  
2115 *   assemble_saturation_rhs_cell_term(saturation_fe_values,
2116 *   darcy_fe_values,
2117 *   global_max_u_F_prime,
2118 *   global_S_variation,
2119 *   local_dof_indices);
2120 *  
2121 *   for (const auto &face : cell->face_iterators())
2122 *   if (face->at_boundary())
2123 *   {
2124 *   darcy_fe_face_values.reinit(darcy_cell, face);
2125 *   saturation_fe_face_values.reinit(cell, face);
2126 *   assemble_saturation_rhs_boundary_term(saturation_fe_face_values,
2127 *   darcy_fe_face_values,
2128 *   local_dof_indices);
2129 *   }
2130 *   }
2131 *  
2132 *   saturation_rhs.compress(VectorOperation::add);
2133 *   }
2134 *  
2135 *  
2136 *  
2137 * @endcode
2138 *
2139 *
2140 * <a name="step_43-TwoPhaseFlowProblemdimassemble_saturation_rhs_cell_term"></a>
2141 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_cell_term</h4>
2142 *
2143
2144 *
2145 * This function takes care of integrating the cell terms of the right hand
2146 * side of the saturation equation, and then assembling it into the global
2147 * right hand side vector. Given the discussion in the introduction, the
2148 * form of these contributions is clear. The only tricky part is getting the
2149 * artificial viscosity and all that is necessary to compute it. The first
2150 * half of the function is devoted to this task.
2151 *
2152
2153 *
2154 * The last part of the function is copying the local contributions into the
2155 * global vector with position specified in local_dof_indices.
2156 *
2157 * @code
2158 *   template <int dim>
2159 *   void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_cell_term(
2160 *   const FEValues<dim> &saturation_fe_values,
2161 *   const FEValues<dim> &darcy_fe_values,
2162 *   const double global_max_u_F_prime,
2163 *   const double global_S_variation,
2164 *   const std::vector<types::global_dof_index> &local_dof_indices)
2165 *   {
2166 *   const unsigned int dofs_per_cell = saturation_fe_values.dofs_per_cell;
2167 *   const unsigned int n_q_points = saturation_fe_values.n_quadrature_points;
2168 *  
2169 *   std::vector<double> old_saturation_solution_values(n_q_points);
2170 *   std::vector<double> old_old_saturation_solution_values(n_q_points);
2171 *   std::vector<Tensor<1, dim>> old_grad_saturation_solution_values(n_q_points);
2172 *   std::vector<Tensor<1, dim>> old_old_grad_saturation_solution_values(
2173 *   n_q_points);
2174 *   std::vector<Vector<double>> present_darcy_solution_values(
2175 *   n_q_points, Vector<double>(dim + 1));
2176 *  
2177 *   saturation_fe_values.get_function_values(old_saturation_solution,
2178 *   old_saturation_solution_values);
2179 *   saturation_fe_values.get_function_values(
2180 *   old_old_saturation_solution, old_old_saturation_solution_values);
2181 *   saturation_fe_values.get_function_gradients(
2182 *   old_saturation_solution, old_grad_saturation_solution_values);
2183 *   saturation_fe_values.get_function_gradients(
2184 *   old_old_saturation_solution, old_old_grad_saturation_solution_values);
2185 *   darcy_fe_values.get_function_values(darcy_solution,
2186 *   present_darcy_solution_values);
2187 *  
2188 *   const double nu =
2189 *   compute_viscosity(old_saturation_solution_values,
2190 *   old_old_saturation_solution_values,
2191 *   old_grad_saturation_solution_values,
2192 *   old_old_grad_saturation_solution_values,
2193 *   present_darcy_solution_values,
2194 *   global_max_u_F_prime,
2195 *   global_S_variation,
2196 *   saturation_fe_values.get_cell()->diameter());
2197 *  
2198 *   Vector<double> local_rhs(dofs_per_cell);
2199 *  
2200 *   for (unsigned int q = 0; q < n_q_points; ++q)
2201 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2202 *   {
2203 *   const double old_s = old_saturation_solution_values[q];
2204 *   Tensor<1, dim> present_u;
2205 *   for (unsigned int d = 0; d < dim; ++d)
2206 *   present_u[d] = present_darcy_solution_values[q](d);
2207 *  
2208 *   const double phi_i_s = saturation_fe_values.shape_value(i, q);
2209 *   const Tensor<1, dim> grad_phi_i_s =
2210 *   saturation_fe_values.shape_grad(i, q);
2211 *  
2212 *   local_rhs(i) +=
2213 *   (time_step * fractional_flow(old_s, viscosity) * present_u *
2214 *   grad_phi_i_s -
2215 *   time_step * nu * old_grad_saturation_solution_values[q] *
2216 *   grad_phi_i_s +
2217 *   porosity * old_s * phi_i_s) *
2218 *   saturation_fe_values.JxW(q);
2219 *   }
2220 *  
2221 *   saturation_constraints.distribute_local_to_global(local_rhs,
2222 *   local_dof_indices,
2223 *   saturation_rhs);
2224 *   }
2225 *  
2226 *  
2227 * @endcode
2228 *
2229 *
2230 * <a name="step_43-TwoPhaseFlowProblemdimassemble_saturation_rhs_boundary_term"></a>
2231 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_boundary_term</h4>
2232 *
2233
2234 *
2235 * The next function is responsible for the boundary integral terms in the
2236 * right hand side form of the saturation equation. For these, we have to
2237 * compute the upwinding flux on the global boundary faces, i.e. we impose
2238 * Dirichlet boundary conditions weakly only on inflow parts of the global
2239 * boundary. As before, this has been described in @ref step_21 "step-21" so we refrain
2240 * from giving more descriptions about that.
2241 *
2242 * @code
2243 *   template <int dim>
2244 *   void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_boundary_term(
2245 *   const FEFaceValues<dim> &saturation_fe_face_values,
2246 *   const FEFaceValues<dim> &darcy_fe_face_values,
2247 *   const std::vector<types::global_dof_index> &local_dof_indices)
2248 *   {
2249 *   const unsigned int dofs_per_cell = saturation_fe_face_values.dofs_per_cell;
2250 *   const unsigned int n_face_q_points =
2251 *   saturation_fe_face_values.n_quadrature_points;
2252 *  
2253 *   Vector<double> local_rhs(dofs_per_cell);
2254 *  
2255 *   std::vector<double> old_saturation_solution_values_face(n_face_q_points);
2256 *   std::vector<Vector<double>> present_darcy_solution_values_face(
2257 *   n_face_q_points, Vector<double>(dim + 1));
2258 *   std::vector<double> neighbor_saturation(n_face_q_points);
2259 *  
2260 *   saturation_fe_face_values.get_function_values(
2261 *   old_saturation_solution, old_saturation_solution_values_face);
2262 *   darcy_fe_face_values.get_function_values(
2263 *   darcy_solution, present_darcy_solution_values_face);
2264 *  
2265 *   SaturationBoundaryValues<dim> saturation_boundary_values;
2266 *   saturation_boundary_values.value_list(
2267 *   saturation_fe_face_values.get_quadrature_points(), neighbor_saturation);
2268 *  
2269 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
2270 *   {
2271 *   Tensor<1, dim> present_u_face;
2272 *   for (unsigned int d = 0; d < dim; ++d)
2273 *   present_u_face[d] = present_darcy_solution_values_face[q](d);
2274 *  
2275 *   const double normal_flux =
2276 *   present_u_face * saturation_fe_face_values.normal_vector(q);
2277 *  
2278 *   const bool is_outflow_q_point = (normal_flux >= 0);
2279 *  
2280 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2281 *   local_rhs(i) -=
2282 *   time_step * normal_flux *
2283 *   fractional_flow((is_outflow_q_point == true ?
2284 *   old_saturation_solution_values_face[q] :
2285 *   neighbor_saturation[q]),
2286 *   viscosity) *
2287 *   saturation_fe_face_values.shape_value(i, q) *
2288 *   saturation_fe_face_values.JxW(q);
2289 *   }
2290 *   saturation_constraints.distribute_local_to_global(local_rhs,
2291 *   local_dof_indices,
2292 *   saturation_rhs);
2293 *   }
2294 *  
2295 *  
2296 * @endcode
2297 *
2298 *
2299 * <a name="step_43-TwoPhaseFlowProblemdimsolve"></a>
2300 * <h3>TwoPhaseFlowProblem<dim>::solve</h3>
2301 *
2302
2303 *
2304 * This function implements the operator splitting algorithm, i.e. in each
2305 * time step it either re-computes the solution of the Darcy system or
2306 * extrapolates velocity/pressure from previous time steps, then determines
2307 * the size of the time step, and then updates the saturation variable. The
2308 * implementation largely follows similar code in @ref step_31 "step-31". It is, next to
2309 * the run() function, the central one in this program.
2310 *
2311
2312 *
2313 * At the beginning of the function, we ask whether to solve the
2314 * pressure-velocity part by evaluating the a posteriori criterion (see the
2315 * following function). If necessary, we will solve the pressure-velocity
2316 * part using the GMRES solver with the Schur complement block
2317 * preconditioner as is described in the introduction.
2318 *
2319 * @code
2320 *   template <int dim>
2321 *   void TwoPhaseFlowProblem<dim>::solve()
2322 *   {
2323 *   const bool solve_for_pressure_and_velocity =
2324 *   determine_whether_to_solve_for_pressure_and_velocity();
2325 *  
2326 *   if (solve_for_pressure_and_velocity == true)
2327 *   {
2328 *   std::cout << " Solving Darcy (pressure-velocity) system..."
2329 *   << std::endl;
2330 *  
2331 *   assemble_darcy_system();
2332 *   build_darcy_preconditioner();
2333 *  
2334 *   {
2335 *   const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
2336 *   PreconditionType>
2337 *   mp_inverse(darcy_preconditioner_matrix.block(1, 1),
2338 *   *bottom_right_preconditioner);
2339 *  
2340 *   const LinearSolvers::BlockSchurPreconditioner<PreconditionType,
2341 *   PreconditionType>
2342 *   preconditioner(darcy_matrix, mp_inverse, *top_left_preconditioner);
2343 *  
2344 *   SolverControl solver_control(darcy_matrix.m(),
2345 *   1e-16 * darcy_rhs.l2_norm());
2346 *  
2347 *   SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
2348 *   solver_control,
2349 *   SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(
2350 *   100));
2351 *  
2352 *   for (unsigned int i = 0; i < darcy_solution.size(); ++i)
2353 *   if (darcy_constraints.is_constrained(i))
2354 *   darcy_solution(i) = 0;
2355 *  
2356 *   gmres.solve(darcy_matrix, darcy_solution, darcy_rhs, preconditioner);
2357 *  
2358 *   darcy_constraints.distribute(darcy_solution);
2359 *  
2360 *   std::cout << " ..." << solver_control.last_step()
2361 *   << " GMRES iterations." << std::endl;
2362 *   }
2363 *  
2364 *   {
2365 *   second_last_computed_darcy_solution = last_computed_darcy_solution;
2366 *   last_computed_darcy_solution = darcy_solution;
2367 *  
2368 *   saturation_matching_last_computed_darcy_solution =
2369 *   saturation_solution;
2370 *   }
2371 *   }
2372 * @endcode
2373 *
2374 * On the other hand, if we have decided that we don't want to compute the
2375 * solution of the Darcy system for the current time step, then we need to
2376 * simply extrapolate the previous two Darcy solutions to the same time as
2377 * we would have computed the velocity/pressure at. We do a simple linear
2378 * extrapolation, i.e. given the current length @f$dt@f$ of the macro time
2379 * step from the time when we last computed the Darcy solution to now
2380 * (given by <code>current_macro_time_step</code>), and @f$DT@f$ the length of
2381 * the last macro time step (given by <code>old_macro_time_step</code>),
2382 * then we get @f$u^\ast = u_p + dt \frac{u_p-u_{pp}}{DT} = (1+dt/DT)u_p -
2383 * dt/DT u_{pp}@f$, where @f$u_p@f$ and @f$u_{pp}@f$ are the last two computed Darcy
2384 * solutions. We can implement this formula using just two lines of code.
2385 *
2386
2387 *
2388 * Note that the algorithm here only works if we have at least two
2389 * previously computed Darcy solutions from which we can extrapolate to
2390 * the current time, and this is ensured by requiring re-computation of
2391 * the Darcy solution for the first 2 time steps.
2392 *
2393 * @code
2394 *   else
2395 *   {
2396 *   darcy_solution = last_computed_darcy_solution;
2397 *   darcy_solution.sadd(1 + current_macro_time_step / old_macro_time_step,
2398 *   -current_macro_time_step / old_macro_time_step,
2399 *   second_last_computed_darcy_solution);
2400 *   }
2401 *  
2402 *  
2403 * @endcode
2404 *
2405 * With the so computed velocity vector, compute the optimal time step
2406 * based on the CFL criterion discussed in the introduction...
2407 *
2408 * @code
2409 *   {
2410 *   old_time_step = time_step;
2411 *  
2412 *   const double max_u_F_prime = get_max_u_F_prime();
2413 *   if (max_u_F_prime > 0)
2414 *   time_step = porosity * GridTools::minimal_cell_diameter(triangulation) /
2415 *   saturation_degree / max_u_F_prime / 50;
2416 *   else
2417 *   time_step = end_time - time;
2418 *   }
2419 *  
2420 *  
2421 *  
2422 * @endcode
2423 *
2424 * ...and then also update the length of the macro time steps we use while
2425 * we're dealing with time step sizes. In particular, this involves: (i)
2426 * If we have just recomputed the Darcy solution, then the length of the
2427 * previous macro time step is now fixed and the length of the current
2428 * macro time step is, up to now, simply the length of the current (micro)
2429 * time step. (ii) If we have not recomputed the Darcy solution, then the
2430 * length of the current macro time step has just grown by
2431 * <code>time_step</code>.
2432 *
2433 * @code
2434 *   if (solve_for_pressure_and_velocity == true)
2435 *   {
2436 *   old_macro_time_step = current_macro_time_step;
2437 *   current_macro_time_step = time_step;
2438 *   }
2439 *   else
2440 *   current_macro_time_step += time_step;
2441 *  
2442 * @endcode
2443 *
2444 * The last step in this function is to recompute the saturation solution
2445 * based on the velocity field we've just obtained. This naturally happens
2446 * in every time step, and we don't skip any of these computations. At the
2447 * end of computing the saturation, we project back into the allowed
2448 * interval @f$[0,1]@f$ to make sure our solution remains physical.
2449 *
2450 * @code
2451 *   {
2452 *   std::cout << " Solving saturation transport equation..." << std::endl;
2453 *  
2454 *   assemble_saturation_system();
2455 *  
2456 *   SolverControl solver_control(saturation_matrix.m(),
2457 *   1e-16 * saturation_rhs.l2_norm());
2458 *   SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
2459 *  
2460 *   #ifdef DEAL_II_TRILINOS_WITH_EPETRA
2461 *   PreconditionType preconditioner;
2462 *   #else
2463 *   PreconditionType preconditioner("FAST_IC");
2464 *   #endif
2465 *   preconditioner.initialize(saturation_matrix);
2466 *  
2467 *   cg.solve(saturation_matrix,
2468 *   saturation_solution,
2469 *   saturation_rhs,
2470 *   preconditioner);
2471 *  
2472 *   saturation_constraints.distribute(saturation_solution);
2473 *   project_back_saturation();
2474 *  
2475 *   std::cout << " ..." << solver_control.last_step()
2476 *   << " CG iterations." << std::endl;
2477 *   }
2478 *   }
2479 *  
2480 *  
2481 * @endcode
2482 *
2483 *
2484 * <a name="step_43-TwoPhaseFlowProblemdimrefine_mesh"></a>
2485 * <h3>TwoPhaseFlowProblem<dim>::refine_mesh</h3>
2486 *
2487
2488 *
2489 * The next function does the refinement and coarsening of the mesh. It does
2490 * its work in three blocks: (i) Compute refinement indicators by looking at
2491 * the gradient of a solution vector extrapolated linearly from the previous
2492 * two using the respective sizes of the time step (or taking the only
2493 * solution we have if this is the first time step). (ii) Flagging those
2494 * cells for refinement and coarsening where the gradient is larger or
2495 * smaller than a certain threshold, preserving minimal and maximal levels
2496 * of mesh refinement. (iii) Transferring the solution from the old to the
2497 * new mesh. None of this is particularly difficult.
2498 *
2499 * @code
2500 *   template <int dim>
2501 *   void TwoPhaseFlowProblem<dim>::refine_mesh(const unsigned int min_grid_level,
2502 *   const unsigned int max_grid_level)
2503 *   {
2504 *   Vector<double> refinement_indicators(triangulation.n_active_cells());
2505 *   {
2506 *   const QMidpoint<dim> quadrature_formula;
2507 *   FEValues<dim> fe_values(saturation_fe,
2508 *   quadrature_formula,
2509 *   update_gradients);
2510 *   std::vector<Tensor<1, dim>> grad_saturation(1);
2511 *  
2512 *   TrilinosWrappers::MPI::Vector extrapolated_saturation_solution(
2513 *   saturation_solution);
2514 *   if (timestep_number != 0)
2515 *   extrapolated_saturation_solution.sadd((1. + time_step / old_time_step),
2516 *   time_step / old_time_step,
2517 *   old_saturation_solution);
2518 *  
2519 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2520 *   {
2521 *   fe_values.reinit(cell);
2522 *   fe_values.get_function_gradients(extrapolated_saturation_solution,
2523 *   grad_saturation);
2524 *  
2525 *   refinement_indicators(cell->active_cell_index()) =
2526 *   grad_saturation[0].norm();
2527 *   }
2528 *   }
2529 *  
2530 *   {
2531 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2532 *   {
2533 *   cell->clear_coarsen_flag();
2534 *   cell->clear_refine_flag();
2535 *  
2536 *   if ((static_cast<unsigned int>(cell->level()) < max_grid_level) &&
2537 *   (std::fabs(refinement_indicators(cell->active_cell_index())) >
2538 *   saturation_refinement_threshold))
2539 *   cell->set_refine_flag();
2540 *   else if ((static_cast<unsigned int>(cell->level()) >
2541 *   min_grid_level) &&
2542 *   (std::fabs(
2543 *   refinement_indicators(cell->active_cell_index())) <
2544 *   0.5 * saturation_refinement_threshold))
2545 *   cell->set_coarsen_flag();
2546 *   }
2547 *   }
2548 *  
2549 *   triangulation.prepare_coarsening_and_refinement();
2550 *  
2551 *   {
2552 *   std::vector<TrilinosWrappers::MPI::Vector> x_saturation(3);
2553 *   x_saturation[0] = saturation_solution;
2554 *   x_saturation[1] = old_saturation_solution;
2555 *   x_saturation[2] = saturation_matching_last_computed_darcy_solution;
2556 *  
2557 *   std::vector<TrilinosWrappers::MPI::BlockVector> x_darcy(2);
2558 *   x_darcy[0] = last_computed_darcy_solution;
2559 *   x_darcy[1] = second_last_computed_darcy_solution;
2560 *  
2561 *   SolutionTransfer<dim, TrilinosWrappers::MPI::Vector> saturation_soltrans(
2562 *   saturation_dof_handler);
2563 *  
2564 *   SolutionTransfer<dim, TrilinosWrappers::MPI::BlockVector> darcy_soltrans(
2565 *   darcy_dof_handler);
2566 *  
2567 *  
2568 *   triangulation.prepare_coarsening_and_refinement();
2569 *   saturation_soltrans.prepare_for_coarsening_and_refinement(x_saturation);
2570 *  
2571 *   darcy_soltrans.prepare_for_coarsening_and_refinement(x_darcy);
2572 *  
2573 *   triangulation.execute_coarsening_and_refinement();
2574 *   setup_dofs();
2575 *  
2576 *   std::vector<TrilinosWrappers::MPI::Vector> tmp_saturation(3);
2577 *   tmp_saturation[0].reinit(saturation_solution);
2578 *   tmp_saturation[1].reinit(saturation_solution);
2579 *   tmp_saturation[2].reinit(saturation_solution);
2580 *   saturation_soltrans.interpolate(tmp_saturation);
2581 *  
2582 *   saturation_solution = tmp_saturation[0];
2583 *   old_saturation_solution = tmp_saturation[1];
2584 *   saturation_matching_last_computed_darcy_solution = tmp_saturation[2];
2585 *  
2586 *   saturation_constraints.distribute(saturation_solution);
2587 *   saturation_constraints.distribute(old_saturation_solution);
2588 *   saturation_constraints.distribute(
2589 *   saturation_matching_last_computed_darcy_solution);
2590 *  
2591 *   std::vector<TrilinosWrappers::MPI::BlockVector> tmp_darcy(2);
2592 *   tmp_darcy[0].reinit(darcy_solution);
2593 *   tmp_darcy[1].reinit(darcy_solution);
2594 *   darcy_soltrans.interpolate(tmp_darcy);
2595 *  
2596 *   last_computed_darcy_solution = tmp_darcy[0];
2597 *   second_last_computed_darcy_solution = tmp_darcy[1];
2598 *  
2599 *   darcy_constraints.distribute(last_computed_darcy_solution);
2600 *   darcy_constraints.distribute(second_last_computed_darcy_solution);
2601 *  
2602 *   rebuild_saturation_matrix = true;
2603 *   }
2604 *   }
2605 *  
2606 *  
2607 *  
2608 * @endcode
2609 *
2610 *
2611 * <a name="step_43-TwoPhaseFlowProblemdimoutput_results"></a>
2612 * <h3>TwoPhaseFlowProblem<dim>::output_results</h3>
2613 *
2614
2615 *
2616 * This function generates graphical output. It is in essence a copy of the
2617 * implementation in @ref step_31 "step-31".
2618 *
2619 * @code
2620 *   template <int dim>
2621 *   void TwoPhaseFlowProblem<dim>::output_results() const
2622 *   {
2623 *   const FESystem<dim> joint_fe(darcy_fe, 1, saturation_fe, 1);
2624 *   DoFHandler<dim> joint_dof_handler(triangulation);
2625 *   joint_dof_handler.distribute_dofs(joint_fe);
2626 *   Assert(joint_dof_handler.n_dofs() ==
2627 *   darcy_dof_handler.n_dofs() + saturation_dof_handler.n_dofs(),
2628 *   ExcInternalError());
2629 *  
2630 *   Vector<double> joint_solution(joint_dof_handler.n_dofs());
2631 *  
2632 *   {
2633 *   std::vector<types::global_dof_index> local_joint_dof_indices(
2634 *   joint_fe.n_dofs_per_cell());
2635 *   std::vector<types::global_dof_index> local_darcy_dof_indices(
2636 *   darcy_fe.n_dofs_per_cell());
2637 *   std::vector<types::global_dof_index> local_saturation_dof_indices(
2638 *   saturation_fe.n_dofs_per_cell());
2639 *  
2640 *   auto joint_cell = joint_dof_handler.begin_active();
2641 *   const auto joint_endc = joint_dof_handler.end();
2642 *   auto darcy_cell = darcy_dof_handler.begin_active();
2643 *   auto saturation_cell = saturation_dof_handler.begin_active();
2644 *  
2645 *   for (; joint_cell != joint_endc;
2646 *   ++joint_cell, ++darcy_cell, ++saturation_cell)
2647 *   {
2648 *   joint_cell->get_dof_indices(local_joint_dof_indices);
2649 *   darcy_cell->get_dof_indices(local_darcy_dof_indices);
2650 *   saturation_cell->get_dof_indices(local_saturation_dof_indices);
2651 *  
2652 *   for (unsigned int i = 0; i < joint_fe.n_dofs_per_cell(); ++i)
2653 *   if (joint_fe.system_to_base_index(i).first.first == 0)
2654 *   {
2655 *   Assert(joint_fe.system_to_base_index(i).second <
2656 *   local_darcy_dof_indices.size(),
2657 *   ExcInternalError());
2658 *   joint_solution(local_joint_dof_indices[i]) = darcy_solution(
2659 *   local_darcy_dof_indices[joint_fe.system_to_base_index(i)
2660 *   .second]);
2661 *   }
2662 *   else
2663 *   {
2664 *   Assert(joint_fe.system_to_base_index(i).first.first == 1,
2665 *   ExcInternalError());
2666 *   Assert(joint_fe.system_to_base_index(i).second <
2667 *   local_darcy_dof_indices.size(),
2668 *   ExcInternalError());
2669 *   joint_solution(local_joint_dof_indices[i]) =
2670 *   saturation_solution(
2671 *   local_saturation_dof_indices
2672 *   [joint_fe.system_to_base_index(i).second]);
2673 *   }
2674 *   }
2675 *   }
2676 *   std::vector<std::string> joint_solution_names(dim, "velocity");
2677 *   joint_solution_names.emplace_back("pressure");
2678 *   joint_solution_names.emplace_back("saturation");
2679 *  
2680 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
2681 *   data_component_interpretation(
2682 *   dim, DataComponentInterpretation::component_is_part_of_vector);
2683 *   data_component_interpretation.push_back(
2684 *   DataComponentInterpretation::component_is_scalar);
2685 *   data_component_interpretation.push_back(
2686 *   DataComponentInterpretation::component_is_scalar);
2687 *  
2688 *   DataOut<dim> data_out;
2689 *  
2690 *   data_out.attach_dof_handler(joint_dof_handler);
2691 *   data_out.add_data_vector(joint_solution,
2692 *   joint_solution_names,
2693 *   DataOut<dim>::type_dof_data,
2694 *   data_component_interpretation);
2695 *  
2696 *   data_out.build_patches();
2697 *  
2698 *   std::string filename =
2699 *   "solution-" + Utilities::int_to_string(timestep_number, 5) + ".vtu";
2700 *   std::ofstream output(filename);
2701 *   data_out.write_vtu(output);
2702 *   }
2703 *  
2704 *  
2705 *  
2706 * @endcode
2707 *
2708 *
2709 * <a name="step_43-Toolfunctions"></a>
2710 * <h3>Tool functions</h3>
2711 *
2712
2713 *
2714 *
2715 * <a name="step_43-TwoPhaseFlowProblemdimdetermine_whether_to_solve_for_pressure_and_velocity"></a>
2716 * <h4>TwoPhaseFlowProblem<dim>::determine_whether_to_solve_for_pressure_and_velocity</h4>
2717 *
2718
2719 *
2720 * This function implements the a posteriori criterion for adaptive operator
2721 * splitting. The function is relatively straightforward given the way we
2722 * have implemented other functions above and given the formula for the
2723 * criterion derived in the paper.
2724 *
2725
2726 *
2727 * If one decides that one wants the original IMPES method in which the
2728 * Darcy equation is solved in every time step, then this can be achieved by
2729 * setting the threshold value <code>AOS_threshold</code> (with a default of
2730 * @f$5.0@f$) to zero, thereby forcing the function to always return true.
2731 *
2732
2733 *
2734 * Finally, note that the function returns true unconditionally for the
2735 * first two time steps to ensure that we have always solved the Darcy
2736 * system at least twice when skipping its solution, thereby allowing us to
2737 * extrapolate the velocity from the last two solutions in
2738 * <code>solve()</code>.
2739 *
2740 * @code
2741 *   template <int dim>
2742 *   bool TwoPhaseFlowProblem<
2743 *   dim>::determine_whether_to_solve_for_pressure_and_velocity() const
2744 *   {
2745 *   if (timestep_number <= 2)
2746 *   return true;
2747 *  
2748 *   const QGauss<dim> quadrature_formula(saturation_degree + 2);
2749 *   const unsigned int n_q_points = quadrature_formula.size();
2750 *  
2751 *   FEValues<dim> fe_values(saturation_fe,
2752 *   quadrature_formula,
2753 *   update_values | update_quadrature_points);
2754 *  
2755 *   std::vector<double> old_saturation_after_solving_pressure(n_q_points);
2756 *   std::vector<double> present_saturation(n_q_points);
2757 *  
2758 *   std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
2759 *  
2760 *   double max_global_aop_indicator = 0.0;
2761 *  
2762 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2763 *   {
2764 *   double max_local_mobility_reciprocal_difference = 0.0;
2765 *   double max_local_permeability_inverse_l1_norm = 0.0;
2766 *  
2767 *   fe_values.reinit(cell);
2768 *   fe_values.get_function_values(
2769 *   saturation_matching_last_computed_darcy_solution,
2770 *   old_saturation_after_solving_pressure);
2771 *   fe_values.get_function_values(saturation_solution, present_saturation);
2772 *  
2773 *   k_inverse.value_list(fe_values.get_quadrature_points(),
2774 *   k_inverse_values);
2775 *  
2776 *   for (unsigned int q = 0; q < n_q_points; ++q)
2777 *   {
2778 *   const double mobility_reciprocal_difference = std::fabs(
2779 *   mobility_inverse(present_saturation[q], viscosity) -
2780 *   mobility_inverse(old_saturation_after_solving_pressure[q],
2781 *   viscosity));
2782 *  
2783 *   max_local_mobility_reciprocal_difference =
2784 *   std::max(max_local_mobility_reciprocal_difference,
2785 *   mobility_reciprocal_difference);
2786 *  
2787 *   max_local_permeability_inverse_l1_norm =
2788 *   std::max(max_local_permeability_inverse_l1_norm,
2789 *   l1_norm(k_inverse_values[q]));
2790 *   }
2791 *  
2792 *   max_global_aop_indicator =
2793 *   std::max(max_global_aop_indicator,
2794 *   (max_local_mobility_reciprocal_difference *
2795 *   max_local_permeability_inverse_l1_norm));
2796 *   }
2797 *  
2798 *   return (max_global_aop_indicator > AOS_threshold);
2799 *   }
2800 *  
2801 *  
2802 *  
2803 * @endcode
2804 *
2805 *
2806 * <a name="step_43-TwoPhaseFlowProblemdimproject_back_saturation"></a>
2807 * <h4>TwoPhaseFlowProblem<dim>::project_back_saturation</h4>
2808 *
2809
2810 *
2811 * The next function simply makes sure that the saturation values always
2812 * remain within the physically reasonable range of @f$[0,1]@f$. While the
2813 * continuous equations guarantee that this is so, the discrete equations
2814 * don't. However, if we allow the discrete solution to escape this range we
2815 * get into trouble because terms like @f$F(S)@f$ and @f$F'(S)@f$ will produce
2816 * unreasonable results (e.g. @f$F'(S)<0@f$ for @f$S<0@f$, which would imply that
2817 * the wetting fluid phase flows <i>against</i> the direction of the bulk
2818 * fluid velocity)). Consequently, at the end of each time step, we simply
2819 * project the saturation field back into the physically reasonable region.
2820 *
2821 * @code
2822 *   template <int dim>
2823 *   void TwoPhaseFlowProblem<dim>::project_back_saturation()
2824 *   {
2825 *   for (unsigned int i = 0; i < saturation_solution.size(); ++i)
2826 *   if (saturation_solution(i) < 0.2)
2827 *   saturation_solution(i) = 0.2;
2828 *   else if (saturation_solution(i) > 1)
2829 *   saturation_solution(i) = 1;
2830 *   }
2831 *  
2832 *  
2833 *  
2834 * @endcode
2835 *
2836 *
2837 * <a name="step_43-TwoPhaseFlowProblemdimget_max_u_F_prime"></a>
2838 * <h4>TwoPhaseFlowProblem<dim>::get_max_u_F_prime</h4>
2839 *
2840
2841 *
2842 * Another simpler helper function: Compute the maximum of the total
2843 * velocity times the derivative of the fraction flow function, i.e.,
2844 * compute @f$\|\mathbf{u} F'(S)\|_{L_\infty(\Omega)}@f$. This term is used in
2845 * both the computation of the time step as well as in normalizing the
2846 * entropy-residual term in the artificial viscosity.
2847 *
2848 * @code
2849 *   template <int dim>
2850 *   double TwoPhaseFlowProblem<dim>::get_max_u_F_prime() const
2851 *   {
2852 *   const QGauss<dim> quadrature_formula(darcy_degree + 2);
2853 *   const unsigned int n_q_points = quadrature_formula.size();
2854 *  
2855 *   FEValues<dim> darcy_fe_values(darcy_fe, quadrature_formula, update_values);
2856 *   FEValues<dim> saturation_fe_values(saturation_fe,
2857 *   quadrature_formula,
2858 *   update_values);
2859 *  
2860 *   std::vector<Vector<double>> darcy_solution_values(n_q_points,
2861 *   Vector<double>(dim + 1));
2862 *   std::vector<double> saturation_values(n_q_points);
2863 *  
2864 *   double max_velocity_times_dF_dS = 0;
2865 *  
2866 *   auto cell = darcy_dof_handler.begin_active();
2867 *   const auto endc = darcy_dof_handler.end();
2868 *   auto saturation_cell = saturation_dof_handler.begin_active();
2869 *   for (; cell != endc; ++cell, ++saturation_cell)
2870 *   {
2871 *   darcy_fe_values.reinit(cell);
2872 *   saturation_fe_values.reinit(saturation_cell);
2873 *  
2874 *   darcy_fe_values.get_function_values(darcy_solution,
2875 *   darcy_solution_values);
2876 *   saturation_fe_values.get_function_values(old_saturation_solution,
2877 *   saturation_values);
2878 *  
2879 *   for (unsigned int q = 0; q < n_q_points; ++q)
2880 *   {
2881 *   Tensor<1, dim> velocity;
2882 *   for (unsigned int i = 0; i < dim; ++i)
2883 *   velocity[i] = darcy_solution_values[q](i);
2884 *  
2885 *   const double dF_dS =
2886 *   fractional_flow_derivative(saturation_values[q], viscosity);
2887 *  
2888 *   max_velocity_times_dF_dS =
2889 *   std::max(max_velocity_times_dF_dS, velocity.norm() * dF_dS);
2890 *   }
2891 *   }
2892 *  
2893 *   return max_velocity_times_dF_dS;
2894 *   }
2895 *  
2896 *  
2897 * @endcode
2898 *
2899 *
2900 * <a name="step_43-TwoPhaseFlowProblemdimget_extrapolated_saturation_range"></a>
2901 * <h4>TwoPhaseFlowProblem<dim>::get_extrapolated_saturation_range</h4>
2902 *
2903
2904 *
2905 * For computing the stabilization term, we need to know the range of the
2906 * saturation variable. Unlike in @ref step_31 "step-31", this range is trivially bounded
2907 * by the interval @f$[0,1]@f$ but we can do a bit better by looping over a
2908 * collection of quadrature points and seeing what the values are there. If
2909 * we can, i.e., if there are at least two timesteps around, we can even
2910 * take the values extrapolated to the next time step.
2911 *
2912
2913 *
2914 * As before, the function is taken with minimal modifications from @ref step_31 "step-31".
2915 *
2916 * @code
2917 *   template <int dim>
2918 *   std::pair<double, double>
2919 *   TwoPhaseFlowProblem<dim>::get_extrapolated_saturation_range() const
2920 *   {
2921 *   const QGauss<dim> quadrature_formula(saturation_degree + 2);
2922 *   const unsigned int n_q_points = quadrature_formula.size();
2923 *  
2924 *   FEValues<dim> fe_values(saturation_fe, quadrature_formula, update_values);
2925 *   std::vector<double> old_saturation_values(n_q_points);
2926 *   std::vector<double> old_old_saturation_values(n_q_points);
2927 *  
2928 *   if (timestep_number != 0)
2929 *   {
2930 *   double min_saturation = std::numeric_limits<double>::max(),
2931 *   max_saturation = std::numeric_limits<double>::lowest();
2932 *  
2933 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2934 *   {
2935 *   fe_values.reinit(cell);
2936 *   fe_values.get_function_values(old_saturation_solution,
2937 *   old_saturation_values);
2938 *   fe_values.get_function_values(old_old_saturation_solution,
2939 *   old_old_saturation_values);
2940 *  
2941 *   for (unsigned int q = 0; q < n_q_points; ++q)
2942 *   {
2943 *   const double saturation =
2944 *   (1. + time_step / old_time_step) * old_saturation_values[q] -
2945 *   time_step / old_time_step * old_old_saturation_values[q];
2946 *  
2947 *   min_saturation = std::min(min_saturation, saturation);
2948 *   max_saturation = std::max(max_saturation, saturation);
2949 *   }
2950 *   }
2951 *  
2952 *   return std::make_pair(min_saturation, max_saturation);
2953 *   }
2954 *   else
2955 *   {
2956 *   double min_saturation = std::numeric_limits<double>::max(),
2957 *   max_saturation = std::numeric_limits<double>::lowest();
2958 *  
2959 *   for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2960 *   {
2961 *   fe_values.reinit(cell);
2962 *   fe_values.get_function_values(old_saturation_solution,
2963 *   old_saturation_values);
2964 *  
2965 *   for (unsigned int q = 0; q < n_q_points; ++q)
2966 *   {
2967 *   const double saturation = old_saturation_values[q];
2968 *  
2969 *   min_saturation = std::min(min_saturation, saturation);
2970 *   max_saturation = std::max(max_saturation, saturation);
2971 *   }
2972 *   }
2973 *  
2974 *   return std::make_pair(min_saturation, max_saturation);
2975 *   }
2976 *   }
2977 *  
2978 *  
2979 *  
2980 * @endcode
2981 *
2982 *
2983 * <a name="step_43-TwoPhaseFlowProblemdimcompute_viscosity"></a>
2984 * <h4>TwoPhaseFlowProblem<dim>::compute_viscosity</h4>
2985 *
2986
2987 *
2988 * The final tool function is used to compute the artificial viscosity on a
2989 * given cell. This isn't particularly complicated if you have the formula
2990 * for it in front of you, and looking at the implementation in @ref step_31 "step-31". The
2991 * major difference to that tutorial program is that the velocity here is
2992 * not simply @f$\mathbf u@f$ but @f$\mathbf u F'(S)@f$ and some of the formulas
2993 * need to be adjusted accordingly.
2994 *
2995 * @code
2996 *   template <int dim>
2997 *   double TwoPhaseFlowProblem<dim>::compute_viscosity(
2998 *   const std::vector<double> &old_saturation,
2999 *   const std::vector<double> &old_old_saturation,
3000 *   const std::vector<Tensor<1, dim>> &old_saturation_grads,
3001 *   const std::vector<Tensor<1, dim>> &old_old_saturation_grads,
3002 *   const std::vector<Vector<double>> &present_darcy_values,
3003 *   const double global_max_u_F_prime,
3004 *   const double global_S_variation,
3005 *   const double cell_diameter) const
3006 *   {
3007 *   const double beta = .4 * dim;
3008 *   const double alpha = 1;
3009 *  
3010 *   if (global_max_u_F_prime == 0)
3011 *   return 5e-3 * cell_diameter;
3012 *  
3013 *   const unsigned int n_q_points = old_saturation.size();
3014 *  
3015 *   double max_residual = 0;
3016 *   double max_velocity_times_dF_dS = 0;
3017 *  
3018 *   const bool use_dF_dS = true;
3019 *  
3020 *   for (unsigned int q = 0; q < n_q_points; ++q)
3021 *   {
3022 *   Tensor<1, dim> u;
3023 *   for (unsigned int d = 0; d < dim; ++d)
3024 *   u[d] = present_darcy_values[q](d);
3025 *  
3026 *   const double dS_dt = porosity *
3027 *   (old_saturation[q] - old_old_saturation[q]) /
3028 *   old_time_step;
3029 *  
3030 *   const double dF_dS = fractional_flow_derivative(
3031 *   (old_saturation[q] + old_old_saturation[q]) / 2.0, viscosity);
3032 *  
3033 *   const double u_grad_S =
3034 *   u * dF_dS * (old_saturation_grads[q] + old_old_saturation_grads[q]) /
3035 *   2.0;
3036 *  
3037 *   const double residual =
3038 *   std::abs((dS_dt + u_grad_S) *
3039 *   std::pow((old_saturation[q] + old_old_saturation[q]) / 2,
3040 *   alpha - 1.));
3041 *  
3042 *   max_residual = std::max(residual, max_residual);
3043 *   max_velocity_times_dF_dS =
3044 *   std::max(std::sqrt(u * u) * (use_dF_dS ? std::max(dF_dS, 1.) : 1),
3045 *   max_velocity_times_dF_dS);
3046 *   }
3047 *  
3048 *   const double c_R = 1.0;
3049 *   const double global_scaling = c_R * porosity *
3050 *   (global_max_u_F_prime)*global_S_variation /
3051 *   std::pow(global_Omega_diameter, alpha - 2.);
3052 *  
3053 *   return (beta *
3054 *   (max_velocity_times_dF_dS)*std::min(cell_diameter,
3055 *   std::pow(cell_diameter, alpha) *
3056 *   max_residual /
3057 *   global_scaling));
3058 *   }
3059 *  
3060 *  
3061 * @endcode
3062 *
3063 *
3064 * <a name="step_43-TwoPhaseFlowProblemdimrun"></a>
3065 * <h3>TwoPhaseFlowProblem<dim>::run</h3>
3066 *
3067
3068 *
3069 * This function is, besides <code>solve()</code>, the primary function of
3070 * this program as it controls the time iteration as well as when the
3071 * solution is written into output files and when to do mesh refinement.
3072 *
3073
3074 *
3075 * With the exception of the startup code that loops back to the beginning
3076 * of the function through the <code>goto start_time_iteration</code> label,
3077 * everything should be relatively straightforward. In any case, it mimics
3078 * the corresponding function in @ref step_31 "step-31".
3079 *
3080 * @code
3081 *   template <int dim>
3082 *   void TwoPhaseFlowProblem<dim>::run()
3083 *   {
3084 *   const unsigned int initial_refinement = (dim == 2 ? 5 : 2);
3085 *   const unsigned int n_pre_refinement_steps = (dim == 2 ? 3 : 2);
3086 *  
3087 *  
3088 *   GridGenerator::hyper_cube(triangulation, 0, 1);
3089 *   triangulation.refine_global(initial_refinement);
3090 *   global_Omega_diameter = GridTools::diameter(triangulation);
3091 *  
3092 *   setup_dofs();
3093 *  
3094 *   unsigned int pre_refinement_step = 0;
3095 *  
3096 *   start_time_iteration:
3097 *  
3098 *   VectorTools::project(saturation_dof_handler,
3099 *   saturation_constraints,
3100 *   QGauss<dim>(saturation_degree + 2),
3101 *   SaturationInitialValues<dim>(),
3102 *   old_saturation_solution);
3103 *  
3104 *   time_step = old_time_step = 0;
3105 *   current_macro_time_step = old_macro_time_step = 0;
3106 *  
3107 *   time = 0;
3108 *  
3109 *   do
3110 *   {
3111 *   std::cout << "Timestep " << timestep_number << ": t=" << time
3112 *   << ", dt=" << time_step << std::endl;
3113 *  
3114 *   solve();
3115 *  
3116 *   std::cout << std::endl;
3117 *  
3118 *   if (timestep_number % 200 == 0)
3119 *   output_results();
3120 *  
3121 *   if (timestep_number % 25 == 0)
3122 *   refine_mesh(initial_refinement,
3123 *   initial_refinement + n_pre_refinement_steps);
3124 *  
3125 *   if ((timestep_number == 0) &&
3126 *   (pre_refinement_step < n_pre_refinement_steps))
3127 *   {
3128 *   ++pre_refinement_step;
3129 *   goto start_time_iteration;
3130 *   }
3131 *  
3132 *   time += time_step;
3133 *   ++timestep_number;
3134 *  
3135 *   old_old_saturation_solution = old_saturation_solution;
3136 *   old_saturation_solution = saturation_solution;
3137 *   }
3138 *   while (time <= end_time);
3139 *   }
3140 *   } // namespace Step43
3141 *  
3142 *  
3143 *  
3144 * @endcode
3145 *
3146 *
3147 * <a name="step_43-Thecodemaincodefunction"></a>
3148 * <h3>The <code>main()</code> function</h3>
3149 *
3150
3151 *
3152 * The main function looks almost the same as in all other programs. The need
3153 * to initialize the MPI subsystem for a program that uses Trilinos -- even
3154 * for programs that do not actually run in parallel -- is explained in
3155 * @ref step_31 "step-31".
3156 *
3157 * @code
3158 *   int main(int argc, char *argv[])
3159 *   {
3160 *   try
3161 *   {
3162 *   using namespace dealii;
3163 *   using namespace Step43;
3164 *  
3165 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(
3166 *   argc, argv, numbers::invalid_unsigned_int);
3167 *  
3168 * @endcode
3169 *
3170 * This program can only be run in serial. Otherwise, throw an exception.
3171 *
3172 * @code
3173 *   AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
3174 *   ExcMessage(
3175 *   "This program can only be run in serial, use ./step-43"));
3176 *  
3177 *   TwoPhaseFlowProblem<2> two_phase_flow_problem(1);
3178 *   two_phase_flow_problem.run();
3179 *   }
3180 *   catch (std::exception &exc)
3181 *   {
3182 *   std::cerr << std::endl
3183 *   << std::endl
3184 *   << "----------------------------------------------------"
3185 *   << std::endl;
3186 *   std::cerr << "Exception on processing: " << std::endl
3187 *   << exc.what() << std::endl
3188 *   << "Aborting!" << std::endl
3189 *   << "----------------------------------------------------"
3190 *   << std::endl;
3191 *  
3192 *   return 1;
3193 *   }
3194 *   catch (...)
3195 *   {
3196 *   std::cerr << std::endl
3197 *   << std::endl
3198 *   << "----------------------------------------------------"
3199 *   << std::endl;
3200 *   std::cerr << "Unknown exception!" << std::endl
3201 *   << "Aborting!" << std::endl
3202 *   << "----------------------------------------------------"
3203 *   << std::endl;
3204 *   return 1;
3205 *   }
3206 *  
3207 *   return 0;
3208 *   }
3209 * @endcode
3210@anchor step_43-ResultsSection
3211<a name="step_43-Results"></a><h1>Results</h1>
3212
3213
3214
3215The output of this program is not really much different from that of
3216@ref step_21 "step-21": it solves the same problem, after all. Of more importance are
3217quantitative metrics such as the accuracy of the solution as well as
3218the time needed to compute it. These are documented in detail in the
3219two publications listed at the top of this page and we won't repeat
3220them here.
3221
3222That said, no tutorial program is complete without a couple of good
3223pictures, so here is some output of a run in 3d:
3224
3225<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
3226 <tr>
3227 <td align="center">
3228 <img src="https://dealii.org/images/steps/developer/step-43.3d.velocity.png" alt="">
3229 <p align="center">
3230 Velocity vectors of flow through the porous medium with random
3231 permeability model. Streaming paths of high permeability and resulting
3232 high velocity are clearly visible.
3233 </p>
3234 </td>
3235 <td align="center">
3236 <img src="https://dealii.org/images/steps/developer/step-43.3d.streamlines.png" alt="">
3237 <p align="center">
3238 Streamlines colored by the saturation along the streamline path. Blue
3239 streamlines indicate low saturations, i.e., the flow along these
3240 streamlines must be slow or else more fluid would have been
3241 transported along them. On the other hand, green paths indicate high
3242 velocities since the fluid front has already reached further into the
3243 domain.
3244 </p>
3245 </td>
3246 </tr>
3247 <tr>
3248 <td align="center">
3249 <img src="https://dealii.org/images/steps/developer/step-43.3d.saturation.png" alt="">
3250 <p align="center">
3251 Streamlines with a volume rendering of the saturation, showing how far
3252 the fluid front has advanced at this time.
3253 </p>
3254 </td>
3255 <td align="center">
3256 <img src="https://dealii.org/images/steps/developer/step-43.3d.mesh.png" alt="">
3257 <p align="center">
3258 Surface of the mesh showing the adaptive refinement along the front.
3259 </p>
3260 </td>
3261 </tr>
3262</table>
3263
3264
3265<a name="step-43-extensions"></a>
3266<a name="step_43-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
3267
3268
3269The primary objection one may have to this program is that it is still too
3270slow: 3d computations on reasonably fine meshes are simply too expensive to be
3271done routinely and with reasonably quick turn-around. This is similar to the
3272situation we were in when we wrote @ref step_31 "step-31", from which this program has taken
3273much inspiration. The solution is similar as it was there as well: We need to
3274parallelize the program in a way similar to how we derived @ref step_32 "step-32" out of
3275@ref step_31 "step-31". In fact, all of the techniques used in @ref step_32 "step-32" would be transferable
3276to this program as well, making the program run on dozens or hundreds of
3277processors immediately.
3278
3279A different direction is to make the program more relevant to many other
3280porous media applications. Specifically, one avenue is to go to the primary
3281user of porous media flow simulators, namely the oil industry. There,
3282applications in this area are dominated by multiphase flow (i.e., more than
3283the two phases we have here), and the reactions they may have with each other
3284(or any other way phases may exchange mass, such as through dissolution in and
3285bubbling out of gas from the oil phase). Furthermore, the presence of gas
3286often leads to compressibility effects of the fluid. Jointly, these effects
3287are typically formulated in the widely-used "black oil model". True reactions
3288between multiple phases also play a role in oil reservoir modeling when
3289considering controlled burns of oil in the reservoir to raise pressure and
3290temperature. These are much more complex problems, though, and left for future
3291projects.
3292
3293Finally, from a mathematical perspective, we have derived the
3294criterion for re-computing the velocity/pressure solution at a given
3295time step under the assumption that we want to compare the solution we
3296would get at the current time step with that computed the last time we
3297actually solved this system. However, in the program, whenever we did
3298not re-compute the solution, we didn't just use the previously
3299computed solution but instead extrapolated from the previous two times
3300we solved the system. Consequently, the criterion was pessimistically
3301stated: what we should really compare is the solution we would get at
3302the current time step with the extrapolated one. Re-stating the
3303theorem in this regard is left as an exercise.
3304
3305There are also other ways to extend the mathematical foundation of
3306this program; for example, one may say that it isn't the velocity we
3307care about, but in fact the saturation. Thus, one may ask whether the
3308criterion we use here to decide whether @f$\mathbf u@f$ needs to be
3309recomputed is appropriate; one may, for example, suggest that it is
3310also important to decide whether (and by how much) a wrong velocity
3311field in fact affects the solution of the saturation equation. This
3312would then naturally lead to a sensitivity analysis.
3313
3314From an algorithmic viewpoint, we have here used a criterion for refinement
3315that is often used in engineering, namely by looking at the gradient of
3316the solution. However, if you inspect the solution, you will find that
3317it quickly leads to refinement almost everywhere, even in regions where it
3318is clearly not necessary: frequently used therefore does not need to imply
3319that it is a useful criterion to begin with. On the other hand, replacing
3320this criterion by a different and better one should not be very difficult.
3321For example, the KellyErrorEstimator class used in many other programs
3322should certainly be applicable to the current problem as well.
3323 *
3324 *
3325<a name="step_43-PlainProg"></a>
3326<h1> The plain program</h1>
3327@include "step-43.cc"
3328*/
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10398
begin
*  *  iterator begin()
Definition ALE_Finite_Strain_Plasticity.h:10397
 
*  *  *  struct InterferenceTaperTransform *  
Definition ALE_Finite_Strain_Plasticity.h:7883
AffineConstraints::condense
void condense(SparsityPattern &sparsity) const
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Function::vector_value
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Point
Definition point.h:111
TensorFunction::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
second
Point< 2 > second
Definition grid_out.cc:4640
first
Point< 2 > first
Definition grid_out.cc:4639
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
MeshWorker::loop
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:562
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition dof_tools_constraints.cc:3094
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition dof_tools_sparsity.cc:55
DoFTools::make_zero_boundary_constraints
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
Definition dof_tools_constraints.cc:4810
complete_index_set
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1187
data
std::vector< index_type > data
Definition mpi.cc:734
CGALWrappers::K
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
Definition intersections.cc:62
DoFRenumbering::component_wise
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition dof_renumbering.cc:670
DoFRenumbering::Cuthill_McKee
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Definition dof_renumbering.cc:370
DoFRenumbering::random
void random(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_renumbering.cc:2268
DoFTools::count_dofs_per_fe_block
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition dof_tools.cc:2366
FETools::extrapolate
void extrapolate(const DoFHandler< dim, spacedim > &dof1, const InVector &z1, const DoFHandler< dim, spacedim > &dof2, OutVector &z2)
GridTools::minimal_cell_diameter
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< spacedim >()))
Definition grid_tools_geometry.cc:406
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria)
Definition grid_tools_geometry.cc:98
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:56
LAPACKSupport::N
constexpr char N
Definition lapack_support.h:155
LAPACKSupport::zero
constexpr types::blas_int zero
Definition lapack_support.h:179
LAPACKSupport::one
constexpr types::blas_int one
Definition lapack_support.h:183
MatrixTools::apply_boundary_values
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Definition matrix_tools.cc:63
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:210
Patterns::internal::escape
std::string escape(const std::string &input, const PatternBase::OutputStyle style)
Definition patterns.cc:49
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
PlasticityLab::component
*  *  *  ScaleZFunction< dim, Number, components >::ScaleZFunction *  component(component)
Definition ALE_Finite_Strain_Plasticity.h:8885
WorkStream::internal::tbb_no_coloring::run
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition work_stream.h:517
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition p4est_wrappers.cc:50
std::exp
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6465
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6225
Generated by doxygen 1.9.8