811 *
const unsigned int = 0)
const override
822 * <a name=
"Pressureboundaryvalues"></a>
823 * <h4>Pressure boundary
values</h4>
827 * The next are pressure boundary
values. As mentioned in the introduction,
828 * we choose a linear pressure field:
832 *
class PressureBoundaryValues :
public Function<dim>
835 * PressureBoundaryValues()
840 *
const unsigned int = 0) const override
851 * <a name=
"Saturationboundaryvalues"></a>
852 * <h4>Saturation boundary
values</h4>
856 * Then we also need boundary
values on the inflow portions of the
857 * boundary. The question whether something is an inflow part is decided
858 * when assembling the right hand side, we only have to provide a functional
859 * description of the boundary
values. This is as explained in the
864 *
class SaturationBoundaryValues :
public Function<dim>
867 * SaturationBoundaryValues()
872 *
const unsigned int = 0) const override
886 * <a name=
"Initialdata"></a>
887 * <h4>Initial data</h4>
891 * Finally, we need
initial data. In reality, we only need
initial data
for
892 * the saturation, but we are lazy, so we will later, before the
first time
893 * step, simply
interpolate the entire solution
for the previous time step
894 * from a function that contains all vector components.
898 * We therefore simply create a function that returns
zero in all
899 * components. We
do that by simply forward every function to the
901 *
this program where we presently use the <code>InitialValues</code>
class?
902 * Because
this way it is simpler to later go back and choose a different
907 *
class InitialValues :
public Function<dim>
915 *
const unsigned int component = 0) const override
932 * <a name=
"Theinversepermeabilitytensor"></a>
933 * <h3>The inverse permeability tensor</h3>
937 * As announced in the introduction, we implement two different permeability
938 * tensor fields. Each of them we put into a
namespace of its own, so that
939 * it will be easy later to replace use of
one by the other in the code.
944 * <a name=
"Singlecurvingcrackpermeability"></a>
945 * <h4>Single curving crack permeability</h4>
949 * The
first function
for the permeability was the
one that models a single
950 * curving crack. It was already used at the
end of @ref step_20
"step-20", and its
951 * functional form is given in the introduction of the present tutorial
952 * program. As in some previous programs, we have to declare a (seemingly
953 * unnecessary)
default constructor of the KInverse
class to avoid warnings
954 * from some compilers:
957 *
namespace SingleCurvingCrack
974 *
for (
unsigned int p = 0; p < points.size(); ++p)
978 *
const double distance_to_flowline =
981 *
const double permeability =
986 *
for (
unsigned int d = 0;
d < dim; ++
d)
987 *
values[p][
d][
d] = 1. / permeability;
997 * <a name=
"Randommediumpermeability"></a>
998 * <h4>Random medium permeability</h4>
1002 * This function does as announced in the introduction, i.e. it creates an
1003 * overlay of exponentials at
random places. There is
one thing worth
1004 * considering
for this class. The issue centers around the problem that the
1005 *
class creates the centers of the exponentials using a
random function. If
1006 * we therefore created the centers each time we create an
object of the
1007 * present type, we would get a different list of centers each time. That
's
1008 * not what we expect from classes of this type: they should reliably
1009 * represent the same function.
1013 * The solution to this problem is to make the list of centers a static
1014 * member variable of this class, i.e. there exists exactly one such
1015 * variable for the entire program, rather than for each object of this
1016 * type. That's exactly what we are going to
do.
1020 * The next problem, however, is that we need a way to initialize
this
1021 * variable. Since
this variable is initialized at the beginning of the
1022 * program, we can
't use a regular member function for that since there may
1023 * not be an object of this type around at the time. The C++ standard
1024 * therefore says that only non-member and static member functions can be
1025 * used to initialize a static variable. We use the latter possibility by
1026 * defining a function <code>get_centers</code> that computes the list of
1027 * center points when called.
1031 * Note that this class works just fine in both 2d and 3d, with the only
1032 * difference being that we use more points in 3d: by experimenting we find
1033 * that we need more exponentials in 3d than in 2d (we have more ground to
1034 * cover, after all, if we want to keep the distance between centers roughly
1035 * equal), so we choose 40 in 2d and 100 in 3d. For any other dimension, the
1036 * function does presently not know what to do so simply throws an exception
1037 * indicating exactly this.
1040 * namespace RandomMedium
1042 * template <int dim>
1043 * class KInverse : public TensorFunction<2, dim>
1047 * : TensorFunction<2, dim>()
1051 * value_list(const std::vector<Point<dim>> &points,
1052 * std::vector<Tensor<2, dim>> & values) const override
1054 * Assert(points.size() == values.size(),
1055 * ExcDimensionMismatch(points.size(), values.size()));
1057 * for (unsigned int p = 0; p < points.size(); ++p)
1059 * values[p].clear();
1061 * double permeability = 0;
1062 * for (unsigned int i = 0; i < centers.size(); ++i)
1063 * permeability += std::exp(-(points[p] - centers[i]).norm_square() /
1066 * const double normalized_permeability =
1067 * std::min(std::max(permeability, 0.01), 4.);
1069 * for (unsigned int d = 0; d < dim; ++d)
1070 * values[p][d][d] = 1. / normalized_permeability;
1075 * static std::vector<Point<dim>> centers;
1077 * static std::vector<Point<dim>> get_centers()
1079 * const unsigned int N =
1080 * (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
1082 * std::vector<Point<dim>> centers_list(N);
1083 * for (unsigned int i = 0; i < N; ++i)
1084 * for (unsigned int d = 0; d < dim; ++d)
1085 * centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
1087 * return centers_list;
1093 * template <int dim>
1094 * std::vector<Point<dim>>
1095 * KInverse<dim>::centers = KInverse<dim>::get_centers();
1096 * } // namespace RandomMedium
1103 * <a name="Theinversemobilityandsaturationfunctions"></a>
1104 * <h3>The inverse mobility and saturation functions</h3>
1108 * There are two more pieces of data that we need to describe, namely the
1109 * inverse mobility function and the saturation curve. Their form is also
1110 * given in the introduction:
1113 * double mobility_inverse(const double S, const double viscosity)
1115 * return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
1118 * double fractional_flow(const double S, const double viscosity)
1120 * return S * S / (S * S + viscosity * (1 - S) * (1 - S));
1128 * <a name="Linearsolversandpreconditioners"></a>
1129 * <h3>Linear solvers and preconditioners</h3>
1133 * The linear solvers we use are also completely analogous to the ones used
1134 * in @ref step_20 "step-20". The following classes are therefore copied verbatim from
1135 * there. Note that the classes here are not only copied from
1136 * @ref step_20 "step-20", but also duplicate classes in deal.II. In a future version of this
1137 * example, they should be replaced by an efficient method, though. There is a
1138 * single change: if the size of a linear system is small, i.e. when the mesh
1139 * is very coarse, then it is sometimes not sufficient to set a maximum of
1140 * <code>src.size()</code> CG iterations before the solver in the
1141 * <code>vmult()</code> function converges. (This is, of course, a result of
1142 * numerical round-off, since we know that on paper, the CG method converges
1143 * in at most <code>src.size()</code> steps.) As a consequence, we set the
1144 * maximum number of iterations equal to the maximum of the size of the linear
1148 * template <class MatrixType>
1149 * class InverseMatrix : public Subscriptor
1152 * InverseMatrix(const MatrixType &m)
1156 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1158 * SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
1159 * 1e-8 * src.l2_norm());
1160 * SolverCG<Vector<double>> cg(solver_control);
1164 * cg.solve(*matrix, dst, src, PreconditionIdentity());
1168 * const SmartPointer<const MatrixType> matrix;
1173 * class SchurComplement : public Subscriptor
1176 * SchurComplement(const BlockSparseMatrix<double> & A,
1177 * const InverseMatrix<SparseMatrix<double>> &Minv)
1178 * : system_matrix(&A)
1179 * , m_inverse(&Minv)
1180 * , tmp1(A.block(0, 0).m())
1181 * , tmp2(A.block(0, 0).m())
1184 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1186 * system_matrix->block(0, 1).vmult(tmp1, src);
1187 * m_inverse->vmult(tmp2, tmp1);
1188 * system_matrix->block(1, 0).vmult(dst, tmp2);
1192 * const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
1193 * const SmartPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
1195 * mutable Vector<double> tmp1, tmp2;
1200 * class ApproximateSchurComplement : public Subscriptor
1203 * ApproximateSchurComplement(const BlockSparseMatrix<double> &A)
1204 * : system_matrix(&A)
1205 * , tmp1(A.block(0, 0).m())
1206 * , tmp2(A.block(0, 0).m())
1209 * void vmult(Vector<double> &dst, const Vector<double> &src) const
1211 * system_matrix->block(0, 1).vmult(tmp1, src);
1212 * system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
1213 * system_matrix->block(1, 0).vmult(dst, tmp2);
1217 * const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
1219 * mutable Vector<double> tmp1, tmp2;
1227 * <a name="codeTwoPhaseFlowProblemcodeclassimplementation"></a>
1228 * <h3><code>TwoPhaseFlowProblem</code> class implementation</h3>
1232 * Here now the implementation of the main class. Much of it is actually
1233 * copied from @ref step_20 "step-20", so we won't comment on it in much detail. You should
1234 *
try to get familiar with that program
first, then most of what is
1235 * happening here should be mostly clear.
1240 * <a name=
"TwoPhaseFlowProblemTwoPhaseFlowProblem"></a>
1241 * <h4>TwoPhaseFlowProblem::TwoPhaseFlowProblem</h4>
1245 * First
for the constructor. We use @f$RT_k \times DQ_k \times DQ_k@f$
1246 * spaces. For initializing the
DiscreteTime object, we don
't set the time
1247 * step size in the constructor because we don't have its
value yet.
1248 * The time step size is initially
set to
zero, but it will be computed
1249 * before it is needed to increment time, as described in a subsection of
1250 * the introduction. The time
object internally prevents itself from being
1251 * incremented when @f$dt = 0@f$, forcing us to
set a non-
zero desired size
for
1252 * @f$dt@f$ before advancing time.
1255 *
template <
int dim>
1256 * TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem(
const unsigned int degree)
1265 * , n_refinement_steps(5)
1275 * <a name=
"TwoPhaseFlowProblemmake_grid_and_dofs"></a>
1276 * <h4>TwoPhaseFlowProblem::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. It does
1282 * all the same things as in @ref step_20
"step-20", just now
for three components instead
1286 *
template <
int dim>
1287 *
void TwoPhaseFlowProblem<dim>::make_grid_and_dofs()
1292 * dof_handler.distribute_dofs(fe);
1295 *
const std::vector<types::global_dof_index> dofs_per_component =
1297 *
const unsigned int n_u = dofs_per_component[0],
1298 * n_p = dofs_per_component[dim],
1299 * n_s = dofs_per_component[dim + 1];
1301 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
1303 * <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1304 * <<
" (" << n_u <<
'+' << n_p <<
'+' << n_s <<
')' << std::endl
1307 *
const unsigned int n_couplings = dof_handler.max_couplings_between_dofs();
1309 * sparsity_pattern.reinit(3, 3);
1310 * sparsity_pattern.block(0, 0).reinit(n_u, n_u, n_couplings);
1311 * sparsity_pattern.block(1, 0).reinit(n_p, n_u, n_couplings);
1312 * sparsity_pattern.block(2, 0).reinit(n_s, n_u, n_couplings);
1313 * sparsity_pattern.block(0, 1).reinit(n_u, n_p, n_couplings);
1314 * sparsity_pattern.block(1, 1).reinit(n_p, n_p, n_couplings);
1315 * sparsity_pattern.block(2, 1).reinit(n_s, n_p, n_couplings);
1316 * sparsity_pattern.block(0, 2).reinit(n_u, n_s, n_couplings);
1317 * sparsity_pattern.block(1, 2).reinit(n_p, n_s, n_couplings);
1318 * sparsity_pattern.block(2, 2).reinit(n_s, n_s, n_couplings);
1320 * sparsity_pattern.collect_sizes();
1323 * sparsity_pattern.compress();
1326 * system_matrix.reinit(sparsity_pattern);
1329 * solution.reinit(3);
1330 * solution.block(0).reinit(n_u);
1331 * solution.block(1).reinit(n_p);
1332 * solution.block(2).reinit(n_s);
1333 * solution.collect_sizes();
1335 * old_solution.reinit(3);
1336 * old_solution.block(0).reinit(n_u);
1337 * old_solution.block(1).reinit(n_p);
1338 * old_solution.block(2).reinit(n_s);
1339 * old_solution.collect_sizes();
1341 * system_rhs.reinit(3);
1342 * system_rhs.block(0).reinit(n_u);
1343 * system_rhs.block(1).reinit(n_p);
1344 * system_rhs.block(2).reinit(n_s);
1345 * system_rhs.collect_sizes();
1352 * <a name=
"TwoPhaseFlowProblemassemble_system"></a>
1353 * <h4>TwoPhaseFlowProblem::assemble_system</h4>
1357 * This is the function that assembles the linear system, or at least
1358 * everything except the (1,3) block that depends on the still-unknown
1359 * velocity computed during this time step (we deal with this in
1360 * <code>assemble_rhs_S</code>). Much of it is again as in @ref step_20 "step-20", but we
1361 * have to deal with some nonlinearity this time. However, the top of the
1362 * function is pretty much as usual (note that we
set matrix and right hand
1363 * side to
zero at the beginning — something we didn't have to do for
1364 * stationary problems since there we use each
matrix object only once and
1365 * it is empty at the beginning anyway).
1369 * Note that in its present form, the function uses the permeability
1370 * implemented in the RandomMedium::KInverse class. Switching to the single
1371 * curved crack permeability function is as simple as just changing the
1375 * template <
int dim>
1376 *
void TwoPhaseFlowProblem<dim>::assemble_system()
1378 * system_matrix = 0;
1382 *
QGauss<dim - 1> face_quadrature_formula(degree + 2);
1385 * quadrature_formula,
1389 * face_quadrature_formula,
1394 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1396 *
const unsigned int n_q_points = quadrature_formula.size();
1397 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
1402 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1404 *
const PressureRightHandSide<dim> pressure_right_hand_side;
1405 *
const PressureBoundaryValues<dim> pressure_boundary_values;
1406 *
const RandomMedium::KInverse<dim> k_inverse;
1408 * std::vector<double> pressure_rhs_values(n_q_points);
1409 * std::vector<double> boundary_values(n_face_q_points);
1410 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1412 * std::vector<Vector<double>> old_solution_values(n_q_points,
1414 * std::vector<std::vector<Tensor<1, dim>>> old_solution_grads(
1421 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1423 * fe_values.reinit(cell);
1429 * Here
's the first significant difference: We have to get the values
1430 * of the saturation function of the previous time step at the
1431 * quadrature points. To this end, we can use the
1432 * FEValues::get_function_values (previously already used in @ref step_9 "step-9",
1433 * @ref step_14 "step-14" and @ref step_15 "step-15"), a function that takes a solution vector and
1434 * returns a list of function values at the quadrature points of the
1435 * present cell. In fact, it returns the complete vector-valued
1436 * solution at each quadrature point, i.e. not only the saturation but
1437 * also the velocities and pressure:
1440 * fe_values.get_function_values(old_solution, old_solution_values);
1444 * Then we also have to get the values of the pressure right hand side
1445 * and of the inverse permeability tensor at the quadrature points:
1448 * pressure_right_hand_side.value_list(fe_values.get_quadrature_points(),
1449 * pressure_rhs_values);
1450 * k_inverse.value_list(fe_values.get_quadrature_points(),
1451 * k_inverse_values);
1455 * With all this, we can now loop over all the quadrature points and
1456 * shape functions on this cell and assemble those parts of the matrix
1457 * and right hand side that we deal with in this function. The
1458 * individual terms in the contributions should be self-explanatory
1459 * given the explicit form of the bilinear form stated in the
1463 * for (unsigned int q = 0; q < n_q_points; ++q)
1464 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1466 * const double old_s = old_solution_values[q](dim + 1);
1468 * const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
1469 * const double div_phi_i_u = fe_values[velocities].divergence(i, q);
1470 * const double phi_i_p = fe_values[pressure].value(i, q);
1471 * const double phi_i_s = fe_values[saturation].value(i, q);
1473 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1475 * const Tensor<1, dim> phi_j_u =
1476 * fe_values[velocities].value(j, q);
1477 * const double div_phi_j_u =
1478 * fe_values[velocities].divergence(j, q);
1479 * const double phi_j_p = fe_values[pressure].value(j, q);
1480 * const double phi_j_s = fe_values[saturation].value(j, q);
1482 * local_matrix(i, j) +=
1483 * (phi_i_u * k_inverse_values[q] *
1484 * mobility_inverse(old_s, viscosity) * phi_j_u -
1485 * div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u +
1486 * phi_i_s * phi_j_s) *
1491 * (-phi_i_p * pressure_rhs_values[q]) * fe_values.JxW(q);
1497 * Next, we also have to deal with the pressure boundary values. This,
1498 * again is as in @ref step_20 "step-20":
1501 * for (const auto &face : cell->face_iterators())
1502 * if (face->at_boundary())
1504 * fe_face_values.reinit(cell, face);
1506 * pressure_boundary_values.value_list(
1507 * fe_face_values.get_quadrature_points(), boundary_values);
1509 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1510 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1512 * const Tensor<1, dim> phi_i_u =
1513 * fe_face_values[velocities].value(i, q);
1516 * -(phi_i_u * fe_face_values.normal_vector(q) *
1517 * boundary_values[q] * fe_face_values.JxW(q));
1523 * The final step in the loop over all cells is to transfer local
1524 * contributions into the global matrix and right hand side vector:
1527 * cell->get_dof_indices(local_dof_indices);
1528 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1529 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1530 * system_matrix.add(local_dof_indices[i],
1531 * local_dof_indices[j],
1532 * local_matrix(i, j));
1534 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1535 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1542 * So much for assembly of matrix and right hand side. Note that we do not
1543 * have to interpolate and apply boundary values since they have all been
1544 * taken care of in the weak form already.
1552 * <a name="TwoPhaseFlowProblemassemble_rhs_S"></a>
1553 * <h4>TwoPhaseFlowProblem::assemble_rhs_S</h4>
1557 * As explained in the introduction, we can only evaluate the right hand
1558 * side of the saturation equation once the velocity has been computed. We
1559 * therefore have this separate function to this end.
1562 * template <int dim>
1563 * void TwoPhaseFlowProblem<dim>::assemble_rhs_S()
1565 * QGauss<dim> quadrature_formula(degree + 2);
1566 * QGauss<dim - 1> face_quadrature_formula(degree + 2);
1567 * FEValues<dim> fe_values(fe,
1568 * quadrature_formula,
1569 * update_values | update_gradients |
1570 * update_quadrature_points | update_JxW_values);
1571 * FEFaceValues<dim> fe_face_values(fe,
1572 * face_quadrature_formula,
1573 * update_values | update_normal_vectors |
1574 * update_quadrature_points |
1575 * update_JxW_values);
1576 * FEFaceValues<dim> fe_face_values_neighbor(fe,
1577 * face_quadrature_formula,
1580 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1581 * const unsigned int n_q_points = quadrature_formula.size();
1582 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1584 * Vector<double> local_rhs(dofs_per_cell);
1586 * std::vector<Vector<double>> old_solution_values(n_q_points,
1587 * Vector<double>(dim + 2));
1588 * std::vector<Vector<double>> old_solution_values_face(n_face_q_points,
1589 * Vector<double>(dim +
1591 * std::vector<Vector<double>> old_solution_values_face_neighbor(
1592 * n_face_q_points, Vector<double>(dim + 2));
1593 * std::vector<Vector<double>> present_solution_values(n_q_points,
1594 * Vector<double>(dim +
1596 * std::vector<Vector<double>> present_solution_values_face(
1597 * n_face_q_points, Vector<double>(dim + 2));
1599 * std::vector<double> neighbor_saturation(n_face_q_points);
1600 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1602 * SaturationBoundaryValues<dim> saturation_boundary_values;
1604 * const FEValuesExtractors::Scalar saturation(dim + 1);
1606 * for (const auto &cell : dof_handler.active_cell_iterators())
1609 * fe_values.reinit(cell);
1611 * fe_values.get_function_values(old_solution, old_solution_values);
1612 * fe_values.get_function_values(solution, present_solution_values);
1616 * First for the cell terms. These are, following the formulas in the
1617 * introduction, @f$(S^n,\sigma)-(F(S^n) \mathbf{v}^{n+1},\nabla
1618 * \sigma)@f$, where @f$\sigma@f$ is the saturation component of the test
1622 * for (unsigned int q = 0; q < n_q_points; ++q)
1623 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1625 * const double old_s = old_solution_values[q](dim + 1);
1626 * Tensor<1, dim> present_u;
1627 * for (unsigned int d = 0; d < dim; ++d)
1628 * present_u[d] = present_solution_values[q](d);
1630 * const double phi_i_s = fe_values[saturation].value(i, q);
1631 * const Tensor<1, dim> grad_phi_i_s =
1632 * fe_values[saturation].gradient(i, q);
1635 * (time.get_next_step_size() * fractional_flow(old_s, viscosity) *
1636 * present_u * grad_phi_i_s +
1637 * old_s * phi_i_s) *
1643 * Secondly, we have to deal with the flux parts on the face
1644 * boundaries. This was a bit more involved because we first have to
1645 * determine which are the influx and outflux parts of the cell
1646 * boundary. If we have an influx boundary, we need to evaluate the
1647 * saturation on the other side of the face (or the boundary values,
1648 * if we are at the boundary of the domain).
1652 * All this is a bit tricky, but has been explained in some detail
1653 * already in @ref step_9 "step-9". Take a look there how this is supposed to work!
1656 * for (const auto face_no : cell->face_indices())
1658 * fe_face_values.reinit(cell, face_no);
1660 * fe_face_values.get_function_values(old_solution,
1661 * old_solution_values_face);
1662 * fe_face_values.get_function_values(solution,
1663 * present_solution_values_face);
1665 * if (cell->at_boundary(face_no))
1666 * saturation_boundary_values.value_list(
1667 * fe_face_values.get_quadrature_points(), neighbor_saturation);
1670 * const auto neighbor = cell->neighbor(face_no);
1671 * const unsigned int neighbor_face =
1672 * cell->neighbor_of_neighbor(face_no);
1674 * fe_face_values_neighbor.reinit(neighbor, neighbor_face);
1676 * fe_face_values_neighbor.get_function_values(
1677 * old_solution, old_solution_values_face_neighbor);
1679 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1680 * neighbor_saturation[q] =
1681 * old_solution_values_face_neighbor[q](dim + 1);
1685 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1687 * Tensor<1, dim> present_u_face;
1688 * for (unsigned int d = 0; d < dim; ++d)
1689 * present_u_face[d] = present_solution_values_face[q](d);
1691 * const double normal_flux =
1692 * present_u_face * fe_face_values.normal_vector(q);
1694 * const bool is_outflow_q_point = (normal_flux >= 0);
1696 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1698 * time.get_next_step_size() * normal_flux *
1699 * fractional_flow((is_outflow_q_point == true ?
1700 * old_solution_values_face[q](dim + 1) :
1701 * neighbor_saturation[q]),
1703 * fe_face_values[saturation].value(i, q) *
1704 * fe_face_values.JxW(q);
1708 * cell->get_dof_indices(local_dof_indices);
1709 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1710 * system_rhs(local_dof_indices[i]) += local_rhs(i);
1719 * <a name="TwoPhaseFlowProblemsolve"></a>
1720 * <h4>TwoPhaseFlowProblem::solve</h4>
1724 * After all these preparations, we finally solve the linear system for
1725 * velocity and pressure in the same way as in @ref step_20 "step-20". After that, we have
1726 * to deal with the saturation equation (see below):
1729 * template <int dim>
1730 * void TwoPhaseFlowProblem<dim>::solve()
1732 * const InverseMatrix<SparseMatrix<double>> m_inverse(
1733 * system_matrix.block(0, 0));
1734 * Vector<double> tmp(solution.block(0).size());
1735 * Vector<double> schur_rhs(solution.block(1).size());
1736 * Vector<double> tmp2(solution.block(2).size());
1741 * First the pressure, using the pressure Schur complement of the first
1746 * m_inverse.vmult(tmp, system_rhs.block(0));
1747 * system_matrix.block(1, 0).vmult(schur_rhs, tmp);
1748 * schur_rhs -= system_rhs.block(1);
1751 * SchurComplement schur_complement(system_matrix, m_inverse);
1753 * ApproximateSchurComplement approximate_schur_complement(system_matrix);
1755 * InverseMatrix<ApproximateSchurComplement> preconditioner(
1756 * approximate_schur_complement);
1759 * SolverControl solver_control(solution.block(1).size(),
1760 * 1e-12 * schur_rhs.l2_norm());
1761 * SolverCG<Vector<double>> cg(solver_control);
1763 * cg.solve(schur_complement, solution.block(1), schur_rhs, preconditioner);
1765 * std::cout << " " << solver_control.last_step()
1766 * << " CG Schur complement iterations for pressure." << std::endl;
1775 * system_matrix.block(0, 1).vmult(tmp, solution.block(1));
1777 * tmp += system_rhs.block(0);
1779 * m_inverse.vmult(solution.block(0), tmp);
1784 * Finally, we have to take care of the saturation equation. The first
1785 * business we have here is to determine the time step using the formula
1786 * in the introduction. Knowing the shape of our domain and that we
1787 * created the mesh by regular subdivision of cells, we can compute the
1788 * diameter of each of our cells quite easily (in fact we use the linear
1789 * extensions in coordinate directions of the cells, not the
1790 * diameter). Note that we will learn a more general way to do this in
1791 * @ref step_24 "step-24", where we use the GridTools::minimal_cell_diameter function.
1795 * The maximal velocity we compute using a helper function to compute the
1796 * maximal velocity defined below, and with all this we can evaluate our
1797 * new time step length. We use the method
1798 * DiscreteTime::set_desired_next_time_step() to suggest the new
1799 * calculated value of the time step to the DiscreteTime object. In most
1800 * cases, the time object uses the exact provided value to increment time.
1801 * It some case, the step size may be modified further by the time object.
1802 * For example, if the calculated time increment overshoots the end time,
1803 * it is truncated accordingly.
1806 * time.set_desired_next_step_size(std::pow(0.5, double(n_refinement_steps)) /
1807 * get_maximal_velocity());
1811 * The next step is to assemble the right hand side, and then to pass
1812 * everything on for solution. At the end, we project back saturations
1813 * onto the physically reasonable range:
1818 * SolverControl solver_control(system_matrix.block(2, 2).m(),
1819 * 1e-8 * system_rhs.block(2).l2_norm());
1820 * SolverCG<Vector<double>> cg(solver_control);
1821 * cg.solve(system_matrix.block(2, 2),
1822 * solution.block(2),
1823 * system_rhs.block(2),
1824 * PreconditionIdentity());
1826 * project_back_saturation();
1828 * std::cout << " " << solver_control.last_step()
1829 * << " CG iterations for saturation." << std::endl;
1833 * old_solution = solution;
1840 * <a name="TwoPhaseFlowProblemoutput_results"></a>
1841 * <h4>TwoPhaseFlowProblem::output_results</h4>
1845 * There is nothing surprising here. Since the program will do a lot of time
1846 * steps, we create an output file only every fifth time step and skip all
1847 * other time steps at the top of the file already.
1851 * When creating file names for output close to the bottom of the function,
1852 * we convert the number of the time step to a string representation that
1853 * is padded by leading zeros to four digits. We do this because this way
1854 * all output file names have the same length, and consequently sort well
1855 * when creating a directory listing.
1858 * template <int dim>
1859 * void TwoPhaseFlowProblem<dim>::output_results() const
1861 * if (time.get_step_number() % 5 != 0)
1864 * std::vector<std::string> solution_names;
1868 * solution_names = {"u", "v", "p", "S"};
1872 * solution_names = {"u", "v", "w", "p", "S"};
1876 * Assert(false, ExcNotImplemented());
1879 * DataOut<dim> data_out;
1881 * data_out.attach_dof_handler(dof_handler);
1882 * data_out.add_data_vector(solution, solution_names);
1884 * data_out.build_patches(degree + 1);
1886 * std::ofstream output("solution-" +
1887 * Utilities::int_to_string(time.get_step_number(), 4) +
1889 * data_out.write_vtk(output);
1897 * <a name="TwoPhaseFlowProblemproject_back_saturation"></a>
1898 * <h4>TwoPhaseFlowProblem::project_back_saturation</h4>
1902 * In this function, we simply run over all saturation degrees of freedom
1903 * and make sure that if they should have left the physically reasonable
1904 * range, that they be reset to the interval @f$[0,1]@f$. To do this, we only
1905 * have to loop over all saturation components of the solution vector; these
1906 * are stored in the block 2 (block 0 are the velocities, block 1 are the
1911 * It may be instructive to note that this function almost never triggers
1912 * when the time step is chosen as mentioned in the introduction. However,
1913 * if we choose the timestep only slightly larger, we get plenty of values
1914 * outside the proper range. Strictly speaking, the function is therefore
1915 * unnecessary if we choose the time step small enough. In a sense, the
1916 * function is therefore only a safety device to avoid situations where our
1917 * entire solution becomes unphysical because individual degrees of freedom
1918 * have become unphysical a few time steps earlier.
1921 * template <int dim>
1922 * void TwoPhaseFlowProblem<dim>::project_back_saturation()
1924 * for (unsigned int i = 0; i < solution.block(2).size(); ++i)
1925 * if (solution.block(2)(i) < 0)
1926 * solution.block(2)(i) = 0;
1927 * else if (solution.block(2)(i) > 1)
1928 * solution.block(2)(i) = 1;
1935 * <a name="TwoPhaseFlowProblemget_maximal_velocity"></a>
1936 * <h4>TwoPhaseFlowProblem::get_maximal_velocity</h4>
1940 * The following function is used in determining the maximal allowable time
1941 * step. What it does is to loop over all quadrature points in the domain
1942 * and find what the maximal magnitude of the velocity is.
1945 * template <int dim>
1946 * double TwoPhaseFlowProblem<dim>::get_maximal_velocity() const
1948 * QGauss<dim> quadrature_formula(degree + 2);
1949 * const unsigned int n_q_points = quadrature_formula.size();
1951 * FEValues<dim> fe_values(fe, quadrature_formula, update_values);
1952 * std::vector<Vector<double>> solution_values(n_q_points,
1953 * Vector<double>(dim + 2));
1954 * double max_velocity = 0;
1956 * for (const auto &cell : dof_handler.active_cell_iterators())
1958 * fe_values.reinit(cell);
1959 * fe_values.get_function_values(solution, solution_values);
1961 * for (unsigned int q = 0; q < n_q_points; ++q)
1963 * Tensor<1, dim> velocity;
1964 * for (unsigned int i = 0; i < dim; ++i)
1965 * velocity[i] = solution_values[q](i);
1967 * max_velocity = std::max(max_velocity, velocity.norm());
1971 * return max_velocity;
1978 * <a name="TwoPhaseFlowProblemrun"></a>
1979 * <h4>TwoPhaseFlowProblem::run</h4>
1983 * This is the final function of our main class. Its brevity speaks for
1984 * itself. There are only two points worth noting: First, the function
1985 * projects the initial values onto the finite element space at the
1986 * beginning; the VectorTools::project function doing this requires an
1987 * argument indicating the hanging node constraints. We have none in this
1988 * program (we compute on a uniformly refined mesh), but the function
1989 * requires the argument anyway, of course. So we have to create a
1990 * constraint object. In its original state, constraint objects are
1991 * unsorted, and have to be sorted (using the AffineConstraints::close
1992 * function) before they can be used. This is what we do here, and which is
1999 * The
second point worth mentioning is that we only compute the length of
2000 * the present time step in the middle of solving the linear system
2001 * corresponding to each time step. We can therefore output the present
2002 * time of a time step only at the
end of the time step.
2004 * inside the
loop. Since we are reporting the time and dt after we
2005 * increment it, we have to
call the method
2007 *
DiscreteTime::get_next_step_size(). After many steps, when the simulation
2008 * reaches the
end time, the last dt is chosen by the
DiscreteTime class in
2009 * such a way that the last step finishes exactly at the
end time.
2012 * template <
int dim>
2013 *
void TwoPhaseFlowProblem<dim>::
run()
2015 * make_grid_and_dofs();
2019 * constraints.
close();
2024 * InitialValues<dim>(),
2030 * std::cout <<
"Timestep " << time.get_step_number() + 1 << std::endl;
2032 * assemble_system();
2038 * time.advance_time();
2039 * std::cout <<
" Now at t=" << time.get_current_time()
2040 * <<
", dt=" << time.get_previous_step_size() <<
'.'
2044 *
while (time.is_at_end() ==
false);
2052 * <a name=
"Thecodemaincodefunction"></a>
2053 * <h3>The <code>main</code> function</h3>
2057 * That
's it. In the main function, we pass the degree of the finite element
2058 * space to the constructor of the TwoPhaseFlowProblem object. Here, we use
2059 * zero-th degree elements, i.e. @f$RT_0\times DQ_0 \times DQ_0@f$. The rest is as
2060 * in all the other programs.
2067 * using namespace Step21;
2069 * TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2070 * two_phase_flow_problem.run();
2072 * catch (std::exception &exc)
2074 * std::cerr << std::endl
2076 * << "----------------------------------------------------"
2078 * std::cerr << "Exception on processing: " << std::endl
2079 * << exc.what() << std::endl
2080 * << "Aborting!" << std::endl
2081 * << "----------------------------------------------------"
2088 * std::cerr << std::endl
2090 * << "----------------------------------------------------"
2092 * std::cerr << "Unknown exception!" << std::endl
2093 * << "Aborting!" << std::endl
2094 * << "----------------------------------------------------"
2102<a name="Results"></a><h1>Results</h1>
2105The code as presented here does not actually compute the results
2106found on the web page. The reason is, that even on a decent
2107computer it runs more than a day. If you want to reproduce these
2108results, modify the end time of the DiscreteTime object to `250` within the
2109constructor of TwoPhaseFlowProblem.
2111If we run the program, we get the following kind of output:
2113Number of active cells: 1024
2114Number of degrees of freedom: 4160 (2112+1024+1024)
2117 22 CG Schur complement iterations for pressure.
2118 1 CG iterations for saturation.
2119 Now at t=0.0326742, dt=0.0326742.
2122 17 CG Schur complement iterations for pressure.
2123 1 CG iterations for saturation.
2124 Now at t=0.0653816, dt=0.0327074.
2127 17 CG Schur complement iterations for pressure.
2128 1 CG iterations for saturation.
2129 Now at t=0.0980651, dt=0.0326836.
2133As we can see, the time step is pretty much constant right from the start,
2134which indicates that the velocities in the domain are not strongly dependent
2135on changes in saturation, although they certainly are through the factor
2136@f$\lambda(S)@f$ in the pressure equation.
2138Our second observation is that the number of CG iterations needed to solve the
2139pressure Schur complement equation drops from 22 to 17 between the first and
2140the second time step (in fact, it remains around 17 for the rest of the
2141computations). The reason is actually simple: Before we solve for the pressure
2142during a time step, we don't reset the <code>solution</code> variable to
2143zero. The pressure (and the other variables) therefore have the previous time
2144step
's values at the time we get into the CG solver. Since the velocities and
2145pressures don't change very much as computations progress, the previous time
2146step
's pressure is actually a good initial guess for this time step's
2147pressure. Consequently, the number of iterations we need once we have computed
2148the pressure once is significantly reduced.
2150The
final observation concerns the number of iterations needed to solve
for
2151the saturation, i.e.
one. This shouldn
't surprise us too much: the matrix we
2152have to solve with is the mass matrix. However, this is the mass matrix for
2153the @f$DGQ_0@f$ element of piecewise constants where no element couples with the
2154degrees of freedom on neighboring cells. The matrix is therefore a diagonal
2155one, and it is clear that we should be able to invert this matrix in a single
2159With all this, here are a few movies that show how the saturation progresses
2160over time. First, this is for the single crack model, as implemented in the
2161<code>SingleCurvingCrack::KInverse</code> class:
2163<img src="https://www.dealii.org/images/steps/developer/step-21.centerline.gif" alt="">
2165As can be seen, the water rich fluid snakes its way mostly along the
2166high-permeability zone in the middle of the domain, whereas the rest of the
2167domain is mostly impermeable. This and the next movie are generated using
2168<code>n_refinement_steps=7</code>, leading to a @f$128\times 128@f$ mesh with some
216916,000 cells and about 66,000 unknowns in total.
2172The second movie shows the saturation for the random medium model of class
2173<code>RandomMedium::KInverse</code>, where we have randomly distributed
2174centers of high permeability and fluid hops from one of these zones to
2177<img src="https://www.dealii.org/images/steps/developer/step-21.random2d.gif" alt="">
2180Finally, here is the same situation in three space dimensions, on a mesh with
2181<code>n_refinement_steps=5</code>, which produces a mesh of some 32,000 cells
2182and 167,000 degrees of freedom:
2184<img src="https://www.dealii.org/images/steps/developer/step-21.random3d.gif" alt="">
2186To repeat these computations, all you have to do is to change the line
2188 TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
2190in the main function to
2192 TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
2194The visualization uses a cloud technique, where the saturation is indicated by
2195colored but transparent clouds for each cell. This way, one can also see
2196somewhat what happens deep inside the domain. A different way of visualizing
2197would have been to show isosurfaces of the saturation evolving over
2198time. There are techniques to plot isosurfaces transparently, so that one can
2199see several of them at the same time like the layers of an onion.
2201So why don't we show such isosurfaces? The problem lies in the way isosurfaces
2202are computed: they require that the field to be visualized is continuous, so
2203that the isosurfaces can be generated by following contours at least across a
2204single cell. However, our saturation field is piecewise constant and
2205discontinuous. If we wanted to plot an isosurface
for a saturation @f$S=0.5@f$,
2206chances would be that there is no single
point in the domain where that
2207saturation is actually attained. If we had to define isosurfaces in that
2208context at all, we would have to take the interfaces between cells, where
one
2209of the two adjacent cells has a saturation greater than and the other cell a
2210saturation less than 0.5. However, it appears that most visualization programs
2211are not equipped to
do this kind of transformation.
2214<a name=
"extensions"></a>
2215<a name=
"Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
2218There are a number of areas where
this program can be improved. Three of them
2219are listed below. All of them are, in fact, addressed in a tutorial program
2220that forms the continuation of the current
one: @ref step_43
"step-43".
2223<a name=
"Solvers"></a><h4>Solvers</h4>
2226At present, the program is not particularly fast: the 2
d random medium
2227computation took about a day
for the 1,000 or so time steps. The corresponding
22283
d computation took almost two days
for 800 time steps. The reason why it
2229isn
't faster than this is twofold. First, we rebuild the entire matrix in
2230every time step, although some parts such as the @f$B@f$, @f$B^T@f$, and @f$M^S@f$ blocks
2233Second, we could do a lot better with the solver and
2234preconditioners. Presently, we solve the Schur complement @f$B^TM^u(S)^{-1}B@f$
2235with a CG method, using @f$[B^T (\textrm{diag}(M^u(S)))^{-1} B]^{-1}@f$ as a
2236preconditioner. Applying this preconditioner is expensive, since it involves
2237solving a linear system each time. This may have been appropriate for @ref
2238step_20 "step-20", where we have to solve the entire problem only
2239once. However, here we have to solve it hundreds of times, and in such cases
2240it is worth considering a preconditioner that is more expensive to set up the
2241first time, but cheaper to apply later on.
2243One possibility would be to realize that the matrix we use as preconditioner,
2244@f$B^T (\textrm{diag}(M^u(S)))^{-1} B@f$ is still sparse, and symmetric on top of
2245that. If one looks at the flow field evolve over time, we also see that while
2246@f$S@f$ changes significantly over time, the pressure hardly does and consequently
2247@f$B^T (\textrm{diag}(M^u(S)))^{-1} B \approx B^T (\textrm{diag}(M^u(S^0)))^{-1}
2248B@f$. In other words, the matrix for the first time step should be a good
2249preconditioner also for all later time steps. With a bit of
2250back-and-forthing, it isn't hard to actually get a representation of it as a
2252a sparse incomplete Cholesky decomposition. To form this decomposition is
2253expensive, but we have to do it only once in the
first time step, and can then
2254use it as a cheap preconditioner in the future. We could do better even by
2256a complete decomposition of the
matrix, which should yield an even better
2259Finally, why use the approximation @f$B^
T (\textrm{diag}(M^u(S)))^{-1} B@f$ to
2260precondition @f$B^
T M^u(S)^{-1} B@f$? The latter
matrix, after all, is the mixed
2261form of the Laplace
operator on the pressure space,
for which we use linear
2262elements. We could therefore build a separate
matrix @f$A^p@f$ on the side that
2263directly corresponds to the non-mixed formulation of the Laplacian,
for
2264example
using the bilinear form @f$(\mathbf{K}\lambda(S^n) \nabla
2265\varphi_i,\nabla\varphi_j)@f$. We could then form an incomplete or complete
2266decomposition of
this non-mixed
matrix and use it as a preconditioner of the
2269Using such techniques, it can reasonably be expected that the solution process
2270will be faster by at least an order of magnitude.
2273<a name=
"Timestepping"></a><h4>Time stepping</h4>
2276In the introduction we have identified the time step restriction
2278 \triangle t_{n+1} \le \frac h{|\mathbf{u}^{n+1}(\mathbf{x})|}
2280that has to hold globally, i.e.
for all @f$\mathbf x@f$. After discretization, we
2281satisfy it by choosing
2283 \triangle t_{n+1} = \frac {\min_K h_K}{\max_{\mathbf{x}}|\mathbf{u}^{n+1}(\mathbf{x})|}.
2286This restriction on the time step is somewhat annoying: the finer we make the
2287mesh the smaller the time step; in other words, we get punished twice: each
2288time step is more expensive to solve and we have to
do more time steps.
2290This is particularly annoying since the majority of the additional work is
2291spent solving the implicit part of the equations, i.e. the pressure-velocity
2292system, whereas it is the hyperbolic transport equation
for the saturation
2293that imposes the time step restriction.
2295To avoid
this bottleneck, people have invented a number of approaches. For
2296example, they may only re-compute the pressure-velocity field every few time
2297steps (or,
if you want, use different time step sizes
for the
2298pressure/velocity and saturation equations). This keeps the time step
2299restriction on the cheap
explicit part
while it makes the solution of the
2300implicit part less frequent. Experiments in
this direction are
2301certainly worthwhile;
one starting
point for such an approach is the paper by
2302Zhangxin Chen, Guanren Huan and Baoyan Li: <i>An improved IMPES method
for
2303two-phase flow in porous media</i>, Transport in Porous Media, 54 (2004),
2304pp. 361—376. There are certainly many other papers on
this topic as well, but
2305this one happened to land on our desk a
while back.
2309<a name=
"Adaptivity"></a><h4>Adaptivity</h4>
2312Adaptivity would also clearly help. Looking at the movies,
one clearly sees
2313that most of the action is confined to a relatively small part of the domain
2314(
this particularly obvious
for the saturation, but also holds
for the
2315velocities and pressures). Adaptivity can therefore be expected to keep the
2316necessary number of degrees of freedom low, or alternatively increase the
2319On the other hand, adaptivity
for time dependent problems is not a trivial
2320thing: we would have to change the mesh every few time steps, and we would
2321have to transport our present solution to the next mesh every time we change
2323insurmountable obstacles, but they
do require some additional coding and more
2324than we felt comfortable was worth packing into
this tutorial program.
2327<a name=
"PlainProg"></a>
2328<h1> The plain program</h1>
2329@include
"step-21.cc"
double get_previous_step_size() const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &return_value) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
@ 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.
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &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)
Expression fabs(const Expression &x)
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)
static const types::blas_int zero
@ matrix
Contents is actually a matrix.
static const 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)
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 typename identity< Iterator >::type &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)
int(&) functions(const void *v1, const void *v2)
::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