1083 *
const unsigned int )
const
1090 *
template <
int dim>
1092 * PressureBoundaryValues<dim>::value(
const Point<dim> &p,
1093 *
const unsigned int )
const
1095 *
return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
1096 * alpha * p[0] * p[0] * p[0] / 6);
1101 *
template <
int dim>
1102 *
void ExactSolution<dim>::vector_value(
const Point<dim> &p,
1107 *
values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2;
1108 *
values(1) = alpha * p[0] * p[1];
1109 *
values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
1110 * alpha * p[0] * p[0] * p[0] / 6);
1118 * <a name=
"Theinversepermeabilitytensor"></a>
1119 * <h3>The inverse permeability tensor</h3>
1123 * In addition to the other equation data, we also want to use a
1124 * permeability tensor, or better -- because
this is all that appears in the
1125 * weak form -- the inverse of the permeability tensor,
1126 * <code>KInverse</code>. For the purpose of verifying the exactness of the
1127 * solution and determining convergence orders,
this tensor is more in the
1132 * However, a spatially varying permeability tensor is indispensable in
1133 * real-life porous media flow simulations, and we would like to use the
1134 * opportunity to demonstrate the technique to use tensor valued
functions.
1138 * Possibly unsurprisingly, deal.II also has a base
class not only for
1140 * base
class) but also
for functions that
return tensors of fixed dimension
1141 * and rank, the <code>
TensorFunction</code>
template. Here, the function
1142 * under consideration returns a dim-by-dim
matrix, i.e. a tensor of rank 2
1143 * and dimension <code>dim</code>. We then choose the
template arguments of
1144 * the base
class appropriately.
1148 * The
interface that the <code>
TensorFunction</code>
class provides is
1149 * essentially equivalent to the <code>
Function</code>
class. In particular,
1150 * there exists a <code>value_list</code> function that takes a list of
1151 * points at which to evaluate the function, and returns the
values of the
1152 * function in the
second argument, a list of tensors:
1155 *
template <
int dim>
1171 * The implementation is less interesting. As in previous examples, we add a
1172 *
check to the beginning of the
class to make sure that the sizes of input
1173 * and output parameters are the same (see @ref step_5
"step-5" for a discussion of
this
1174 * technique). Then we
loop over all evaluation points, and
for each one
1179 * There is an oddity at the top of the function (the
1180 * `(
void)points;` statement) that is worth discussing. The values
1181 * we put into the output `values` array does not actually depend
1182 * on the `points` arrays of coordinates at which the function is
1183 * evaluated. In other words, the `points` argument is in fact
1184 * unused, and we could have just not given it a name
if we had
1185 * wanted. But we want to use the `points`
object for checking
1186 * that the `values`
object has the correct size. The problem is
1188 * that expands to nothing; the compiler will then complain that
1189 * the `points`
object is unused. The idiomatic approach to
1190 * silencing
this warning is to have a statement that evaluates
1191 * (reads) variable but doesn
't actually do anything: That's what
1192 * `(
void)points;` does: It reads from `points`, and then casts
1193 * the result of the read to `
void`, i.e.,
nothing. This statement
1194 * is, in other words, completely pointless and implies no actual
1195 * action except to explain to the compiler that yes,
this
1196 * variable is in fact used even in release mode. (In debug mode,
1198 * from the variable, and so the funny statement would not be
1199 * necessary in debug mode.)
1202 *
template <
int dim>
1203 *
void KInverse<dim>::value_list(
const std::vector<
Point<dim>> &points,
1209 *
for (
auto &value :
values)
1219 * <a name=
"MixedLaplaceProblemclassimplementation"></a>
1220 * <h3>MixedLaplaceProblem
class implementation</h3>
1225 * <a name=
"MixedLaplaceProblemMixedLaplaceProblem"></a>
1226 * <h4>MixedLaplaceProblem::MixedLaplaceProblem</h4>
1230 * In the constructor of
this class, we
first store the
value that was
1231 * passed in concerning the degree of the finite elements we shall use (a
1232 * degree of zero,
for example, means to use RT(0) and DG(0)), and then
1233 * construct the vector valued element belonging to the space @f$X_h@f$ described
1234 * in the introduction. The rest of the constructor is as in the early
1235 * tutorial programs.
1239 * The only thing worth describing here is the constructor
call of the
1240 * <code>fe</code> variable. The <code>
FESystem</code>
class to which this
1241 * variable belongs has a number of different constructors that all refer to
1242 * binding simpler elements together into one larger element. In the present
1243 *
case, we want to couple a single RT(degree) element with a single
1244 * DQ(degree) element. The constructor to <code>
FESystem</code> that does
1245 *
this requires us to specify
first the
first base element (the
1247 * of copies
for this base element, and then similarly the kind and number
1248 * of <code>
FE_DGQ</code> elements. Note that the Raviart-Thomas element
1249 * already has <code>dim</code> vector components, so that the coupled
1250 * element will have <code>dim+1</code> vector components, the
first
1251 * <code>dim</code> of which correspond to the velocity variable whereas the
1252 * last one corresponds to the pressure.
1256 * It is also worth comparing the way we constructed
this element from its
1257 * base elements, with the way we have done so in @ref step_8
"step-8": there, we have
1258 * built it as <code>fe (
FE_Q@<dim@>(1), dim)</code>, i.e. we have simply
1259 * used <code>dim</code> copies of the <code>
FE_Q(1)</code> element, one
1260 *
copy for the displacement in each coordinate direction.
1263 *
template <
int dim>
1264 * MixedLaplaceProblem<dim>::MixedLaplaceProblem(
const unsigned int degree)
1275 * <a name=
"MixedLaplaceProblemmake_grid_and_dofs"></a>
1276 * <h4>MixedLaplaceProblem::make_grid_and_dofs</h4>
1280 * This next function starts out with well-known
functions calls that create
1281 * and
refine a mesh, and then associate degrees of freedom with it:
1284 *
template <
int dim>
1285 *
void MixedLaplaceProblem<dim>::make_grid_and_dofs()
1290 * dof_handler.distribute_dofs(fe);
1294 * However, then things become different. As mentioned in the
1295 * introduction, we want to subdivide the
matrix into blocks corresponding
1296 * to the two different kinds of variables, velocity and pressure. To
this
1297 *
end, we
first have to make sure that the indices corresponding to
1298 * velocities and pressures are not intermingled: First all velocity
1299 * degrees of freedom, then all pressure DoFs. This way, the global
matrix
1300 * separates nicely into a @f$2 \times 2@f$ system. To achieve
this, we have to
1301 * renumber degrees of freedom based on their vector component, an
1302 * operation that conveniently is already implemented:
1309 * The next thing is that we want to figure out the sizes of these blocks
1310 * so that we can allocate an appropriate amount of space. To
this end, we
1312 * counts how many shape functions are non-zero for a particular vector
1313 * component. We have <code>dim+1</code> vector components, and
1314 *
DoFTools::count_dofs_per_fe_component() will count how many shape
1315 * functions belong to each of these components.
1319 * There is one problem here. As described in the documentation of that
1320 * function, it <i>wants</i> to put the number of @f$x@f$-velocity shape
1321 * functions into <code>dofs_per_component[0]</code>, the number of
1322 * @f$y@f$-velocity shape functions into <code>dofs_per_component[1]</code>
1323 * (and similar in 3d), and the number of pressure shape functions into
1324 * <code>dofs_per_component[dim]</code>. But, the Raviart-Thomas element
1325 * is special in that it is non-@ref GlossPrimitive "primitive", i.e.,
1326 * for Raviart-Thomas elements all velocity shape functions
1327 * are nonzero in all components. In other words, the function cannot
1328 * distinguish between @f$x@f$ and @f$y@f$ velocity functions because there
1329 * <i>is</i> no such distinction. It therefore puts the overall number
1330 * of velocity into each of <code>dofs_per_component[c]</code>,
1331 * @f$0\le c\le \text{dim}@f$. On the other hand, the number
1332 * of pressure variables equals the number of shape
functions that are
1333 *
nonzero in the dim-th component.
1337 * Using
this knowledge, we can get the number of velocity shape
1338 *
functions from any of the first <code>dim</code> elements of
1339 * <code>dofs_per_component</code>, and then use
this below to initialize
1340 * the vector and
matrix block sizes, as well as create output.
1344 * @note If you find
this concept difficult to understand, you may
1346 * instead, as we do in the corresponding piece of code in @ref step_22 "step-22".
1347 * You might also want to read up on the difference between
1348 * @ref GlossBlock "blocks" and @ref GlossComponent "components"
1354 *
const unsigned int n_u = dofs_per_component[0],
1355 * n_p = dofs_per_component[dim];
1357 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
1361 * <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1362 * <<
" (" << n_u <<
'+' << n_p <<
')' << std::endl;
1366 * The next task is to allocate a sparsity pattern
for the
matrix that we
1367 * will create. We use a compressed sparsity pattern like in the previous
1368 * steps, but as <code>system_matrix</code> is a block
matrix we use the
1371 * four blocks in a @f$2 \times 2@f$ pattern. The blocks
' sizes depend on
1372 * <code>n_u</code> and <code>n_p</code>, which hold the number of velocity
1373 * and pressure variables.
1376 * const std::vector<types::global_dof_index> block_sizes = {n_u, n_p};
1377 * BlockDynamicSparsityPattern dsp(block_sizes, block_sizes);
1378 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
1382 * We use the compressed block sparsity pattern in the same way as the
1383 * non-block version to create the sparsity pattern and then the system
1387 * sparsity_pattern.copy_from(dsp);
1388 * system_matrix.reinit(sparsity_pattern);
1392 * Then we have to resize the solution and right hand side vectors in
1393 * exactly the same way as the block compressed sparsity pattern:
1396 * solution.reinit(block_sizes);
1397 * system_rhs.reinit(block_sizes);
1404 * <a name="MixedLaplaceProblemassemble_system"></a>
1405 * <h4>MixedLaplaceProblem::assemble_system</h4>
1409 * Similarly, the function that assembles the linear system has mostly been
1410 * discussed already in the introduction to this example. At its top, what
1411 * happens are all the usual steps, with the addition that we do not only
1412 * allocate quadrature and <code>FEValues</code> objects for the cell terms,
1413 * but also for face terms. After that, we define the usual abbreviations
1414 * for variables, and the allocate space for the local matrix and right hand
1415 * side contributions, and the array that holds the global numbers of the
1416 * degrees of freedom local to the present cell.
1419 * template <int dim>
1420 * void MixedLaplaceProblem<dim>::assemble_system()
1422 * QGauss<dim> quadrature_formula(degree + 2);
1423 * QGauss<dim - 1> face_quadrature_formula(degree + 2);
1425 * FEValues<dim> fe_values(fe,
1426 * quadrature_formula,
1427 * update_values | update_gradients |
1428 * update_quadrature_points | update_JxW_values);
1429 * FEFaceValues<dim> fe_face_values(fe,
1430 * face_quadrature_formula,
1431 * update_values | update_normal_vectors |
1432 * update_quadrature_points |
1433 * update_JxW_values);
1435 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1436 * const unsigned int n_q_points = quadrature_formula.size();
1437 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1439 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1440 * Vector<double> local_rhs(dofs_per_cell);
1442 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1446 * The next step is to declare objects that represent the source term,
1447 * pressure boundary value, and coefficient in the equation. In addition
1448 * to these objects that represent continuous functions, we also need
1449 * arrays to hold their values at the quadrature points of individual
1450 * cells (or faces, for the boundary values). Note that in the case of the
1451 * coefficient, the array has to be one of matrices.
1454 * const PrescribedSolution::RightHandSide<dim> right_hand_side;
1455 * const PrescribedSolution::PressureBoundaryValues<dim>
1456 * pressure_boundary_values;
1457 * const PrescribedSolution::KInverse<dim> k_inverse;
1459 * std::vector<double> rhs_values(n_q_points);
1460 * std::vector<double> boundary_values(n_face_q_points);
1461 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1465 * Finally, we need a couple of extractors that we will use to get at the
1466 * velocity and pressure components of vector-valued shape
1467 * functions. Their function and use is described in detail in the @ref
1468 * vector_valued report. Essentially, we will use them as subscripts on
1469 * the FEValues objects below: the FEValues object describes all vector
1470 * components of shape functions, while after subscription, it will only
1471 * refer to the velocities (a set of <code>dim</code> components starting
1472 * at component zero) or the pressure (a scalar component located at
1473 * position <code>dim</code>):
1476 * const FEValuesExtractors::Vector velocities(0);
1477 * const FEValuesExtractors::Scalar pressure(dim);
1481 * With all this in place, we can go on with the loop over all cells. The
1482 * body of this loop has been discussed in the introduction, and will not
1483 * be commented any further here:
1486 * for (const auto &cell : dof_handler.active_cell_iterators())
1488 * fe_values.reinit(cell);
1492 * right_hand_side.value_list(fe_values.get_quadrature_points(),
1494 * k_inverse.value_list(fe_values.get_quadrature_points(),
1495 * k_inverse_values);
1497 * for (unsigned int q = 0; q < n_q_points; ++q)
1498 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1500 * const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1501 * const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1502 * const double phi_i_p = fe_values[pressure].value(i, q);
1504 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1506 * const Tensor<1, dim> phi_j_u =
1507 * fe_values[velocities].value(j, q);
1508 * const double div_phi_j_u =
1509 * fe_values[velocities].divergence(j, q);
1510 * const double phi_j_p = fe_values[pressure].value(j, q);
1512 * local_matrix(i, j) +=
1513 * (phi_i_u * k_inverse_values[q] * phi_j_u
1514 * - phi_i_p * div_phi_j_u
1515 * - div_phi_i_u * phi_j_p)
1516 * * fe_values.JxW(q);
1519 * local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
1522 * for (const auto &face : cell->face_iterators())
1523 * if (face->at_boundary())
1525 * fe_face_values.reinit(cell, face);
1527 * pressure_boundary_values.value_list(
1528 * fe_face_values.get_quadrature_points(), boundary_values);
1530 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1531 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1532 * local_rhs(i) += -(fe_face_values[velocities].value(i, q) *
1533 * fe_face_values.normal_vector(q) *
1534 * boundary_values[q] *
1535 * fe_face_values.JxW(q));
1540 * The final step in the loop over all cells is to transfer local
1541 * contributions into the global matrix and right hand side
1542 * vector. Note that we use exactly the same interface as in previous
1543 * examples, although we now use block matrices and vectors instead of
1544 * the regular ones. In other words, to the outside world, block
1545 * objects have the same interface as matrices and vectors, but they
1546 * additionally allow to access individual blocks.
1549 * cell->get_dof_indices(local_dof_indices);
1550 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1551 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1552 * system_matrix.add(local_dof_indices[i],
1553 * local_dof_indices[j],
1554 * local_matrix(i, j));
1555 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1556 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1564 * <a name="Implementationoflinearsolversandpreconditioners"></a>
1565 * <h3>Implementation of linear solvers and preconditioners</h3>
1569 * The linear solvers and preconditioners we use in this example have
1570 * been discussed in significant detail already in the introduction. We
1571 * will therefore not discuss the rationale for our approach here any
1572 * more, but rather only comment on some remaining implementational
1578 * <a name="MixedLaplacesolve"></a>
1579 * <h4>MixedLaplace::solve</h4>
1583 * As already outlined in the introduction, the solve function consists
1584 * essentially of two steps. First, we have to form the first equation
1585 * involving the Schur complement and solve for the pressure (component 1
1586 * of the solution). Then, we can reconstruct the velocities from the
1587 * second equation (component 0 of the solution).
1590 * template <int dim>
1591 * void MixedLaplaceProblem<dim>::solve()
1595 * As a first step we declare references to all block components of the
1596 * matrix, the right hand side and the solution vector that we will
1600 * const auto &M = system_matrix.block(0, 0);
1601 * const auto &B = system_matrix.block(0, 1);
1603 * const auto &F = system_rhs.block(0);
1604 * const auto &G = system_rhs.block(1);
1606 * auto &U = solution.block(0);
1607 * auto &P = solution.block(1);
1611 * Then, we will create corresponding LinearOperator objects and create
1612 * the <code>op_M_inv</code> operator:
1615 * const auto op_M = linear_operator(M);
1616 * const auto op_B = linear_operator(B);
1618 * ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10);
1619 * SolverCG<Vector<double>> solver_M(reduction_control_M);
1620 * PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
1622 * preconditioner_M.initialize(M);
1624 * const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
1628 * This allows us to declare the Schur complement <code>op_S</code> and
1629 * the approximate Schur complement <code>op_aS</code>:
1632 * const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
1633 * const auto op_aS =
1634 * transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
1638 * We now create a preconditioner out of <code>op_aS</code> that
1639 * applies a fixed number of 30 (inexpensive) CG iterations:
1642 * IterationNumberControl iteration_number_control_aS(30, 1.e-18);
1643 * SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
1645 * const auto preconditioner_S =
1646 * inverse_operator(op_aS, solver_aS, PreconditionIdentity());
1650 * Now on to the first equation. The right hand side of it is
1651 * @f$B^TM^{-1}F-G@f$, which is what we compute in the first few lines. We
1652 * then solve the first equation with a CG solver and the
1653 * preconditioner we just declared.
1656 * const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
1658 * SolverControl solver_control_S(2000, 1.e-12);
1659 * SolverCG<Vector<double>> solver_S(solver_control_S);
1661 * const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
1663 * P = op_S_inv * schur_rhs;
1665 * std::cout << solver_control_S.last_step()
1666 * << " CG Schur complement iterations to obtain convergence."
1671 * After we have the pressure, we can compute the velocity. The equation
1672 * reads @f$MU=-BP+F@f$, and we solve it by first computing the right hand
1673 * side, and then multiplying it with the object that represents the
1674 * inverse of the @ref GlossMassMatrix "mass matrix":
1677 * U = op_M_inv * (F - op_B * P);
1684 * <a name="MixedLaplaceProblemclassimplementationcontinued"></a>
1685 * <h3>MixedLaplaceProblem class implementation (continued)</h3>
1690 * <a name="MixedLaplacecompute_errors"></a>
1691 * <h4>MixedLaplace::compute_errors</h4>
1695 * After we have dealt with the linear solver and preconditioners, we
1696 * continue with the implementation of our main class. In particular, the
1697 * next task is to compute the errors in our numerical solution, in both the
1698 * pressures as well as velocities.
1702 * To compute errors in the solution, we have already introduced the
1703 * <code>VectorTools::integrate_difference</code> function in @ref step_7 "step-7" and
1704 * @ref step_11 "step-11". However, there we only dealt with scalar solutions, whereas here
1705 * we have a vector-valued solution with components that even denote
1706 * different quantities and may have different orders of convergence (this
1707 * isn't the
case here, by choice of the used finite elements, but is
1708 * frequently the
case in mixed finite element applications). What we
1709 * therefore have to
do is to `mask
' the components that we are interested
1710 * in. This is easily done: the
1711 * <code>VectorTools::integrate_difference</code> function takes as one of its
1712 * arguments a pointer to a weight function (the parameter defaults to the
1713 * null pointer, meaning unit weights). What we have to do is to pass
1714 * a function object that equals one in the components we are interested in,
1715 * and zero in the other ones. For example, to compute the pressure error,
1716 * we should pass a function that represents the constant vector with a unit
1717 * value in component <code>dim</code>, whereas for the velocity the
1718 * constant vector should be one in the first <code>dim</code> components,
1719 * and zero in the location of the pressure.
1723 * In deal.II, the <code>ComponentSelectFunction</code> does exactly this:
1724 * it wants to know how many vector components the function it is to
1725 * represent should have (in our case this would be <code>dim+1</code>, for
1726 * the joint velocity-pressure space) and which individual or range of
1727 * components should be equal to one. We therefore define two such masks at
1728 * the beginning of the function, following by an object representing the
1729 * exact solution and a vector in which we will store the cellwise errors as
1730 * computed by <code>integrate_difference</code>:
1733 * template <int dim>
1734 * void MixedLaplaceProblem<dim>::compute_errors() const
1736 * const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
1737 * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
1740 * PrescribedSolution::ExactSolution<dim> exact_solution;
1741 * Vector<double> cellwise_errors(triangulation.n_active_cells());
1745 * As already discussed in @ref step_7 "step-7", we have to realize that it is
1746 * impossible to integrate the errors exactly. All we can do is
1747 * approximate this integral using quadrature. This actually presents a
1748 * slight twist here: if we naively chose an object of type
1749 * <code>QGauss@<dim@>(degree+1)</code> as one may be inclined to do (this
1750 * is what we used for integrating the linear system), one realizes that
1751 * the error is very small and does not follow the expected convergence
1752 * curves at all. What is happening is that for the mixed finite elements
1753 * used here, the Gauss points happen to be superconvergence points in
1754 * which the pointwise error is much smaller (and converges with higher
1755 * order) than anywhere else. These are therefore not particularly good
1756 * points for integration. To avoid this problem, we simply use a
1757 * trapezoidal rule and iterate it <code>degree+2</code> times in each
1758 * coordinate direction (again as explained in @ref step_7 "step-7"):
1761 * QTrapezoid<1> q_trapez;
1762 * QIterated<dim> quadrature(q_trapez, degree + 2);
1766 * With this, we can then let the library compute the errors and output
1767 * them to the screen:
1770 * VectorTools::integrate_difference(dof_handler,
1775 * VectorTools::L2_norm,
1777 * const double p_l2_error =
1778 * VectorTools::compute_global_error(triangulation,
1780 * VectorTools::L2_norm);
1782 * VectorTools::integrate_difference(dof_handler,
1787 * VectorTools::L2_norm,
1789 * const double u_l2_error =
1790 * VectorTools::compute_global_error(triangulation,
1792 * VectorTools::L2_norm);
1794 * std::cout << "Errors: ||e_p||_L2 = " << p_l2_error
1795 * << ", ||e_u||_L2 = " << u_l2_error << std::endl;
1802 * <a name="MixedLaplaceoutput_results"></a>
1803 * <h4>MixedLaplace::output_results</h4>
1807 * The last interesting function is the one in which we generate graphical
1808 * output. Note that all velocity components get the same solution name
1809 * "u". Together with using
1810 * DataComponentInterpretation::component_is_part_of_vector this will
1811 * cause DataOut<dim>::write_vtu() to generate a vector representation of
1812 * the individual velocity components, see @ref step_22 "step-22" or the
1813 * @ref VVOutput "Generating graphical output"
1815 * @ref vector_valued
1816 * module for more information. Finally, it seems inappropriate for higher
1817 * order elements to only show a single bilinear quadrilateral per cell in
1818 * the graphical output. We therefore generate patches of size
1819 * (degree+1)x(degree+1) to capture the full information content of the
1820 * solution. See the @ref step_7 "step-7" tutorial program for more information on this.
1823 * template <int dim>
1824 * void MixedLaplaceProblem<dim>::output_results() const
1826 * std::vector<std::string> solution_names(dim, "u");
1827 * solution_names.emplace_back("p");
1828 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1829 * interpretation(dim,
1830 * DataComponentInterpretation::component_is_part_of_vector);
1831 * interpretation.push_back(DataComponentInterpretation::component_is_scalar);
1833 * DataOut<dim> data_out;
1834 * data_out.add_data_vector(dof_handler,
1839 * data_out.build_patches(degree + 1);
1841 * std::ofstream output("solution.vtu");
1842 * data_out.write_vtu(output);
1850 * <a name="MixedLaplacerun"></a>
1851 * <h4>MixedLaplace::run</h4>
1855 * This is the final function of our main class. It's only job is to call
1856 * the other functions in their natural order:
1859 * template <
int dim>
1860 * void MixedLaplaceProblem<dim>::
run()
1862 * make_grid_and_dofs();
1863 * assemble_system();
1874 * <a name=
"Thecodemaincodefunction"></a>
1875 * <h3>The <code>main</code> function</h3>
1879 * The main function we stole from @ref step_6
"step-6" instead of @ref step_4
"step-4". It is almost
1880 *
equal to the one in @ref step_6
"step-6" (apart from the changed
class names, of course),
1881 * the only exception is that we pass the degree of the finite element space
1882 * to the constructor of the mixed Laplace problem (here, we use zero-th order
1890 *
using namespace Step20;
1892 *
const unsigned int fe_degree = 0;
1893 * MixedLaplaceProblem<2> mixed_laplace_problem(fe_degree);
1894 * mixed_laplace_problem.run();
1896 *
catch (std::exception &exc)
1898 * std::cerr << std::endl
1900 * <<
"----------------------------------------------------"
1902 * std::cerr <<
"Exception on processing: " << std::endl
1903 * << exc.what() << std::endl
1904 * <<
"Aborting!" << std::endl
1905 * <<
"----------------------------------------------------"
1912 * std::cerr << std::endl
1914 * <<
"----------------------------------------------------"
1916 * std::cerr <<
"Unknown exception!" << std::endl
1917 * <<
"Aborting!" << std::endl
1918 * <<
"----------------------------------------------------"
1926<a name=
"Results"></a><h1>Results</h1>
1929<a name=
"Outputoftheprogramandgraphicalvisualization"></a><h3>Output of the program and graphical visualization</h3>
1933If we
run the program as is, we get
this output
for the @f$32\times 32@f$
1934mesh we use (
for a total of 1024 cells with 1024 pressure degrees of
1935freedom since we use piecewise constants, and 2112 velocities because
1936the Raviart-Thomas element defines one degree per freedom per face and
1937there are @f$1024 + 32 = 1056@f$ faces
parallel to the @f$x@f$-axis and the same
1938number
parallel to the @f$y@f$-axis):
1941[ 66%] Built target step-20
1942Scanning dependencies of target
run
1943[100%] Run step-20 with Release configuration
1944Number of active cells: 1024
1945Total number of cells: 1365
1946Number of degrees of freedom: 3136 (2112+1024)
194724 CG Schur complement iterations to obtain convergence.
1948Errors: ||e_p||_L2 = 0.0445032, ||e_u||_L2 = 0.010826
1949[100%] Built target
run
1952The fact that the number of iterations is so small, of course, is due to
1953the good (but expensive!) preconditioner we have developed. To get
1954confidence in the solution, let us take a look at it. The following three
1955images show (from left to right) the x-velocity, the y-velocity, and the
1958<table style=
"width:60%" align=
"center">
1960 <td><img src=
"https://www.dealii.org/images/steps/developer/step-20.u_new.jpg" width=
"400" alt=
""></td>
1961 <td><img src=
"https://www.dealii.org/images/steps/developer/step-20.v_new.jpg" width=
"400" alt=
""></td>
1962 <td><img src=
"https://www.dealii.org/images/steps/developer/step-20.p_new.jpg" width=
"400" alt=
""></td>
1968Let us start with the pressure: it is highest at the left and lowest at the
1969right, so flow will be from left to right. In addition, though hardly visible
1970in the graph, we have chosen the pressure field such that the flow left-right
1971flow
first channels towards the
center and then outward again. Consequently,
1972the x-velocity has to increase to get the flow through the narrow part,
1973something that can easily be seen in the left image. The middle image
1974represents inward flow in y-direction at the left
end of the domain, and
1975outward flow in y-direction at the right
end of the domain.
1979As an additional remark, note how the x-velocity in the left image is only
1980continuous in x-direction, whereas the y-velocity is continuous in
1981y-direction. The flow fields are discontinuous in the other directions. This
1982very obviously reflects the continuity properties of the Raviart-Thomas
1983elements, which are, in fact, only in the space H(div) and not in the space
1984@f$H^1@f$. Finally, the pressure field is completely discontinuous, but
1985that should not surprise given that we have chosen <code>
FE_DGQ(0)</code> as
1986the finite element for that solution component.
1990<a name=
"Convergence"></a><h3>Convergence</h3>
1994The program offers two obvious places where playing and observing convergence
1995is in order: the degree of the finite elements used (passed to the constructor
1996of the <code>MixedLaplaceProblem</code> class from <code>main()</code>), and
1997the refinement
level (determined in
1998<code>MixedLaplaceProblem::make_grid_and_dofs</code>). What one can do is to
1999change these
values and observe the errors computed later on in the course of
2004If one does this, one finds the following pattern for the @f$L_2@f$ error
2005in the pressure variable:
2006<table align=
"center" class=
"doxtable">
2009 <th colspan=
"3" align=
"center">Finite element order</th>
2012 <th>Refinement
level</th>
2018 <th>0</th> <td>1.45344</td> <td>0.0831743</td> <td>0.0235186</td>
2021 <th>1</th> <td>0.715099</td> <td>0.0245341</td> <td>0.00293983</td>
2024 <th>2</th> <td>0.356383</td> <td>0.0063458</td> <td>0.000367478</td>
2027 <th>3</th> <td>0.178055</td> <td>0.00159944</td> <td>4.59349e-05</td>
2030 <th>4</th> <td>0.0890105</td> <td>0.000400669</td> <td>5.74184e-06</td>
2033 <th>5</th> <td>0.0445032</td> <td>0.000100218</td> <td>7.17799e-07</td>
2036 <th>6</th> <td>0.0222513</td> <td>2.50576e-05</td> <td>9.0164e-08</td>
2039 <th></th> <th>@f$O(h)@f$</th> <th>@f$O(h^2)@f$</th> <th>@f$O(h^3)@f$</th>
2043The theoretically expected convergence orders are very nicely reflected by the
2044experimentally observed ones indicated in the last row of the table.
2048One can make the same experiment with the @f$L_2@f$ error
2049in the velocity variables:
2050<table align=
"center" class=
"doxtable">
2053 <th colspan=
"3" align=
"center">Finite element order</th>
2056 <th>Refinement
level</th>
2062 <th>0</th> <td>0.367423</td> <td>0.127657</td> <td>5.10388e-14</td>
2065 <th>1</th> <td>0.175891</td> <td>0.0319142</td> <td>9.04414e-15</td>
2068 <th>2</th> <td>0.0869402</td> <td>0.00797856</td> <td>1.23723e-14</td>
2071 <th>3</th> <td>0.0433435</td> <td>0.00199464</td> <td>1.86345e-07</td>
2074 <th>4</th> <td>0.0216559</td> <td>0.00049866</td> <td>2.72566e-07</td>
2077 <th>5</th> <td>0.010826</td> <td>0.000124664</td> <td>3.57141e-07</td>
2080 <th>6</th> <td>0.00541274</td> <td>3.1166e-05</td> <td>4.46124e-07</td>
2083 <th></th> <td>@f$O(h)@f$</td> <td>@f$O(h^2)@f$</td> <td>@f$O(h^3)@f$</td>
2086The result concerning the convergence order is the same here.
2090<a name=
"extensions"></a>
2091<a name=
"Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2094<a name=
"Morerealisticpermeabilityfields"></a><h4>More realistic permeability fields</h4>
2097Realistic flow computations for ground water or oil reservoir simulations will
2098not use a
constant permeability. Here
's a first, rather simple way to change
2099this situation: we use a permeability that decays very rapidly away from a
2100central flowline until it hits a background value of 0.001. This is to mimic
2101the behavior of fluids in sandstone: in most of the domain, the sandstone is
2102homogeneous and, while permeable to fluids, not overly so; on the other stone,
2103the stone has cracked, or faulted, along one line, and the fluids flow much
2104easier along this large crack. Here is how we could implement something like
2109KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2110 std::vector<Tensor<2,dim> > &values) const
2112 AssertDimension (points.size(), values.size());
2114 for (unsigned int p=0; p<points.size(); ++p)
2118 const double distance_to_flowline
2119 = std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
2121 const double permeability = std::max(std::exp(-(distance_to_flowline*
2122 distance_to_flowline)
2126 for (unsigned int d=0; d<dim; ++d)
2127 values[p][d][d] = 1./permeability;
2131Remember that the function returns the inverse of the permeability tensor.
2135With a significantly higher mesh resolution, we can visualize this, here with
2138<table style="width:60%" align="center">
2140 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-wiggle.png" alt=""></td>
2141 <td><img src="https://www.dealii.org/images/steps/developer/step-20.v-wiggle.png" alt=""></td>
2145It is obvious how fluids flow essentially only along the middle line, and not
2150Another possibility would be to use a random permeability field. A simple way
2151to achieve this would be to scatter a number of centers around the domain and
2152then use a permeability field that is the sum of (negative) exponentials for
2153each of these centers. Flow would then try to hop from one center of high
2154permeability to the next one. This is an entirely unscientific attempt at
2155describing a random medium, but one possibility to implement this behavior
2156would look like this:
2159class KInverse : public TensorFunction<2,dim>
2164 virtual void value_list (const std::vector<Point<dim> > &points,
2165 std::vector<Tensor<2,dim> > &values) const;
2168 std::vector<Point<dim> > centers;
2173KInverse<dim>::KInverse ()
2175 const unsigned int N = 40;
2177 for (unsigned int i=0; i<N; ++i)
2178 for (unsigned int d=0; d<dim; ++d)
2179 centers[i][d] = 2.*rand()/RAND_MAX-1;
2185KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
2186 std::vector<Tensor<2,dim> > &values) const
2188 AssertDimension (points.size(), values.size());
2190 for (unsigned int p=0; p<points.size(); ++p)
2194 double permeability = 0;
2195 for (unsigned int i=0; i<centers.size(); ++i)
2196 permeability += std::exp(-(points[p] - centers[i]).norm_square() / (0.1 * 0.1));
2198 const double normalized_permeability
2199 = std::max(permeability, 0.005);
2201 for (unsigned int d=0; d<dim; ++d)
2202 values[p][d][d] = 1./normalized_permeability;
2207A piecewise constant interpolation of the diagonal elements of the
2208inverse of this tensor (i.e., of <code>normalized_permeability</code>)
2211<img src="https://www.dealii.org/images/steps/developer/step-20.k-random.png" alt="">
2214With a permeability field like this, we would get x-velocities and pressures as
2217<table style="width:60%" align="center">
2219 <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-random.png" alt=""></td>
2220 <td><img src="https://www.dealii.org/images/steps/developer/step-20.p-random.png" alt=""></td>
2224We will use these permeability fields again in @ref step_21 "step-21" and @ref step_43 "step-43".
2227<a name="Betterlinearsolvers"></a><h4>Better linear solvers</h4>
2230As mentioned in the introduction, the Schur complement solver used here is not
2231the best one conceivable (nor is it intended to be a particularly good
2232one). Better ones can be found in the literature and can be built using the
2233same block matrix techniques that were introduced here. We pick up on this
2234theme again in @ref step_22 "step-22", where we first build a Schur complement solver for the
2235Stokes equation as we did here, and then in the <a
2236href="step_22.html#improved-solver">Improved Solvers</a> section discuss better
2237ways based on solving the system as a whole but preconditioning based on
2238individual blocks. We will also come back to this in @ref step_43 "step-43".
2241<a name="PlainProg"></a>
2242<h1> The plain program</h1>
2243@include "step-20.cc"
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
__global__ void set(Number *val, const Number s, const size_type N)
#define AssertDimension(dim1, dim2)
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > 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())
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
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.
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()