Reference documentation for deal.II version 9.6.0
\(\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-21.h
Go to the documentation of this file.
1
808,
809 *   const unsigned int /*component*/ = 0) const override
810 *   {
811 *   return 0;
812 *   }
813 *   };
814 *  
815 *  
816 *  
817 * @endcode
818 *
819 *
820 * <a name="step_21-Pressureboundaryvalues"></a>
821 * <h4>Pressure boundary values</h4>
822 *
823
824 *
825 * The next are pressure boundary values. As mentioned in the introduction,
826 * we choose a linear pressure field:
827 *
828 * @code
829 *   template <int dim>
830 *   class PressureBoundaryValues : public Function<dim>
831 *   {
832 *   public:
833 *   PressureBoundaryValues()
834 *   : Function<dim>(1)
835 *   {}
836 *  
837 *   virtual double value(const Point<dim> &p,
838 *   const unsigned int /*component*/ = 0) const override
839 *   {
840 *   return 1 - p[0];
841 *   }
842 *   };
843 *  
844 *  
845 *  
846 * @endcode
847 *
848 *
849 * <a name="step_21-Saturationboundaryvalues"></a>
850 * <h4>Saturation boundary values</h4>
851 *
852
853 *
854 * Then we also need boundary values on the inflow portions of the
855 * boundary. The question whether something is an inflow part is decided
856 * when assembling the right hand side, we only have to provide a functional
857 * description of the boundary values. This is as explained in the
858 * introduction:
859 *
860 * @code
861 *   template <int dim>
862 *   class SaturationBoundaryValues : public Function<dim>
863 *   {
864 *   public:
865 *   SaturationBoundaryValues()
866 *   : Function<dim>(1)
867 *   {}
868 *  
869 *   virtual double value(const Point<dim> &p,
870 *   const unsigned int /*component*/ = 0) const override
871 *   {
872 *   if (p[0] == 0)
873 *   return 1;
874 *   else
875 *   return 0;
876 *   }
877 *   };
878 *  
879 *  
880 *  
881 * @endcode
882 *
883 *
884 * <a name="step_21-Initialdata"></a>
885 * <h4>Initial data</h4>
886 *
887
888 *
889 * Finally, we need initial data. In reality, we only need initial data for
890 * the saturation, but we are lazy, so we will later, before the first time
891 * step, simply interpolate the entire solution for the previous time step
892 * from a function that contains all vector components.
893 *
894
895 *
896 * We therefore simply create a function that returns zero in all
897 * components. We do that by simply forward every function to the
898 * Functions::ZeroFunction class. Why not use that right away in the places of
899 * this program where we presently use the <code>InitialValues</code> class?
900 * Because this way it is simpler to later go back and choose a different
901 * function for initial values.
902 *
903 * @code
904 *   template <int dim>
905 *   class InitialValues : public Function<dim>
906 *   {
907 *   public:
908 *   InitialValues()
909 *   : Function<dim>(dim + 2)
910 *   {}
911 *  
912 *   virtual double value(const Point<dim> &p,
913 *   const unsigned int component = 0) const override
914 *   {
915 *   return Functions::ZeroFunction<dim>(dim + 2).value(p, component);
916 *   }
917 *  
918 *   virtual void vector_value(const Point<dim> &p,
919 *   Vector<double> &values) const override
920 *   {
921 *   Functions::ZeroFunction<dim>(dim + 2).vector_value(p, values);
922 *   }
923 *   };
924 *  
925 *  
926 *  
927 * @endcode
928 *
929 *
930 * <a name="step_21-Theinversepermeabilitytensor"></a>
931 * <h3>The inverse permeability tensor</h3>
932 *
933
934 *
935 * As announced in the introduction, we implement two different permeability
936 * tensor fields. Each of them we put into a namespace of its own, so that
937 * it will be easy later to replace use of one by the other in the code.
938 *
939
940 *
941 *
942 * <a name="step_21-Singlecurvingcrackpermeability"></a>
943 * <h4>Single curving crack permeability</h4>
944 *
945
946 *
947 * The first function for the permeability was the one that models a single
948 * curving crack. It was already used at the end of @ref step_20 "step-20", and its
949 * functional form is given in the introduction of the present tutorial
950 * program. As in some previous programs, we have to declare a (seemingly
951 * unnecessary) default constructor of the KInverse class to avoid warnings
952 * from some compilers:
953 *
954 * @code
955 *   namespace SingleCurvingCrack
956 *   {
957 *   template <int dim>
958 *   class KInverse : public TensorFunction<2, dim>
959 *   {
960 *   public:
961 *   KInverse()
963 *   {}
964 *  
965 *   virtual void
966 *   value_list(const std::vector<Point<dim>> &points,
967 *   std::vector<Tensor<2, dim>> &values) const override
968 *   {
969 *   AssertDimension(points.size(), values.size());
970 *  
971 *   for (unsigned int p = 0; p < points.size(); ++p)
972 *   {
973 *   values[p].clear();
974 *  
975 *   const double distance_to_flowline =
976 *   std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
977 *  
978 *   const double permeability =
979 *   std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
980 *   (0.1 * 0.1)),
981 *   0.01);
982 *  
983 *   for (unsigned int d = 0; d < dim; ++d)
984 *   values[p][d][d] = 1. / permeability;
985 *   }
986 *   }
987 *   };
988 *   } // namespace SingleCurvingCrack
989 *  
990 *  
991 * @endcode
992 *
993 *
994 * <a name="step_21-Randommediumpermeability"></a>
995 * <h4>Random medium permeability</h4>
996 *
997
998 *
999 * This function does as announced in the introduction, i.e. it creates an
1000 * overlay of exponentials at random places. There is one thing worth
1001 * considering for this class. The issue centers around the problem that the
1002 * class creates the centers of the exponentials using a random function. If
1003 * we therefore created the centers each time we create an object of the
1004 * present type, we would get a different list of centers each time. That's
1005 * not what we expect from classes of this type: they should reliably
1006 * represent the same function.
1007 *
1008
1009 *
1010 * The solution to this problem is to make the list of centers a static
1011 * member variable of this class, i.e. there exists exactly one such
1012 * variable for the entire program, rather than for each object of this
1013 * type. That's exactly what we are going to do.
1014 *
1015
1016 *
1017 * The next problem, however, is that we need a way to initialize this
1018 * variable. Since this variable is initialized at the beginning of the
1019 * program, we can't use a regular member function for that since there may
1020 * not be an object of this type around at the time. The C++ standard
1021 * therefore says that only non-member and static member functions can be
1022 * used to initialize a static variable. We use the latter possibility by
1023 * defining a function <code>get_centers</code> that computes the list of
1024 * center points when called.
1025 *
1026
1027 *
1028 * Note that this class works just fine in both 2d and 3d, with the only
1029 * difference being that we use more points in 3d: by experimenting we find
1030 * that we need more exponentials in 3d than in 2d (we have more ground to
1031 * cover, after all, if we want to keep the distance between centers roughly
1032 * equal), so we choose 40 in 2d and 100 in 3d. For any other dimension, the
1033 * function does presently not know what to do so simply throws an exception
1034 * indicating exactly this.
1035 *
1036 * @code
1037 *   namespace RandomMedium
1038 *   {
1039 *   template <int dim>
1040 *   class KInverse : public TensorFunction<2, dim>
1041 *   {
1042 *   public:
1043 *   KInverse()
1044 *   : TensorFunction<2, dim>()
1045 *   {}
1046 *  
1047 *   virtual void
1048 *   value_list(const std::vector<Point<dim>> &points,
1049 *   std::vector<Tensor<2, dim>> &values) const override
1050 *   {
1051 *   AssertDimension(points.size(), values.size());
1052 *  
1053 *   for (unsigned int p = 0; p < points.size(); ++p)
1054 *   {
1055 *   values[p].clear();
1056 *  
1057 *   double permeability = 0;
1058 *   for (unsigned int i = 0; i < centers.size(); ++i)
1059 *   permeability += std::exp(-(points[p] - centers[i]).norm_square() /
1060 *   (0.05 * 0.05));
1061 *  
1062 *   const double normalized_permeability =
1063 *   std::min(std::max(permeability, 0.01), 4.);
1064 *  
1065 *   for (unsigned int d = 0; d < dim; ++d)
1066 *   values[p][d][d] = 1. / normalized_permeability;
1067 *   }
1068 *   }
1069 *  
1070 *   private:
1071 *   static std::vector<Point<dim>> centers;
1072 *  
1073 *   static std::vector<Point<dim>> get_centers()
1074 *   {
1075 *   const unsigned int N =
1076 *   (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
1077 *  
1078 *   std::vector<Point<dim>> centers_list(N);
1079 *   for (unsigned int i = 0; i < N; ++i)
1080 *   for (unsigned int d = 0; d < dim; ++d)
1081 *   centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
1082 *  
1083 *   return centers_list;
1084 *   }
1085 *   };
1086 *  
1087 *  
1088 *  
1089 *   template <int dim>
1090 *   std::vector<Point<dim>> KInverse<dim>::centers =
1091 *   KInverse<dim>::get_centers();
1092 *   } // namespace RandomMedium
1093 *  
1094 *  
1095 *  
1096 * @endcode
1097 *
1098 *
1099 * <a name="step_21-Theinversemobilityandsaturationfunctions"></a>
1100 * <h3>The inverse mobility and saturation functions</h3>
1101 *
1102
1103 *
1104 * There are two more pieces of data that we need to describe, namely the
1105 * inverse mobility function and the saturation curve. Their form is also
1106 * given in the introduction:
1107 *
1108 * @code
1109 *   double mobility_inverse(const double S, const double viscosity)
1110 *   {
1111 *   return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
1112 *   }
1113 *  
1114 *   double fractional_flow(const double S, const double viscosity)
1115 *   {
1116 *   return S * S / (S * S + viscosity * (1 - S) * (1 - S));
1117 *   }
1118 *  
1119 *  
1120 *  
1121 * @endcode
1122 *
1123 *
1124 * <a name="step_21-Linearsolversandpreconditioners"></a>
1125 * <h3>Linear solvers and preconditioners</h3>
1126 *
1127
1128 *
1129 * The linear solvers we use are also completely analogous to the ones used
1130 * in @ref step_20 "step-20". The following classes are therefore copied verbatim from
1131 * there. Note that the classes here are not only copied from
1132 * @ref step_20 "step-20", but also duplicate classes in deal.II. In a future version of this
1133 * example, they should be replaced by an efficient method, though. There is a
1134 * single change: if the size of a linear system is small, i.e. when the mesh
1135 * is very coarse, then it is sometimes not sufficient to set a maximum of
1136 * <code>src.size()</code> CG iterations before the solver in the
1137 * <code>vmult()</code> function converges. (This is, of course, a result of
1138 * numerical round-off, since we know that on paper, the CG method converges
1139 * in at most <code>src.size()</code> steps.) As a consequence, we set the
1140 * maximum number of iterations equal to the maximum of the size of the linear
1141 * system and 200.
1142 *
1143 * @code
1144 *   template <class MatrixType>
1145 *   class InverseMatrix : public Subscriptor
1146 *   {
1147 *   public:
1148 *   InverseMatrix(const MatrixType &m)
1149 *   : matrix(&m)
1150 *   {}
1151 *  
1152 *   void vmult(Vector<double> &dst, const Vector<double> &src) const
1153 *   {
1154 *   SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
1155 *   1e-8 * src.l2_norm());
1156 *   SolverCG<Vector<double>> cg(solver_control);
1157 *  
1158 *   dst = 0;
1159 *  
1160 *   cg.solve(*matrix, dst, src, PreconditionIdentity());
1161 *   }
1162 *  
1163 *   private:
1164 *   const SmartPointer<const MatrixType> matrix;
1165 *   };
1166 *  
1167 *  
1168 *  
1169 *   class SchurComplement : public Subscriptor
1170 *   {
1171 *   public:
1172 *   SchurComplement(const BlockSparseMatrix<double> &A,
1173 *   const InverseMatrix<SparseMatrix<double>> &Minv)
1174 *   : system_matrix(&A)
1175 *   , m_inverse(&Minv)
1176 *   , tmp1(A.block(0, 0).m())
1177 *   , tmp2(A.block(0, 0).m())
1178 *   {}
1179 *  
1180 *   void vmult(Vector<double> &dst, const Vector<double> &src) const
1181 *   {
1182 *   system_matrix->block(0, 1).vmult(tmp1, src);
1183 *   m_inverse->vmult(tmp2, tmp1);
1184 *   system_matrix->block(1, 0).vmult(dst, tmp2);
1185 *   }
1186 *  
1187 *   private:
1188 *   const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
1189 *   const SmartPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
1190 *  
1191 *   mutable Vector<double> tmp1, tmp2;
1192 *   };
1193 *  
1194 *  
1195 *  
1196 *   class ApproximateSchurComplement : public Subscriptor
1197 *   {
1198 *   public:
1199 *   ApproximateSchurComplement(const BlockSparseMatrix<double> &A)
1200 *   : system_matrix(&A)
1201 *   , tmp1(A.block(0, 0).m())
1202 *   , tmp2(A.block(0, 0).m())
1203 *   {}
1204 *  
1205 *   void vmult(Vector<double> &dst, const Vector<double> &src) const
1206 *   {
1207 *   system_matrix->block(0, 1).vmult(tmp1, src);
1208 *   system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
1209 *   system_matrix->block(1, 0).vmult(dst, tmp2);
1210 *   }
1211 *  
1212 *   private:
1213 *   const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
1214 *  
1215 *   mutable Vector<double> tmp1, tmp2;
1216 *   };
1217 *  
1218 *  
1219 *  
1220 * @endcode
1221 *
1222 *
1223 * <a name="step_21-codeTwoPhaseFlowProblemcodeclassimplementation"></a>
1224 * <h3><code>TwoPhaseFlowProblem</code> class implementation</h3>
1225 *
1226
1227 *
1228 * Here now the implementation of the main class. Much of it is actually
1229 * copied from @ref step_20 "step-20", so we won't comment on it in much detail. You should
1230 * try to get familiar with that program first, then most of what is
1231 * happening here should be mostly clear.
1232 *
1233
1234 *
1235 *
1236 * <a name="step_21-TwoPhaseFlowProblemTwoPhaseFlowProblem"></a>
1237 * <h4>TwoPhaseFlowProblem::TwoPhaseFlowProblem</h4>
1238 *
1239
1240 *
1241 * First for the constructor. We use @f$RT_k \times DQ_k \times DQ_k@f$
1242 * spaces. For initializing the DiscreteTime object, we don't set the time
1243 * step size in the constructor because we don't have its value yet.
1244 * The time step size is initially set to zero, but it will be computed
1245 * before it is needed to increment time, as described in a subsection of
1246 * the introduction. The time object internally prevents itself from being
1247 * incremented when @f$dt = 0@f$, forcing us to set a non-zero desired size for
1248 * @f$dt@f$ before advancing time.
1249 *
1250 * @code
1251 *   template <int dim>
1252 *   TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem(const unsigned int degree)
1253 *   : degree(degree)
1254 *   , fe(FE_RaviartThomas<dim>(degree),
1255 *   FE_DGQ<dim>(degree),
1256 *   FE_DGQ<dim>(degree))
1257 *   , dof_handler(triangulation)
1258 *   , n_refinement_steps(5)
1259 *   , time(/*start time*/ 0., /*end time*/ 1.)
1260 *   , viscosity(0.2)
1261 *   {}
1262 *  
1263 *  
1264 *  
1265 * @endcode
1266 *
1267 *
1268 * <a name="step_21-TwoPhaseFlowProblemmake_grid_and_dofs"></a>
1269 * <h4>TwoPhaseFlowProblem::make_grid_and_dofs</h4>
1270 *
1271
1272 *
1273 * This next function starts out with well-known functions calls that create
1274 * and refine a mesh, and then associate degrees of freedom with it. It does
1275 * all the same things as in @ref step_20 "step-20", just now for three components instead
1276 * of two.
1277 *
1278 * @code
1279 *   template <int dim>
1280 *   void TwoPhaseFlowProblem<dim>::make_grid_and_dofs()
1281 *   {
1283 *   triangulation.refine_global(n_refinement_steps);
1284 *  
1285 *   dof_handler.distribute_dofs(fe);
1286 *   DoFRenumbering::component_wise(dof_handler);
1287 *  
1288 *   const std::vector<types::global_dof_index> dofs_per_component =
1289 *   DoFTools::count_dofs_per_fe_component(dof_handler);
1290 *   const unsigned int n_u = dofs_per_component[0],
1291 *   n_p = dofs_per_component[dim],
1292 *   n_s = dofs_per_component[dim + 1];
1293 *  
1294 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
1295 *   << std::endl
1296 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
1297 *   << " (" << n_u << '+' << n_p << '+' << n_s << ')' << std::endl
1298 *   << std::endl;
1299 *  
1300 *   const std::vector<types::global_dof_index> block_sizes = {n_u, n_p, n_s};
1301 *   BlockDynamicSparsityPattern dsp(block_sizes, block_sizes);
1302 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
1303 *  
1304 *   sparsity_pattern.copy_from(dsp);
1305 *   system_matrix.reinit(sparsity_pattern);
1306 *  
1307 *   solution.reinit(block_sizes);
1308 *   old_solution.reinit(block_sizes);
1309 *   system_rhs.reinit(block_sizes);
1310 *   }
1311 *  
1312 *  
1313 * @endcode
1314 *
1315 *
1316 * <a name="step_21-TwoPhaseFlowProblemassemble_system"></a>
1317 * <h4>TwoPhaseFlowProblem::assemble_system</h4>
1318 *
1319
1320 *
1321 * This is the function that assembles the linear system, or at least
1322 * everything except the (1,3) block that depends on the still-unknown
1323 * velocity computed during this time step (we deal with this in
1324 * <code>assemble_rhs_S</code>). Much of it is again as in @ref step_20 "step-20", but we
1325 * have to deal with some nonlinearity this time. However, the top of the
1326 * function is pretty much as usual (note that we set matrix and right hand
1327 * side to zero at the beginning &mdash; something we didn't have to do for
1328 * stationary problems since there we use each matrix object only once and
1329 * it is empty at the beginning anyway).
1330 *
1331
1332 *
1333 * Note that in its present form, the function uses the permeability
1334 * implemented in the RandomMedium::KInverse class. Switching to the single
1335 * curved crack permeability function is as simple as just changing the
1336 * namespace name.
1337 *
1338 * @code
1339 *   template <int dim>
1340 *   void TwoPhaseFlowProblem<dim>::assemble_system()
1341 *   {
1342 *   system_matrix = 0;
1343 *   system_rhs = 0;
1344 *  
1345 *   const QGauss<dim> quadrature_formula(degree + 2);
1346 *   const QGauss<dim - 1> face_quadrature_formula(degree + 2);
1347 *  
1348 *   FEValues<dim> fe_values(fe,
1349 *   quadrature_formula,
1352 *   FEFaceValues<dim> fe_face_values(fe,
1353 *   face_quadrature_formula,
1356 *   update_JxW_values);
1357 *  
1358 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1359 *  
1360 *   const unsigned int n_q_points = quadrature_formula.size();
1361 *   const unsigned int n_face_q_points = face_quadrature_formula.size();
1362 *  
1363 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1364 *   Vector<double> local_rhs(dofs_per_cell);
1365 *  
1366 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1367 *  
1368 *   const PressureRightHandSide<dim> pressure_right_hand_side;
1369 *   const PressureBoundaryValues<dim> pressure_boundary_values;
1370 *   const RandomMedium::KInverse<dim> k_inverse;
1371 *  
1372 *   std::vector<double> pressure_rhs_values(n_q_points);
1373 *   std::vector<double> boundary_values(n_face_q_points);
1374 *   std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1375 *  
1376 *   std::vector<Vector<double>> old_solution_values(n_q_points,
1377 *   Vector<double>(dim + 2));
1378 *  
1379 *   const FEValuesExtractors::Vector velocities(0);
1380 *   const FEValuesExtractors::Scalar pressure(dim);
1381 *   const FEValuesExtractors::Scalar saturation(dim + 1);
1382 *  
1383 *   for (const auto &cell : dof_handler.active_cell_iterators())
1384 *   {
1385 *   fe_values.reinit(cell);
1386 *   local_matrix = 0;
1387 *   local_rhs = 0;
1388 *  
1389 * @endcode
1390 *
1391 * Here's the first significant difference: We have to get the values
1392 * of the saturation function of the previous time step at the
1393 * quadrature points. To this end, we can use the
1394 * FEValues::get_function_values (previously already used in @ref step_9 "step-9",
1395 * @ref step_14 "step-14" and @ref step_15 "step-15"), a function that takes a solution vector and
1396 * returns a list of function values at the quadrature points of the
1397 * present cell. In fact, it returns the complete vector-valued
1398 * solution at each quadrature point, i.e. not only the saturation but
1399 * also the velocities and pressure:
1400 *
1401 * @code
1402 *   fe_values.get_function_values(old_solution, old_solution_values);
1403 *  
1404 * @endcode
1405 *
1406 * Then we also have to get the values of the pressure right hand side
1407 * and of the inverse permeability tensor at the quadrature points:
1408 *
1409 * @code
1410 *   pressure_right_hand_side.value_list(fe_values.get_quadrature_points(),
1411 *   pressure_rhs_values);
1412 *   k_inverse.value_list(fe_values.get_quadrature_points(),
1413 *   k_inverse_values);
1414 *  
1415 * @endcode
1416 *
1417 * With all this, we can now loop over all the quadrature points and
1418 * shape functions on this cell and assemble those parts of the matrix
1419 * and right hand side that we deal with in this function. The
1420 * individual terms in the contributions should be self-explanatory
1421 * given the explicit form of the bilinear form stated in the
1422 * introduction:
1423 *
1424 * @code
1425 *   for (unsigned int q = 0; q < n_q_points; ++q)
1426 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1427 *   {
1428 *   const double old_s = old_solution_values[q](dim + 1);
1429 *  
1430 *   const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1431 *   const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1432 *   const double phi_i_p = fe_values[pressure].value(i, q);
1433 *   const double phi_i_s = fe_values[saturation].value(i, q);
1434 *  
1435 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1436 *   {
1437 *   const Tensor<1, dim> phi_j_u =
1438 *   fe_values[velocities].value(j, q);
1439 *   const double div_phi_j_u =
1440 *   fe_values[velocities].divergence(j, q);
1441 *   const double phi_j_p = fe_values[pressure].value(j, q);
1442 *   const double phi_j_s = fe_values[saturation].value(j, q);
1443 *  
1444 *   local_matrix(i, j) +=
1445 *   (phi_i_u * k_inverse_values[q] *
1446 *   mobility_inverse(old_s, viscosity) * phi_j_u -
1447 *   div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u +
1448 *   phi_i_s * phi_j_s) *
1449 *   fe_values.JxW(q);
1450 *   }
1451 *  
1452 *   local_rhs(i) +=
1453 *   (-phi_i_p * pressure_rhs_values[q]) * fe_values.JxW(q);
1454 *   }
1455 *  
1456 *  
1457 * @endcode
1458 *
1459 * Next, we also have to deal with the pressure boundary values. This,
1460 * again is as in @ref step_20 "step-20":
1461 *
1462 * @code
1463 *   for (const auto &face : cell->face_iterators())
1464 *   if (face->at_boundary())
1465 *   {
1466 *   fe_face_values.reinit(cell, face);
1467 *  
1468 *   pressure_boundary_values.value_list(
1469 *   fe_face_values.get_quadrature_points(), boundary_values);
1470 *  
1471 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1472 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1473 *   {
1474 *   const Tensor<1, dim> phi_i_u =
1475 *   fe_face_values[velocities].value(i, q);
1476 *  
1477 *   local_rhs(i) +=
1478 *   -(phi_i_u * fe_face_values.normal_vector(q) *
1479 *   boundary_values[q] * fe_face_values.JxW(q));
1480 *   }
1481 *   }
1482 *  
1483 * @endcode
1484 *
1485 * The final step in the loop over all cells is to transfer local
1486 * contributions into the global matrix and right hand side vector:
1487 *
1488 * @code
1489 *   cell->get_dof_indices(local_dof_indices);
1490 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1491 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1492 *   system_matrix.add(local_dof_indices[i],
1493 *   local_dof_indices[j],
1494 *   local_matrix(i, j));
1495 *  
1496 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1497 *   system_rhs(local_dof_indices[i]) += local_rhs(i);
1498 *   }
1499 *   }
1500 *  
1501 *  
1502 * @endcode
1503 *
1504 * So much for assembly of matrix and right hand side. Note that we do not
1505 * have to interpolate and apply boundary values since they have all been
1506 * taken care of in the weak form already.
1507 *
1508
1509 *
1510 *
1511
1512 *
1513 *
1514 * <a name="step_21-TwoPhaseFlowProblemassemble_rhs_S"></a>
1515 * <h4>TwoPhaseFlowProblem::assemble_rhs_S</h4>
1516 *
1517
1518 *
1519 * As explained in the introduction, we can only evaluate the right hand
1520 * side of the saturation equation once the velocity has been computed. We
1521 * therefore have this separate function to this end.
1522 *
1523 * @code
1524 *   template <int dim>
1525 *   void TwoPhaseFlowProblem<dim>::assemble_rhs_S()
1526 *   {
1527 *   const QGauss<dim> quadrature_formula(degree + 2);
1528 *   const QGauss<dim - 1> face_quadrature_formula(degree + 2);
1529 *   FEValues<dim> fe_values(fe,
1530 *   quadrature_formula,
1531 *   update_values | update_gradients |
1532 *   update_quadrature_points | update_JxW_values);
1533 *   FEFaceValues<dim> fe_face_values(fe,
1534 *   face_quadrature_formula,
1535 *   update_values | update_normal_vectors |
1536 *   update_quadrature_points |
1537 *   update_JxW_values);
1538 *   FEFaceValues<dim> fe_face_values_neighbor(fe,
1539 *   face_quadrature_formula,
1540 *   update_values);
1541 *  
1542 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1543 *   const unsigned int n_q_points = quadrature_formula.size();
1544 *   const unsigned int n_face_q_points = face_quadrature_formula.size();
1545 *  
1546 *   Vector<double> local_rhs(dofs_per_cell);
1547 *  
1548 *   std::vector<Vector<double>> old_solution_values(n_q_points,
1549 *   Vector<double>(dim + 2));
1550 *   std::vector<Vector<double>> old_solution_values_face(n_face_q_points,
1551 *   Vector<double>(dim +
1552 *   2));
1553 *   std::vector<Vector<double>> old_solution_values_face_neighbor(
1554 *   n_face_q_points, Vector<double>(dim + 2));
1555 *   std::vector<Vector<double>> present_solution_values(n_q_points,
1556 *   Vector<double>(dim +
1557 *   2));
1558 *   std::vector<Vector<double>> present_solution_values_face(
1559 *   n_face_q_points, Vector<double>(dim + 2));
1560 *  
1561 *   std::vector<double> neighbor_saturation(n_face_q_points);
1562 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1563 *  
1564 *   SaturationBoundaryValues<dim> saturation_boundary_values;
1565 *  
1566 *   const FEValuesExtractors::Scalar saturation(dim + 1);
1567 *  
1568 *   for (const auto &cell : dof_handler.active_cell_iterators())
1569 *   {
1570 *   local_rhs = 0;
1571 *   fe_values.reinit(cell);
1572 *  
1573 *   fe_values.get_function_values(old_solution, old_solution_values);
1574 *   fe_values.get_function_values(solution, present_solution_values);
1575 *  
1576 * @endcode
1577 *
1578 * First for the cell terms. These are, following the formulas in the
1579 * introduction, @f$(S^n,\sigma)-(F(S^n) \mathbf{v}^{n+1},\nabla
1580 * \sigma)@f$, where @f$\sigma@f$ is the saturation component of the test
1581 * function:
1582 *
1583 * @code
1584 *   for (unsigned int q = 0; q < n_q_points; ++q)
1585 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1586 *   {
1587 *   const double old_s = old_solution_values[q](dim + 1);
1588 *   Tensor<1, dim> present_u;
1589 *   for (unsigned int d = 0; d < dim; ++d)
1590 *   present_u[d] = present_solution_values[q](d);
1591 *  
1592 *   const double phi_i_s = fe_values[saturation].value(i, q);
1593 *   const Tensor<1, dim> grad_phi_i_s =
1594 *   fe_values[saturation].gradient(i, q);
1595 *  
1596 *   local_rhs(i) +=
1597 *   (time.get_next_step_size() * fractional_flow(old_s, viscosity) *
1598 *   present_u * grad_phi_i_s +
1599 *   old_s * phi_i_s) *
1600 *   fe_values.JxW(q);
1601 *   }
1602 *  
1603 * @endcode
1604 *
1605 * Secondly, we have to deal with the flux parts on the face
1606 * boundaries. This was a bit more involved because we first have to
1607 * determine which are the influx and outflux parts of the cell
1608 * boundary. If we have an influx boundary, we need to evaluate the
1609 * saturation on the other side of the face (or the boundary values,
1610 * if we are at the boundary of the domain).
1611 *
1612
1613 *
1614 * All this is a bit tricky, but has been explained in some detail
1615 * already in @ref step_9 "step-9". Take a look there how this is supposed to work!
1616 *
1617 * @code
1618 *   for (const auto face_no : cell->face_indices())
1619 *   {
1620 *   fe_face_values.reinit(cell, face_no);
1621 *  
1622 *   fe_face_values.get_function_values(old_solution,
1623 *   old_solution_values_face);
1624 *   fe_face_values.get_function_values(solution,
1625 *   present_solution_values_face);
1626 *  
1627 *   if (cell->at_boundary(face_no))
1628 *   saturation_boundary_values.value_list(
1629 *   fe_face_values.get_quadrature_points(), neighbor_saturation);
1630 *   else
1631 *   {
1632 *   const auto neighbor = cell->neighbor(face_no);
1633 *   const unsigned int neighbor_face =
1634 *   cell->neighbor_of_neighbor(face_no);
1635 *  
1636 *   fe_face_values_neighbor.reinit(neighbor, neighbor_face);
1637 *  
1638 *   fe_face_values_neighbor.get_function_values(
1639 *   old_solution, old_solution_values_face_neighbor);
1640 *  
1641 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1642 *   neighbor_saturation[q] =
1643 *   old_solution_values_face_neighbor[q](dim + 1);
1644 *   }
1645 *  
1646 *  
1647 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1648 *   {
1649 *   Tensor<1, dim> present_u_face;
1650 *   for (unsigned int d = 0; d < dim; ++d)
1651 *   present_u_face[d] = present_solution_values_face[q](d);
1652 *  
1653 *   const double normal_flux =
1654 *   present_u_face * fe_face_values.normal_vector(q);
1655 *  
1656 *   const bool is_outflow_q_point = (normal_flux >= 0);
1657 *  
1658 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1659 *   local_rhs(i) -=
1660 *   time.get_next_step_size() * normal_flux *
1661 *   fractional_flow((is_outflow_q_point == true ?
1662 *   old_solution_values_face[q](dim + 1) :
1663 *   neighbor_saturation[q]),
1664 *   viscosity) *
1665 *   fe_face_values[saturation].value(i, q) *
1666 *   fe_face_values.JxW(q);
1667 *   }
1668 *   }
1669 *  
1670 *   cell->get_dof_indices(local_dof_indices);
1671 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1672 *   system_rhs(local_dof_indices[i]) += local_rhs(i);
1673 *   }
1674 *   }
1675 *  
1676 *  
1677 *  
1678 * @endcode
1679 *
1680 *
1681 * <a name="step_21-TwoPhaseFlowProblemsolve"></a>
1682 * <h4>TwoPhaseFlowProblem::solve</h4>
1683 *
1684
1685 *
1686 * After all these preparations, we finally solve the linear system for
1687 * velocity and pressure in the same way as in @ref step_20 "step-20". After that, we have
1688 * to deal with the saturation equation (see below):
1689 *
1690 * @code
1691 *   template <int dim>
1692 *   void TwoPhaseFlowProblem<dim>::solve()
1693 *   {
1694 *   const InverseMatrix<SparseMatrix<double>> m_inverse(
1695 *   system_matrix.block(0, 0));
1696 *   Vector<double> tmp(solution.block(0).size());
1697 *   Vector<double> schur_rhs(solution.block(1).size());
1698 *   Vector<double> tmp2(solution.block(2).size());
1699 *  
1700 *  
1701 * @endcode
1702 *
1703 * First the pressure, using the pressure Schur complement of the first
1704 * two equations:
1705 *
1706 * @code
1707 *   {
1708 *   m_inverse.vmult(tmp, system_rhs.block(0));
1709 *   system_matrix.block(1, 0).vmult(schur_rhs, tmp);
1710 *   schur_rhs -= system_rhs.block(1);
1711 *  
1712 *  
1713 *   SchurComplement schur_complement(system_matrix, m_inverse);
1714 *  
1715 *   ApproximateSchurComplement approximate_schur_complement(system_matrix);
1716 *  
1717 *   InverseMatrix<ApproximateSchurComplement> preconditioner(
1718 *   approximate_schur_complement);
1719 *  
1720 *  
1721 *   SolverControl solver_control(solution.block(1).size(),
1722 *   1e-12 * schur_rhs.l2_norm());
1723 *   SolverCG<Vector<double>> cg(solver_control);
1724 *  
1725 *   cg.solve(schur_complement, solution.block(1), schur_rhs, preconditioner);
1726 *  
1727 *   std::cout << " " << solver_control.last_step()
1728 *   << " CG Schur complement iterations for pressure." << std::endl;
1729 *   }
1730 *  
1731 * @endcode
1732 *
1733 * Now the velocity:
1734 *
1735 * @code
1736 *   {
1737 *   system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1738 *   tmp *= -1;
1739 *   tmp += system_rhs.block(0);
1740 *  
1741 *   m_inverse.vmult(solution.block(0), tmp);
1742 *   }
1743 *  
1744 * @endcode
1745 *
1746 * Finally, we have to take care of the saturation equation. The first
1747 * business we have here is to determine the time step using the formula
1748 * in the introduction. Knowing the shape of our domain and that we
1749 * created the mesh by regular subdivision of cells, we can compute the
1750 * diameter of each of our cells quite easily (in fact we use the linear
1751 * extensions in coordinate directions of the cells, not the
1752 * diameter). Note that we will learn a more general way to do this in
1753 * @ref step_24 "step-24", where we use the GridTools::minimal_cell_diameter function.
1754 *
1755
1756 *
1757 * The maximal velocity we compute using a helper function to compute the
1758 * maximal velocity defined below, and with all this we can evaluate our
1759 * new time step length. We use the method
1760 * DiscreteTime::set_desired_next_time_step() to suggest the new
1761 * calculated value of the time step to the DiscreteTime object. In most
1762 * cases, the time object uses the exact provided value to increment time.
1763 * It some case, the step size may be modified further by the time object.
1764 * For example, if the calculated time increment overshoots the end time,
1765 * it is truncated accordingly.
1766 *
1767 * @code
1768 *   time.set_desired_next_step_size(std::pow(0.5, double(n_refinement_steps)) /
1769 *   get_maximal_velocity());
1770 *  
1771 * @endcode
1772 *
1773 * The next step is to assemble the right hand side, and then to pass
1774 * everything on for solution. At the end, we project back saturations
1775 * onto the physically reasonable range:
1776 *
1777 * @code
1778 *   assemble_rhs_S();
1779 *   {
1780 *   SolverControl solver_control(system_matrix.block(2, 2).m(),
1781 *   1e-8 * system_rhs.block(2).l2_norm());
1782 *   SolverCG<Vector<double>> cg(solver_control);
1783 *   cg.solve(system_matrix.block(2, 2),
1784 *   solution.block(2),
1785 *   system_rhs.block(2),
1786 *   PreconditionIdentity());
1787 *  
1788 *   project_back_saturation();
1789 *  
1790 *   std::cout << " " << solver_control.last_step()
1791 *   << " CG iterations for saturation." << std::endl;
1792 *   }
1793 *  
1794 *  
1795 *   old_solution = solution;
1796 *   }
1797 *  
1798 *  
1799 * @endcode
1800 *
1801 *
1802 * <a name="step_21-TwoPhaseFlowProblemoutput_results"></a>
1803 * <h4>TwoPhaseFlowProblem::output_results</h4>
1804 *
1805
1806 *
1807 * There is nothing surprising here. Since the program will do a lot of time
1808 * steps, we create an output file only every fifth time step and skip all
1809 * other time steps at the top of the file already.
1810 *
1811
1812 *
1813 * When creating file names for output close to the bottom of the function,
1814 * we convert the number of the time step to a string representation that
1815 * is padded by leading zeros to four digits. We do this because this way
1816 * all output file names have the same length, and consequently sort well
1817 * when creating a directory listing.
1818 *
1819 * @code
1820 *   template <int dim>
1821 *   void TwoPhaseFlowProblem<dim>::output_results() const
1822 *   {
1823 *   if (time.get_step_number() % 5 != 0)
1824 *   return;
1825 *  
1826 *   std::vector<std::string> solution_names;
1827 *   switch (dim)
1828 *   {
1829 *   case 2:
1830 *   solution_names = {"u", "v", "p", "S"};
1831 *   break;
1832 *  
1833 *   case 3:
1834 *   solution_names = {"u", "v", "w", "p", "S"};
1835 *   break;
1836 *  
1837 *   default:
1838 *   DEAL_II_NOT_IMPLEMENTED();
1839 *   }
1840 *  
1841 *   DataOut<dim> data_out;
1842 *  
1843 *   data_out.attach_dof_handler(dof_handler);
1844 *   data_out.add_data_vector(solution, solution_names);
1845 *  
1846 *   data_out.build_patches(degree + 1);
1847 *  
1848 *   std::ofstream output("solution-" +
1849 *   Utilities::int_to_string(time.get_step_number(), 4) +
1850 *   ".vtk");
1851 *   data_out.write_vtk(output);
1852 *   }
1853 *  
1854 *  
1855 *  
1856 * @endcode
1857 *
1858 *
1859 * <a name="step_21-TwoPhaseFlowProblemproject_back_saturation"></a>
1860 * <h4>TwoPhaseFlowProblem::project_back_saturation</h4>
1861 *
1862
1863 *
1864 * In this function, we simply run over all saturation degrees of freedom
1865 * and make sure that if they should have left the physically reasonable
1866 * range, that they be reset to the interval @f$[0,1]@f$. To do this, we only
1867 * have to loop over all saturation components of the solution vector; these
1868 * are stored in the block 2 (block 0 are the velocities, block 1 are the
1869 * pressures).
1870 *
1871
1872 *
1873 * It may be instructive to note that this function almost never triggers
1874 * when the time step is chosen as mentioned in the introduction. However,
1875 * if we choose the timestep only slightly larger, we get plenty of values
1876 * outside the proper range. Strictly speaking, the function is therefore
1877 * unnecessary if we choose the time step small enough. In a sense, the
1878 * function is therefore only a safety device to avoid situations where our
1879 * entire solution becomes unphysical because individual degrees of freedom
1880 * have become unphysical a few time steps earlier.
1881 *
1882 * @code
1883 *   template <int dim>
1884 *   void TwoPhaseFlowProblem<dim>::project_back_saturation()
1885 *   {
1886 *   for (unsigned int i = 0; i < solution.block(2).size(); ++i)
1887 *   if (solution.block(2)(i) < 0)
1888 *   solution.block(2)(i) = 0;
1889 *   else if (solution.block(2)(i) > 1)
1890 *   solution.block(2)(i) = 1;
1891 *   }
1892 *  
1893 *  
1894 * @endcode
1895 *
1896 *
1897 * <a name="step_21-TwoPhaseFlowProblemget_maximal_velocity"></a>
1898 * <h4>TwoPhaseFlowProblem::get_maximal_velocity</h4>
1899 *
1900
1901 *
1902 * The following function is used in determining the maximal allowable time
1903 * step. What it does is to loop over all quadrature points in the domain
1904 * and find what the maximal magnitude of the velocity is.
1905 *
1906 * @code
1907 *   template <int dim>
1908 *   double TwoPhaseFlowProblem<dim>::get_maximal_velocity() const
1909 *   {
1910 *   const QGauss<dim> quadrature_formula(degree + 2);
1911 *   const unsigned int n_q_points = quadrature_formula.size();
1912 *  
1913 *   FEValues<dim> fe_values(fe, quadrature_formula, update_values);
1914 *   std::vector<Vector<double>> solution_values(n_q_points,
1915 *   Vector<double>(dim + 2));
1916 *   double max_velocity = 0;
1917 *  
1918 *   for (const auto &cell : dof_handler.active_cell_iterators())
1919 *   {
1920 *   fe_values.reinit(cell);
1921 *   fe_values.get_function_values(solution, solution_values);
1922 *  
1923 *   for (unsigned int q = 0; q < n_q_points; ++q)
1924 *   {
1925 *   Tensor<1, dim> velocity;
1926 *   for (unsigned int i = 0; i < dim; ++i)
1927 *   velocity[i] = solution_values[q](i);
1928 *  
1929 *   max_velocity = std::max(max_velocity, velocity.norm());
1930 *   }
1931 *   }
1932 *  
1933 *   return max_velocity;
1934 *   }
1935 *  
1936 *  
1937 * @endcode
1938 *
1939 *
1940 * <a name="step_21-TwoPhaseFlowProblemrun"></a>
1941 * <h4>TwoPhaseFlowProblem::run</h4>
1942 *
1943
1944 *
1945 * This is the final function of our main class. Its brevity speaks for
1946 * itself. There are only two points worth noting: First, the function
1947 * projects the initial values onto the finite element space at the
1948 * beginning; the VectorTools::project function doing this requires an
1949 * argument indicating the hanging node constraints. We have none in this
1950 * program (we compute on a uniformly refined mesh), but the function
1951 * requires the argument anyway, of course. So we have to create a
1952 * constraint object. In its original state, constraint objects are
1953 * unsorted, and have to be sorted (using the AffineConstraints::close
1954 * function) before they can be used. This is what we do here, and which is
1955 * why we can't simply call the VectorTools::project function with an
1956 * anonymous temporary object <code>AffineConstraints<double>()</code> as the
1957 * second argument.
1958 *
1959
1960 *
1961 * The second point worth mentioning is that we only compute the length of
1962 * the present time step in the middle of solving the linear system
1963 * corresponding to each time step. We can therefore output the present
1964 * time of a time step only at the end of the time step.
1965 * We increment time by calling the method DiscreteTime::advance_time()
1966 * inside the loop. Since we are reporting the time and dt after we
1967 * increment it, we have to call the method
1969 * DiscreteTime::get_next_step_size(). After many steps, when the simulation
1970 * reaches the end time, the last dt is chosen by the DiscreteTime class in
1971 * such a way that the last step finishes exactly at the end time.
1972 *
1973 * @code
1974 *   template <int dim>
1975 *   void TwoPhaseFlowProblem<dim>::run()
1976 *   {
1977 *   make_grid_and_dofs();
1978 *  
1979 *   {
1980 *   AffineConstraints<double> constraints;
1981 *   constraints.close();
1982 *  
1983 *   VectorTools::project(dof_handler,
1984 *   constraints,
1985 *   QGauss<dim>(degree + 2),
1986 *   InitialValues<dim>(),
1987 *   old_solution);
1988 *   }
1989 *  
1990 *   do
1991 *   {
1992 *   std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
1993 *  
1994 *   assemble_system();
1995 *  
1996 *   solve();
1997 *  
1998 *   output_results();
1999 *  
2000 *   time.advance_time();
2001 *   std::cout << " Now at t=" << time.get_current_time()
2002 *   << ", dt=" << time.get_previous_step_size() << '.'
2003 *   << std::endl
2004 *   << std::endl;
2005 *   }
2006 *   while (time.is_at_end() == false);
2007 *   }
2008 *   } // namespace Step21
2009 *  
2010 *  
2011 * @endcode
2012 *
2013 *
2014 * <a name="step_21-Thecodemaincodefunction"></a>
2015 * <h3>The <code>main</code> function</h3>
2016 *
2017
2018 *
2019 * That's it. In the main function, we pass the degree of the finite element
2020 * space to the constructor of the TwoPhaseFlowProblem object. Here, we use
2021 * zero-th degree elements, i.e. @f$RT_0\times DQ_0 \times DQ_0@f$. The rest is as
2022 * in all the other programs.
2023 *
2024 * @code
2025 *   int main()
2026 *   {
2027 *   try
2028 *   {
2029 *   using namespace Step21;
2030 *  
2031 *   TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2032 *   two_phase_flow_problem.run();
2033 *   }
2034 *   catch (std::exception &exc)
2035 *   {
2036 *   std::cerr << std::endl
2037 *   << std::endl
2038 *   << "----------------------------------------------------"
2039 *   << std::endl;
2040 *   std::cerr << "Exception on processing: " << std::endl
2041 *   << exc.what() << std::endl
2042 *   << "Aborting!" << std::endl
2043 *   << "----------------------------------------------------"
2044 *   << std::endl;
2045 *  
2046 *   return 1;
2047 *   }
2048 *   catch (...)
2049 *   {
2050 *   std::cerr << std::endl
2051 *   << std::endl
2052 *   << "----------------------------------------------------"
2053 *   << std::endl;
2054 *   std::cerr << "Unknown exception!" << std::endl
2055 *   << "Aborting!" << std::endl
2056 *   << "----------------------------------------------------"
2057 *   << std::endl;
2058 *   return 1;
2059 *   }
2060 *  
2061 *   return 0;
2062 *   }
2063 * @endcode
2064<a name="step_21-Results"></a><h1>Results</h1>
2065
2066
2067The code as presented here does not actually compute the results
2068found on the web page. The reason is, that even on a decent
2069computer it runs more than a day. If you want to reproduce these
2070results, modify the end time of the DiscreteTime object to `250` within the
2071constructor of TwoPhaseFlowProblem.
2072
2073If we run the program, we get the following kind of output:
2074@code
2075Number of active cells: 1024
2076Number of degrees of freedom: 4160 (2112+1024+1024)
2077
2078Timestep 1
2079 22 CG Schur complement iterations for pressure.
2080 1 CG iterations for saturation.
2081 Now at t=0.0326742, dt=0.0326742.
2082
2083Timestep 2
2084 17 CG Schur complement iterations for pressure.
2085 1 CG iterations for saturation.
2086 Now at t=0.0653816, dt=0.0327074.
2087
2088Timestep 3
2089 17 CG Schur complement iterations for pressure.
2090 1 CG iterations for saturation.
2091 Now at t=0.0980651, dt=0.0326836.
2092
2093...
2094@endcode
2095As we can see, the time step is pretty much constant right from the start,
2096which indicates that the velocities in the domain are not strongly dependent
2097on changes in saturation, although they certainly are through the factor
2098@f$\lambda(S)@f$ in the pressure equation.
2099
2100Our second observation is that the number of CG iterations needed to solve the
2101pressure Schur complement equation drops from 22 to 17 between the first and
2102the second time step (in fact, it remains around 17 for the rest of the
2103computations). The reason is actually simple: Before we solve for the pressure
2104during a time step, we don't reset the <code>solution</code> variable to
2105zero. The pressure (and the other variables) therefore have the previous time
2106step's values at the time we get into the CG solver. Since the velocities and
2107pressures don't change very much as computations progress, the previous time
2108step's pressure is actually a good initial guess for this time step's
2109pressure. Consequently, the number of iterations we need once we have computed
2110the pressure once is significantly reduced.
2111
2112The final observation concerns the number of iterations needed to solve for
2113the saturation, i.e. one. This shouldn't surprise us too much: the matrix we
2114have to solve with is the mass matrix. However, this is the mass matrix for
2115the @f$DGQ_0@f$ element of piecewise constants where no element couples with the
2116degrees of freedom on neighboring cells. The matrix is therefore a diagonal
2117one, and it is clear that we should be able to invert this matrix in a single
2118CG iteration.
2119
2120
2121With all this, here are a few movies that show how the saturation progresses
2122over time. First, this is for the single crack model, as implemented in the
2123<code>SingleCurvingCrack::KInverse</code> class:
2124
2125<img src="https://www.dealii.org/images/steps/developer/step-21.centerline.gif" alt="">
2126
2127As can be seen, the water rich fluid snakes its way mostly along the
2128high-permeability zone in the middle of the domain, whereas the rest of the
2129domain is mostly impermeable. This and the next movie are generated using
2130<code>n_refinement_steps=7</code>, leading to a @f$128\times 128@f$ mesh with some
213116,000 cells and about 66,000 unknowns in total.
2132
2133
2134The second movie shows the saturation for the random medium model of class
2135<code>RandomMedium::KInverse</code>, where we have randomly distributed
2136centers of high permeability and fluid hops from one of these zones to
2137the next:
2138
2139<img src="https://www.dealii.org/images/steps/developer/step-21.random2d.gif" alt="">
2140
2141
2142Finally, here is the same situation in three space dimensions, on a mesh with
2143<code>n_refinement_steps=5</code>, which produces a mesh of some 32,000 cells
2144and 167,000 degrees of freedom:
2145
2146<img src="https://www.dealii.org/images/steps/developer/step-21.random3d.gif" alt="">
2147
2148To repeat these computations, all you have to do is to change the line
2149@code
2150 TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2151@endcode
2152in the main function to
2153@code
2154 TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
2155@endcode
2156The visualization uses a cloud technique, where the saturation is indicated by
2157colored but transparent clouds for each cell. This way, one can also see
2158somewhat what happens deep inside the domain. A different way of visualizing
2159would have been to show isosurfaces of the saturation evolving over
2160time. There are techniques to plot isosurfaces transparently, so that one can
2161see several of them at the same time like the layers of an onion.
2162
2163So why don't we show such isosurfaces? The problem lies in the way isosurfaces
2164are computed: they require that the field to be visualized is continuous, so
2165that the isosurfaces can be generated by following contours at least across a
2166single cell. However, our saturation field is piecewise constant and
2167discontinuous. If we wanted to plot an isosurface for a saturation @f$S=0.5@f$,
2168chances would be that there is no single point in the domain where that
2169saturation is actually attained. If we had to define isosurfaces in that
2170context at all, we would have to take the interfaces between cells, where one
2171of the two adjacent cells has a saturation greater than and the other cell a
2172saturation less than 0.5. However, it appears that most visualization programs
2173are not equipped to do this kind of transformation.
2174
2175
2176<a name="step-21-extensions"></a>
2177<a name="step_21-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2178
2179
2180There are a number of areas where this program can be improved. Three of them
2181are listed below. All of them are, in fact, addressed in a tutorial program
2182that forms the continuation of the current one: @ref step_43 "step-43".
2183
2184
2185<a name="step_21-Solvers"></a><h4>Solvers</h4>
2186
2187
2188At present, the program is not particularly fast: the 2d random medium
2189computation took about a day for the 1,000 or so time steps. The corresponding
21903d computation took almost two days for 800 time steps. The reason why it
2191isn't faster than this is twofold. First, we rebuild the entire matrix in
2192every time step, although some parts such as the @f$B@f$, @f$B^T@f$, and @f$M^S@f$ blocks
2193never change.
2194
2195Second, we could do a lot better with the solver and
2196preconditioners. Presently, we solve the Schur complement @f$B^TM^u(S)^{-1}B@f$
2197with a CG method, using @f$[B^T (\textrm{diag}(M^u(S)))^{-1} B]^{-1}@f$ as a
2198preconditioner. Applying this preconditioner is expensive, since it involves
2199solving a linear system each time. This may have been appropriate for @ref
2200step_20 "step-20", where we have to solve the entire problem only
2201once. However, here we have to solve it hundreds of times, and in such cases
2202it is worth considering a preconditioner that is more expensive to set up the
2203first time, but cheaper to apply later on.
2204
2205One possibility would be to realize that the matrix we use as preconditioner,
2206@f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ is still sparse, and symmetric on top of
2207that. If one looks at the flow field evolve over time, we also see that while
2208@f$S@f$ changes significantly over time, the pressure hardly does and consequently
2209@f$B^T (\textrm{diag}(M^u(S)))^{-1} B \approx B^T (\textrm{diag}(M^u(S^0)))^{-1}
2210B@f$. In other words, the matrix for the first time step should be a good
2211preconditioner also for all later time steps. With a bit of
2212back-and-forthing, it isn't hard to actually get a representation of it as a
2213SparseMatrix object. We could then hand it off to the SparseMIC class to form
2214a sparse incomplete Cholesky decomposition. To form this decomposition is
2215expensive, but we have to do it only once in the first time step, and can then
2216use it as a cheap preconditioner in the future. We could do better even by
2217using the SparseDirectUMFPACK class that produces not only an incomplete, but
2218a complete decomposition of the matrix, which should yield an even better
2219preconditioner.
2220
2221Finally, why use the approximation @f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ to
2222precondition @f$B^T M^u(S)^{-1} B@f$? The latter matrix, after all, is the mixed
2223form of the Laplace operator on the pressure space, for which we use linear
2224elements. We could therefore build a separate matrix @f$A^p@f$ on the side that
2225directly corresponds to the non-mixed formulation of the Laplacian, for
2226example using the bilinear form @f$(\mathbf{K}\lambda(S^n) \nabla
2227\varphi_i,\nabla\varphi_j)@f$. We could then form an incomplete or complete
2228decomposition of this non-mixed matrix and use it as a preconditioner of the
2229mixed form.
2230
2231Using such techniques, it can reasonably be expected that the solution process
2232will be faster by at least an order of magnitude.
2233
2234
2235<a name="step_21-Timestepping"></a><h4>Time stepping</h4>
2236
2237
2238In the introduction we have identified the time step restriction
2239@f[
2240 \triangle t_{n+1} \le \frac h{|\mathbf{u}^{n+1}(\mathbf{x})|}
2241@f]
2242that has to hold globally, i.e. for all @f$\mathbf x@f$. After discretization, we
2243satisfy it by choosing
2244@f[
2245 \triangle t_{n+1} = \frac {\min_K h_K}{\max_{\mathbf{x}}|\mathbf{u}^{n+1}(\mathbf{x})|}.
2246@f]
2247
2248This restriction on the time step is somewhat annoying: the finer we make the
2249mesh the smaller the time step; in other words, we get punished twice: each
2250time step is more expensive to solve and we have to do more time steps.
2251
2252This is particularly annoying since the majority of the additional work is
2253spent solving the implicit part of the equations, i.e. the pressure-velocity
2254system, whereas it is the hyperbolic transport equation for the saturation
2255that imposes the time step restriction.
2256
2257To avoid this bottleneck, people have invented a number of approaches. For
2258example, they may only re-compute the pressure-velocity field every few time
2259steps (or, if you want, use different time step sizes for the
2260pressure/velocity and saturation equations). This keeps the time step
2261restriction on the cheap explicit part while it makes the solution of the
2262implicit part less frequent. Experiments in this direction are
2263certainly worthwhile; one starting point for such an approach is the paper by
2264Zhangxin Chen, Guanren Huan and Baoyan Li: <i>An improved IMPES method for
2265two-phase flow in porous media</i>, Transport in Porous Media, 54 (2004),
2266pp. 361&mdash;376. There are certainly many other papers on this topic as well, but
2267this one happened to land on our desk a while back.
2268
2269
2270
2271<a name="step_21-Adaptivity"></a><h4>Adaptivity</h4>
2272
2273
2274Adaptivity would also clearly help. Looking at the movies, one clearly sees
2275that most of the action is confined to a relatively small part of the domain
2276(this particularly obvious for the saturation, but also holds for the
2277velocities and pressures). Adaptivity can therefore be expected to keep the
2278necessary number of degrees of freedom low, or alternatively increase the
2279accuracy.
2280
2281On the other hand, adaptivity for time dependent problems is not a trivial
2282thing: we would have to change the mesh every few time steps, and we would
2283have to transport our present solution to the next mesh every time we change
2284it (something that the SolutionTransfer class can help with). These are not
2285insurmountable obstacles, but they do require some additional coding and more
2286than we felt comfortable was worth packing into this tutorial program.
2287 *
2288 *
2289<a name="step_21-PlainProg"></a>
2290<h1> The plain program</h1>
2291@include "step-21.cc"
2292*/
double get_previous_step_size() const
void advance_time()
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &return_value) const override
Definition point.h:111
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
TensorFunction(const time_type initial_time=time_type(0.0))
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
#define AssertDimension(dim1, dim2)
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:564
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)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Event initial
Definition event.cc:64
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void random(DoFHandler< dim, spacedim > &dof_handler)
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)
int(& functions)(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation