811 *
const unsigned int = 0)
const override
822 * <a name=
"step_21-Pressureboundaryvalues"></a>
851 * <a name=
"step_21-Saturationboundaryvalues"></a>
886 * <a name=
"step_21-Initialdata"></a>
920 *
virtual void vector_value(
const Point<dim> &p,
932 * <a name=
"step_21-Theinversepermeabilitytensor"></a>
944 * <a name=
"step_21-Singlecurvingcrackpermeability"></a>
968 *
value_list(
const std::vector<
Point<dim>> &points,
973 *
for (
unsigned int p = 0; p < points.size(); ++p)
978 *
std::fabs(points[p][1] - 0.5 - 0.1 *
std::sin(10 * points[p][0]));
985 *
for (
unsigned int d = 0;
d < dim; ++
d)
996 * <a name=
"step_21-Randommediumpermeability"></a>
1007 * not what we expect from classes of this type: they should reliably
1008 * represent the same function.
1012 * The solution to this problem is to make the list of centers a static
1013 * member variable of this class, i.e. there exists exactly one such
1014 * variable for the entire program, rather than for each object of this
1021 *
program,
we can't use a regular member function for that since there may
1022 * not be an object of this type around at the time. The C++ standard
1023 * therefore says that only non-member and static member functions can be
1024 * used to initialize a static variable. We use the latter possibility by
1025 * defining a function <code>get_centers</code> that computes the list of
1026 * center points when called.
1030 * Note that this class works just fine in both 2d and 3d, with the only
1031 * difference being that we use more points in 3d: by experimenting we find
1032 * that we need more exponentials in 3d than in 2d (we have more ground to
1033 * cover, after all, if we want to keep the distance between centers roughly
1034 * equal), so we choose 40 in 2d and 100 in 3d. For any other dimension, the
1035 * function does presently not know what to do so simply throws an exception
1036 * indicating exactly this.
1039 * namespace RandomMedium
1041 * template <int dim>
1042 * class KInverse : public TensorFunction<2, dim>
1046 * : TensorFunction<2, dim>()
1050 * value_list(const std::vector<Point<dim>> &points,
1051 * std::vector<Tensor<2, dim>> &values) const override
1053 * AssertDimension(points.size(), values.size());
1055 * for (unsigned int p = 0; p < points.size(); ++p)
1057 * values[p].clear();
1059 * double permeability = 0;
1060 * for (unsigned int i = 0; i < centers.size(); ++i)
1061 * permeability += std::exp(-(points[p] - centers[i]).norm_square() /
1064 * const double normalized_permeability =
1065 * std::min(std::max(permeability, 0.01), 4.);
1067 * for (unsigned int d = 0; d < dim; ++d)
1068 * values[p][d][d] = 1. / normalized_permeability;
1073 * static std::vector<Point<dim>> centers;
1075 * static std::vector<Point<dim>> get_centers()
1077 * const unsigned int N =
1078 * (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
1080 * std::vector<Point<dim>> centers_list(N);
1081 * for (unsigned int i = 0; i < N; ++i)
1082 * for (unsigned int d = 0; d < dim; ++d)
1083 * centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
1085 * return centers_list;
1091 * template <int dim>
1092 * std::vector<Point<dim>> KInverse<dim>::centers =
1093 * KInverse<dim>::get_centers();
1094 * } // namespace RandomMedium
1101 * <a name="step_21-Theinversemobilityandsaturationfunctions"></a>
1102 * <h3>The inverse mobility and saturation functions</h3>
1106 * There are two more pieces of data that we need to describe, namely the
1107 * inverse mobility function and the saturation curve. Their form is also
1108 * given in the introduction:
1111 * double mobility_inverse(const double S, const double viscosity)
1113 * return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
1116 * double fractional_flow(const double S, const double viscosity)
1118 * return S * S / (S * S + viscosity * (1 - S) * (1 - S));
1126 * <a name="step_21-Linearsolversandpreconditioners"></a>
1127 * <h3>Linear solvers and preconditioners</h3>
1131 * The linear solvers we use are also completely analogous to the ones used
1132 * in @ref step_20 "step-20". The following classes are therefore copied verbatim from
1133 * there. Note that the classes here are not only copied from
1134 * @ref step_20 "step-20", but also duplicate classes in deal.II. In a future version of this
1135 * example, they should be replaced by an efficient method, though. There is a
1136 * single change: if the size of a linear system is small, i.e. when the mesh
1137 * is very coarse, then it is sometimes not sufficient to set a maximum of
1138 * <code>src.size()</code> CG iterations before the solver in the
1139 * <code>vmult()</code> function converges. (This is, of course, a result of
1140 * numerical round-off, since we know that on paper, the CG method converges
1141 * in at most <code>src.size()</code> steps.) As a consequence, we set the
1142 * maximum number of iterations equal to the maximum of the size of the linear
1146 * template <class MatrixType>
1147 * class InverseMatrix : public EnableObserverPointer
1150 * InverseMatrix(const MatrixType &m)
1154 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1156 * SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
1157 * 1e-8 * src.l2_norm());
1158 * SolverCG<Vector<double>> cg(solver_control);
1162 * cg.solve(*matrix, dst, src, PreconditionIdentity());
1166 * const ObserverPointer<const MatrixType> matrix;
1171 * class SchurComplement : public EnableObserverPointer
1174 * SchurComplement(const BlockSparseMatrix<double> &A,
1175 * const InverseMatrix<SparseMatrix<double>> &Minv)
1176 * : system_matrix(&A)
1177 * , m_inverse(&Minv)
1178 * , tmp1(A.block(0, 0).m())
1179 * , tmp2(A.block(0, 0).m())
1182 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1184 * system_matrix->block(0, 1).vmult(tmp1, src);
1185 * m_inverse->vmult(tmp2, tmp1);
1186 * system_matrix->block(1, 0).vmult(dst, tmp2);
1190 * const ObserverPointer<const BlockSparseMatrix<double>> system_matrix;
1191 * const ObserverPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
1193 * mutable Vector<double> tmp1, tmp2;
1198 * class ApproximateSchurComplement : public EnableObserverPointer
1201 * ApproximateSchurComplement(const BlockSparseMatrix<double> &A)
1202 * : system_matrix(&A)
1203 * , tmp1(A.block(0, 0).m())
1204 * , tmp2(A.block(0, 0).m())
1207 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1209 * system_matrix->block(0, 1).vmult(tmp1, src);
1210 * system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
1211 * system_matrix->block(1, 0).vmult(dst, tmp2);
1215 * const ObserverPointer<const BlockSparseMatrix<double>> system_matrix;
1217 * mutable Vector<double> tmp1, tmp2;
1225 * <a name="step_21-codeTwoPhaseFlowProblemcodeclassimplementation"></a>
1226 * <h3><code>TwoPhaseFlowProblem</code> class implementation</h3>
1230 * Here now the implementation of the main class. Much of it is actually
1238 * <a name=
"step_21-TwoPhaseFlowProblemTwoPhaseFlowProblem"></a>
1253 *
template <
int dim>
1260 *
, n_refinement_steps(5)
1270 * <a name=
"step_21-TwoPhaseFlowProblemmake_grid_and_dofs"></a>
1281 *
template <
int dim>
1287 *
dof_handler.distribute_dofs(fe);
1290 *
const std::vector<types::global_dof_index> dofs_per_component =
1292 *
const unsigned int n_u = dofs_per_component[0],
1293 *
n_p = dofs_per_component[dim],
1294 *
n_s = dofs_per_component[dim + 1];
1296 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
1298 *
<<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1299 *
<<
" (" <<
n_u <<
'+' <<
n_p <<
'+' <<
n_s <<
')' << std::endl
1306 *
sparsity_pattern.copy_from(
dsp);
1307 *
system_matrix.reinit(sparsity_pattern);
1318 * <a name=
"step_21-TwoPhaseFlowProblemassemble_system"></a>
1344 *
system_matrix = 0;
1360 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1368 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1385 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1387 *
fe_values.
reinit(cell);
1393 *
Here's the first significant difference: We have to get the values
1394 * of the saturation function of the previous time step at the
1395 * quadrature points. To this end, we can use the
1396 * FEValues::get_function_values (previously already used in @ref step_9 "step-9",
1397 * @ref step_14 "step-14" and @ref step_15 "step-15"), a function that takes a solution vector and
1398 * returns a list of function values at the quadrature points of the
1399 * present cell. In fact, it returns the complete vector-valued
1400 * solution at each quadrature point, i.e. not only the saturation but
1401 * also the velocities and pressure:
1404 * fe_values.get_function_values(old_solution, old_solution_values);
1408 * Then we also have to get the values of the pressure right hand side
1409 * and of the inverse permeability tensor at the quadrature points:
1412 * pressure_right_hand_side.value_list(fe_values.get_quadrature_points(),
1413 * pressure_rhs_values);
1414 * k_inverse.value_list(fe_values.get_quadrature_points(),
1415 * k_inverse_values);
1419 * With all this, we can now loop over all the quadrature points and
1420 * shape functions on this cell and assemble those parts of the matrix
1421 * and right hand side that we deal with in this function. The
1422 * individual terms in the contributions should be self-explanatory
1423 * given the explicit form of the bilinear form stated in the
1427 * for (unsigned int q = 0; q < n_q_points; ++q)
1428 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1430 * const double old_s = old_solution_values[q](dim + 1);
1432 * const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1433 * const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1434 * const double phi_i_p = fe_values[pressure].value(i, q);
1435 * const double phi_i_s = fe_values[saturation].value(i, q);
1437 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1439 * const Tensor<1, dim> phi_j_u =
1440 * fe_values[velocities].value(j, q);
1441 * const double div_phi_j_u =
1442 * fe_values[velocities].divergence(j, q);
1443 * const double phi_j_p = fe_values[pressure].value(j, q);
1444 * const double phi_j_s = fe_values[saturation].value(j, q);
1446 * local_matrix(i, j) +=
1447 * (phi_i_u * k_inverse_values[q] *
1448 * mobility_inverse(old_s, viscosity) * phi_j_u -
1449 * div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u +
1450 * phi_i_s * phi_j_s) *
1455 * (-phi_i_p * pressure_rhs_values[q]) * fe_values.JxW(q);
1461 * Next, we also have to deal with the pressure boundary values. This,
1462 * again is as in @ref step_20 "step-20":
1465 * for (const auto &face : cell->face_iterators())
1466 * if (face->at_boundary())
1468 * fe_face_values.reinit(cell, face);
1470 * pressure_boundary_values.value_list(
1471 * fe_face_values.get_quadrature_points(), boundary_values);
1473 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1474 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1476 * const Tensor<1, dim> phi_i_u =
1477 * fe_face_values[velocities].value(i, q);
1480 * -(phi_i_u * fe_face_values.normal_vector(q) *
1481 * boundary_values[q] * fe_face_values.JxW(q));
1487 * The final step in the loop over all cells is to transfer local
1488 * contributions into the global matrix and right hand side vector:
1491 * cell->get_dof_indices(local_dof_indices);
1492 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1493 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1494 * system_matrix.add(local_dof_indices[i],
1495 * local_dof_indices[j],
1496 * local_matrix(i, j));
1498 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1499 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1506 * So much for assembly of matrix and right hand side. Note that we do not
1507 * have to interpolate and apply boundary values since they have all been
1508 * taken care of in the weak form already.
1516 * <a name="step_21-TwoPhaseFlowProblemassemble_rhs_S"></a>
1517 * <h4>TwoPhaseFlowProblem::assemble_rhs_S</h4>
1521 * As explained in the introduction, we can only evaluate the right hand
1522 * side of the saturation equation once the velocity has been computed. We
1523 * therefore have this separate function to this end.
1526 * template <int dim>
1527 * void TwoPhaseFlowProblem<dim>::assemble_rhs_S()
1529 * const QGauss<dim> quadrature_formula(degree + 2);
1530 * const QGauss<dim - 1> face_quadrature_formula(degree + 2);
1531 * FEValues<dim> fe_values(fe,
1532 * quadrature_formula,
1533 * update_values | update_gradients |
1534 * update_quadrature_points | update_JxW_values);
1535 * FEFaceValues<dim> fe_face_values(fe,
1536 * face_quadrature_formula,
1537 * update_values | update_normal_vectors |
1538 * update_quadrature_points |
1539 * update_JxW_values);
1540 * FEFaceValues<dim> fe_face_values_neighbor(fe,
1541 * face_quadrature_formula,
1544 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1545 * const unsigned int n_q_points = quadrature_formula.size();
1546 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1548 * Vector<double> local_rhs(dofs_per_cell);
1550 * std::vector<Vector<double>> old_solution_values(n_q_points,
1551 * Vector<double>(dim + 2));
1552 * std::vector<Vector<double>> old_solution_values_face(n_face_q_points,
1553 * Vector<double>(dim +
1555 * std::vector<Vector<double>> old_solution_values_face_neighbor(
1556 * n_face_q_points, Vector<double>(dim + 2));
1557 * std::vector<Vector<double>> present_solution_values(n_q_points,
1558 * Vector<double>(dim +
1560 * std::vector<Vector<double>> present_solution_values_face(
1561 * n_face_q_points, Vector<double>(dim + 2));
1563 * std::vector<double> neighbor_saturation(n_face_q_points);
1564 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1566 * SaturationBoundaryValues<dim> saturation_boundary_values;
1568 * const FEValuesExtractors::Scalar saturation(dim + 1);
1570 * for (const auto &cell : dof_handler.active_cell_iterators())
1573 * fe_values.reinit(cell);
1575 * fe_values.get_function_values(old_solution, old_solution_values);
1576 * fe_values.get_function_values(solution, present_solution_values);
1580 * First for the cell terms. These are, following the formulas in the
1581 * introduction, @f$(S^n,\sigma)-(F(S^n) \mathbf{v}^{n+1},\nabla
1582 * \sigma)@f$, where @f$\sigma@f$ is the saturation component of the test
1586 * for (unsigned int q = 0; q < n_q_points; ++q)
1587 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1589 * const double old_s = old_solution_values[q](dim + 1);
1590 * Tensor<1, dim> present_u;
1591 * for (unsigned int d = 0; d < dim; ++d)
1592 * present_u[d] = present_solution_values[q](d);
1594 * const double phi_i_s = fe_values[saturation].value(i, q);
1595 * const Tensor<1, dim> grad_phi_i_s =
1596 * fe_values[saturation].gradient(i, q);
1599 * (time.get_next_step_size() * fractional_flow(old_s, viscosity) *
1600 * present_u * grad_phi_i_s +
1601 * old_s * phi_i_s) *
1607 * Secondly, we have to deal with the flux parts on the face
1608 * boundaries. This was a bit more involved because we first have to
1609 * determine which are the influx and outflux parts of the cell
1610 * boundary. If we have an influx boundary, we need to evaluate the
1611 * saturation on the other side of the face (or the boundary values,
1612 * if we are at the boundary of the domain).
1616 * All this is a bit tricky, but has been explained in some detail
1617 * already in @ref step_9 "step-9". Take a look there how this is supposed to work!
1620 * for (const auto face_no : cell->face_indices())
1622 * fe_face_values.reinit(cell, face_no);
1624 * fe_face_values.get_function_values(old_solution,
1625 * old_solution_values_face);
1626 * fe_face_values.get_function_values(solution,
1627 * present_solution_values_face);
1629 * if (cell->at_boundary(face_no))
1630 * saturation_boundary_values.value_list(
1631 * fe_face_values.get_quadrature_points(), neighbor_saturation);
1634 * const auto neighbor = cell->neighbor(face_no);
1635 * const unsigned int neighbor_face =
1636 * cell->neighbor_of_neighbor(face_no);
1638 * fe_face_values_neighbor.reinit(neighbor, neighbor_face);
1640 * fe_face_values_neighbor.get_function_values(
1641 * old_solution, old_solution_values_face_neighbor);
1643 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1644 * neighbor_saturation[q] =
1645 * old_solution_values_face_neighbor[q](dim + 1);
1649 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1651 * Tensor<1, dim> present_u_face;
1652 * for (unsigned int d = 0; d < dim; ++d)
1653 * present_u_face[d] = present_solution_values_face[q](d);
1655 * const double normal_flux =
1656 * present_u_face * fe_face_values.normal_vector(q);
1658 * const bool is_outflow_q_point = (normal_flux >= 0);
1660 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1662 * time.get_next_step_size() * normal_flux *
1663 * fractional_flow((is_outflow_q_point == true ?
1664 * old_solution_values_face[q](dim + 1) :
1665 * neighbor_saturation[q]),
1667 * fe_face_values[saturation].value(i, q) *
1668 * fe_face_values.JxW(q);
1672 * cell->get_dof_indices(local_dof_indices);
1673 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1674 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1683 * <a name="step_21-TwoPhaseFlowProblemsolve"></a>
1684 * <h4>TwoPhaseFlowProblem::solve</h4>
1688 * After all these preparations, we finally solve the linear system for
1689 * velocity and pressure in the same way as in @ref step_20 "step-20". After that, we have
1690 * to deal with the saturation equation (see below):
1693 * template <int dim>
1694 * void TwoPhaseFlowProblem<dim>::solve()
1696 * const InverseMatrix<SparseMatrix<double>> m_inverse(
1697 * system_matrix.block(0, 0));
1698 * Vector<double> tmp(solution.block(0).size());
1699 * Vector<double> schur_rhs(solution.block(1).size());
1700 * Vector<double> tmp2(solution.block(2).size());
1705 * First the pressure, using the pressure Schur complement of the first
1710 * m_inverse.vmult(tmp, system_rhs.block(0));
1711 * system_matrix.block(1, 0).vmult(schur_rhs, tmp);
1712 * schur_rhs -= system_rhs.block(1);
1715 * SchurComplement schur_complement(system_matrix, m_inverse);
1717 * ApproximateSchurComplement approximate_schur_complement(system_matrix);
1719 * InverseMatrix<ApproximateSchurComplement> preconditioner(
1720 * approximate_schur_complement);
1723 * SolverControl solver_control(solution.block(1).size(),
1724 * 1e-12 * schur_rhs.l2_norm());
1725 * SolverCG<Vector<double>> cg(solver_control);
1727 * cg.solve(schur_complement, solution.block(1), schur_rhs, preconditioner);
1729 * std::cout << " " << solver_control.last_step()
1730 * << " CG Schur complement iterations for pressure." << std::endl;
1739 * system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1741 * tmp += system_rhs.block(0);
1743 * m_inverse.vmult(solution.block(0), tmp);
1748 * Finally, we have to take care of the saturation equation. The first
1749 * business we have here is to determine the time step using the formula
1750 * in the introduction. Knowing the shape of our domain and that we
1751 * created the mesh by regular subdivision of cells, we can compute the
1752 * diameter of each of our cells quite easily (in fact we use the linear
1753 * extensions in coordinate directions of the cells, not the
1754 * diameter). Note that we will learn a more general way to do this in
1755 * @ref step_24 "step-24", where we use the GridTools::minimal_cell_diameter function.
1759 * The maximal velocity we compute using a helper function to compute the
1760 * maximal velocity defined below, and with all this we can evaluate our
1761 * new time step length. We use the method
1762 * DiscreteTime::set_desired_next_time_step() to suggest the new
1763 * calculated value of the time step to the DiscreteTime object. In most
1764 * cases, the time object uses the exact provided value to increment time.
1765 * It some case, the step size may be modified further by the time object.
1766 * For example, if the calculated time increment overshoots the end time,
1767 * it is truncated accordingly.
1770 * time.set_desired_next_step_size(std::pow(0.5, double(n_refinement_steps)) /
1771 * get_maximal_velocity());
1775 * The next step is to assemble the right hand side, and then to pass
1776 * everything on for solution. At the end, we project back saturations
1777 * onto the physically reasonable range:
1782 * SolverControl solver_control(system_matrix.block(2, 2).m(),
1783 * 1e-8 * system_rhs.block(2).l2_norm());
1784 * SolverCG<Vector<double>> cg(solver_control);
1785 * cg.solve(system_matrix.block(2, 2),
1786 * solution.block(2),
1787 * system_rhs.block(2),
1788 * PreconditionIdentity());
1790 * project_back_saturation();
1792 * std::cout << " " << solver_control.last_step()
1793 * << " CG iterations for saturation." << std::endl;
1797 * old_solution = solution;
1804 * <a name="step_21-TwoPhaseFlowProblemoutput_results"></a>
1805 * <h4>TwoPhaseFlowProblem::output_results</h4>
1809 * There is nothing surprising here. Since the program will do a lot of time
1810 * steps, we create an output file only every fifth time step and skip all
1811 * other time steps at the top of the file already.
1815 * When creating file names for output close to the bottom of the function,
1816 * we convert the number of the time step to a string representation that
1817 * is padded by leading zeros to four digits. We do this because this way
1818 * all output file names have the same length, and consequently sort well
1819 * when creating a directory listing.
1822 * template <int dim>
1823 * void TwoPhaseFlowProblem<dim>::output_results() const
1825 * if (time.get_step_number() % 5 != 0)
1828 * std::vector<std::string> solution_names;
1832 * solution_names = {"u", "v", "p", "S"};
1836 * solution_names = {"u", "v", "w", "p", "S"};
1840 * DEAL_II_NOT_IMPLEMENTED();
1843 * DataOut<dim> data_out;
1845 * data_out.attach_dof_handler(dof_handler);
1846 * data_out.add_data_vector(solution, solution_names);
1848 * data_out.build_patches(degree + 1);
1850 * std::ofstream output("solution-" +
1851 * Utilities::int_to_string(time.get_step_number(), 4) +
1853 * data_out.write_vtk(output);
1861 * <a name="step_21-TwoPhaseFlowProblemproject_back_saturation"></a>
1862 * <h4>TwoPhaseFlowProblem::project_back_saturation</h4>
1866 * In this function, we simply run over all saturation degrees of freedom
1867 * and make sure that if they should have left the physically reasonable
1868 * range, that they be reset to the interval @f$[0,1]@f$. To do this, we only
1869 * have to loop over all saturation components of the solution vector; these
1870 * are stored in the block 2 (block 0 are the velocities, block 1 are the
1875 * It may be instructive to note that this function almost never triggers
1876 * when the time step is chosen as mentioned in the introduction. However,
1877 * if we choose the timestep only slightly larger, we get plenty of values
1878 * outside the proper range. Strictly speaking, the function is therefore
1879 * unnecessary if we choose the time step small enough. In a sense, the
1880 * function is therefore only a safety device to avoid situations where our
1881 * entire solution becomes unphysical because individual degrees of freedom
1882 * have become unphysical a few time steps earlier.
1885 * template <int dim>
1886 * void TwoPhaseFlowProblem<dim>::project_back_saturation()
1888 * for (unsigned int i = 0; i < solution.block(2).size(); ++i)
1889 * if (solution.block(2)(i) < 0)
1890 * solution.block(2)(i) = 0;
1891 * else if (solution.block(2)(i) > 1)
1892 * solution.block(2)(i) = 1;
1899 * <a name="step_21-TwoPhaseFlowProblemget_maximal_velocity"></a>
1900 * <h4>TwoPhaseFlowProblem::get_maximal_velocity</h4>
1904 * The following function is used in determining the maximal allowable time
1905 * step. What it does is to loop over all quadrature points in the domain
1906 * and find what the maximal magnitude of the velocity is.
1909 * template <int dim>
1910 * double TwoPhaseFlowProblem<dim>::get_maximal_velocity() const
1912 * const QGauss<dim> quadrature_formula(degree + 2);
1913 * const unsigned int n_q_points = quadrature_formula.size();
1915 * FEValues<dim> fe_values(fe, quadrature_formula, update_values);
1916 * std::vector<Vector<double>> solution_values(n_q_points,
1917 * Vector<double>(dim + 2));
1918 * double max_velocity = 0;
1920 * for (const auto &cell : dof_handler.active_cell_iterators())
1922 * fe_values.reinit(cell);
1923 * fe_values.get_function_values(solution, solution_values);
1925 * for (unsigned int q = 0; q < n_q_points; ++q)
1927 * Tensor<1, dim> velocity;
1928 * for (unsigned int i = 0; i < dim; ++i)
1929 * velocity[i] = solution_values[q](i);
1931 * max_velocity = std::max(max_velocity, velocity.norm());
1935 * return max_velocity;
1942 * <a name="step_21-TwoPhaseFlowProblemrun"></a>
1943 * <h4>TwoPhaseFlowProblem::run</h4>
1947 * This is the final function of our main class. Its brevity speaks for
1948 * itself. There are only two points worth noting: First, the function
1949 * projects the initial values onto the finite element space at the
1950 * beginning; the VectorTools::project function doing this requires an
1951 * argument indicating the hanging node constraints. We have none in this
1952 * program (we compute on a uniformly refined mesh), but the function
1953 * requires the argument anyway, of course. So we have to create a
1954 * constraint object. In its original state, constraint objects are
1955 * unsorted, and have to be sorted (using the AffineConstraints::close
1956 * function) before they can be used. This is what we do here, and which is
1983 *
constraints.close();
1994 *
std::cout <<
"Timestep " << time.get_step_number() + 1 << std::endl;
2002 *
time.advance_time();
2003 *
std::cout <<
" Now at t=" << time.get_current_time()
2004 *
<<
", dt=" << time.get_previous_step_size() <<
'.'
2008 *
while (time.is_at_end() ==
false);
2016 * <a name=
"step_21-Thecodemaincodefunction"></a>
2021 *
That's it. In the main function, we pass the degree of the finite element
2022 * space to the constructor of the TwoPhaseFlowProblem object. Here, we use
2023 * zero-th degree elements, i.e. @f$RT_0\times DQ_0 \times DQ_0@f$. The rest is as
2024 * in all the other programs.
2031 * using namespace Step21;
2033 * TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2034 * two_phase_flow_problem.run();
2036 * catch (std::exception &exc)
2038 * std::cerr << std::endl
2040 * << "----------------------------------------------------"
2042 * std::cerr << "Exception on processing: " << std::endl
2043 * << exc.what() << std::endl
2044 * << "Aborting!" << std::endl
2045 * << "----------------------------------------------------"
2052 * std::cerr << std::endl
2054 * << "----------------------------------------------------"
2056 * std::cerr << "Unknown exception!" << std::endl
2057 * << "Aborting!" << std::endl
2058 * << "----------------------------------------------------"
2066<a name="step_21-Results"></a><h1>Results</h1>
2069The code as presented here does not actually compute the results
2070found on the web page. The reason is, that even on a decent
2071computer it runs more than a day. If you want to reproduce these
2072results, modify the end time of the DiscreteTime object to `250` within the
2073constructor of TwoPhaseFlowProblem.
2075If we run the program, we get the following kind of output:
2077Number of active cells: 1024
2078Number of degrees of freedom: 4160 (2112+1024+1024)
2081 22 CG Schur complement iterations for pressure.
2082 1 CG iterations for saturation.
2083 Now at t=0.0326742, dt=0.0326742.
2086 17 CG Schur complement iterations for pressure.
2087 1 CG iterations for saturation.
2088 Now at t=0.0653816, dt=0.0327074.
2091 17 CG Schur complement iterations for pressure.
2092 1 CG iterations for saturation.
2093 Now at t=0.0980651, dt=0.0326836.
2097As we can see, the time step is pretty much constant right from the start,
2098which indicates that the velocities in the domain are not strongly dependent
2099on changes in saturation, although they certainly are through the factor
2100@f$\lambda(S)@f$ in the pressure equation.
2102Our second observation is that the number of CG iterations needed to solve the
2103pressure Schur complement equation drops from 22 to 17 between the first and
2104the second time step (in fact, it remains around 17 for the rest of the
2105computations). The reason is actually simple: Before we solve for the pressure
2108step
's values at the time we get into the CG solver. Since the velocities and
2110step
's pressure is actually a good initial guess for this time step's
2116have to solve with is the mass matrix. However, this is the mass matrix for
2117the @f$DGQ_0@f$ element of piecewise constants where no element couples with the
2118degrees of freedom on neighboring cells. The matrix is therefore a diagonal
2119one, and it is clear that we should be able to invert this matrix in a single
2123With all this, here are a few movies that show how the saturation progresses
2124over time. First, this is for the single crack model, as implemented in the
2125<code>SingleCurvingCrack::KInverse</code> class:
2127<img src="https://www.dealii.org/images/steps/developer/step-21.centerline.gif" alt="">
2129As can be seen, the water rich fluid snakes its way mostly along the
2130high-permeability zone in the middle of the domain, whereas the rest of the
2131domain is mostly impermeable. This and the next movie are generated using
2132<code>n_refinement_steps=7</code>, leading to a @f$128\times 128@f$ mesh with some
213316,000 cells and about 66,000 unknowns in total.
2136The second movie shows the saturation for the random medium model of class
2137<code>RandomMedium::KInverse</code>, where we have randomly distributed
2138centers of high permeability and fluid hops from one of these zones to
2141<img src="https://www.dealii.org/images/steps/developer/step-21.random2d.gif" alt="">
2144Finally, here is the same situation in three space dimensions, on a mesh with
2145<code>n_refinement_steps=5</code>, which produces a mesh of some 32,000 cells
2146and 167,000 degrees of freedom:
2148<img src="https://www.dealii.org/images/steps/developer/step-21.random3d.gif" alt="">
2150To repeat these computations, all you have to do is to change the line
2152 TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2154in the main function to
2156 TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
2158The visualization uses a cloud technique, where the saturation is indicated by
2159colored but transparent clouds for each cell. This way, one can also see
2160somewhat what happens deep inside the domain. A different way of visualizing
2161would have been to show isosurfaces of the saturation evolving over
2162time. There are techniques to plot isosurfaces transparently, so that one can
2163see several of them at the same time like the layers of an onion.
2178<a name=
"step-21-extensions"></a>
2179<a name=
"step_21-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
2187<a name=
"step_21-Solvers"></a><
h4>Solvers</
h4>
2193isn't faster than this is twofold. First, we rebuild the entire matrix in
2194every time step, although some parts such as the @f$B@f$, @f$B^T@f$, and @f$M^S@f$ blocks
2197Second, we could do a lot better with the solver and
2198preconditioners. Presently, we solve the Schur complement @f$B^TM^u(S)^{-1}B@f$
2199with a CG method, using @f$[B^T (\textrm{diag}(M^u(S)))^{-1} B]^{-1}@f$ as a
2200preconditioner. Applying this preconditioner is expensive, since it involves
2201solving a linear system each time. This may have been appropriate for @ref
2202step_20 "step-20", where we have to solve the entire problem only
2203once. However, here we have to solve it hundreds of times, and in such cases
2204it is worth considering a preconditioner that is more expensive to set up the
2205first time, but cheaper to apply later on.
2207One possibility would be to realize that the matrix we use as preconditioner,
2208@f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ is still sparse, and symmetric on top of
2209that. If one looks at the flow field evolve over time, we also see that while
2210@f$S@f$ changes significantly over time, the pressure hardly does and consequently
2211@f$B^T (\textrm{diag}(M^u(S)))^{-1} B \approx B^T (\textrm{diag}(M^u(S^0)))^{-1}
2212B@f$. In other words, the matrix for the first time step should be a good
2213preconditioner also for all later time steps. With a bit of
2291<a name=
"step_21-PlainProg"></a>
double get_previous_step_size() 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
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#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())
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.
std::vector< index_type > data
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)
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.
constexpr types::blas_int zero
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
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