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