795 *
class SaturationBoundaryValues :
public Function<dim>
798 * SaturationBoundaryValues()
803 *
const unsigned int component = 0)
const override;
811 *
const unsigned int )
const
821 *
class SaturationInitialValues :
public Function<dim>
824 * SaturationInitialValues()
829 *
const unsigned int component = 0)
const override;
839 *
const unsigned int )
const
846 *
void SaturationInitialValues<dim>::vector_value(
const Point<dim> &p,
857 * <a name=
"Permeabilitymodels"></a>
858 * <h3>Permeability models</h3>
862 * In
this tutorial, we still use the two permeability models previously
863 * used in @ref step_21
"step-21" so we again refrain from commenting in detail about them.
866 *
namespace SingleCurvingCrack
883 *
void KInverse<dim>::value_list(
const std::vector<
Point<dim>> &points,
886 *
Assert(points.size() == values.size(),
889 *
for (
unsigned int p = 0; p < points.size(); ++p)
893 *
const double distance_to_flowline =
894 *
std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
896 *
const double permeability =
897 *
std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
901 *
for (
unsigned int d = 0;
d < dim; ++
d)
902 * values[p][
d][
d] = 1. / permeability;
908 *
namespace RandomMedium
923 *
static std::vector<Point<dim>> centers;
929 * std::vector<Point<dim>> KInverse<dim>::centers = []() {
930 *
const unsigned int N =
933 * std::vector<Point<dim>> centers_list(
N);
934 *
for (
unsigned int i = 0; i <
N; ++i)
935 *
for (
unsigned int d = 0;
d < dim; ++
d)
936 * centers_list[i][
d] =
static_cast<double>(rand()) / RAND_MAX;
938 *
return centers_list;
944 *
void KInverse<dim>::value_list(
const std::vector<
Point<dim>> &points,
949 *
for (
unsigned int p = 0; p < points.size(); ++p)
953 *
double permeability = 0;
954 *
for (
unsigned int i = 0; i < centers.size(); ++i)
956 * std::exp(-(points[p] - centers[i]).norm_square() / (0.05 * 0.05));
958 *
const double normalized_permeability =
961 *
for (
unsigned int d = 0;
d < dim; ++
d)
962 * values[p][
d][
d] = 1. / normalized_permeability;
971 * <a name=
"Physicalquantities"></a>
972 * <h3>Physical quantities</h3>
976 * The implementations of all the physical quantities such as total mobility
977 * @f$\lambda_t@f$ and fractional flow of water @f$F@f$ are taken from @ref step_21
"step-21" so
978 * again we don
't have do any comment about them. Compared to @ref step_21 "step-21" we
979 * have added checks that the saturation passed to these functions is in
980 * fact within the physically valid range. Furthermore, given that the
981 * wetting phase moves at speed @f$\mathbf u F'(S)@f$ it is clear that @f$F
'(S)@f$
982 * must be greater or equal to zero, so we assert that as well to make sure
983 * that our calculations to get at the formula for the derivative made
987 * double mobility_inverse(const double S, const double viscosity)
989 * return 1.0 / (1.0 / viscosity * S * S + (1 - S) * (1 - S));
993 * double fractional_flow(const double S, const double viscosity)
995 * Assert((S >= 0) && (S <= 1),
996 * ExcMessage("Saturation is outside its physically valid range."));
998 * return S * S / (S * S + viscosity * (1 - S) * (1 - S));
1002 * double fractional_flow_derivative(const double S, const double viscosity)
1004 * Assert((S >= 0) && (S <= 1),
1005 * ExcMessage("Saturation is outside its physically valid range."));
1007 * const double temp = (S * S + viscosity * (1 - S) * (1 - S));
1009 * const double numerator =
1010 * 2.0 * S * temp - S * S * (2.0 * S - 2.0 * viscosity * (1 - S));
1011 * const double denominator = std::pow(temp, 2.0);
1013 * const double F_prime = numerator / denominator;
1015 * Assert(F_prime >= 0, ExcInternalError());
1024 * <a name="Helperclassesforsolversandpreconditioners"></a>
1025 * <h3>Helper classes for solvers and preconditioners</h3>
1029 * In this first part we define a number of classes that we need in the
1030 * construction of linear solvers and preconditioners. This part is
1031 * essentially the same as that used in @ref step_31 "step-31". The only difference is that
1032 * the original variable name stokes_matrix is replaced by another name
1033 * darcy_matrix to match our problem.
1036 * namespace LinearSolvers
1038 * template <class MatrixType, class PreconditionerType>
1039 * class InverseMatrix : public Subscriptor
1042 * InverseMatrix(const MatrixType & m,
1043 * const PreconditionerType &preconditioner);
1046 * template <typename VectorType>
1047 * void vmult(VectorType &dst, const VectorType &src) const;
1050 * const SmartPointer<const MatrixType> matrix;
1051 * const PreconditionerType & preconditioner;
1055 * template <class MatrixType, class PreconditionerType>
1056 * InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
1057 * const MatrixType & m,
1058 * const PreconditionerType &preconditioner)
1060 * , preconditioner(preconditioner)
1065 * template <class MatrixType, class PreconditionerType>
1066 * template <typename VectorType>
1067 * void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1069 * const VectorType &src) const
1071 * SolverControl solver_control(src.size(), 1e-7 * src.l2_norm());
1072 * SolverCG<VectorType> cg(solver_control);
1078 * cg.solve(*matrix, dst, src, preconditioner);
1080 * catch (std::exception &e)
1082 * Assert(false, ExcMessage(e.what()));
1086 * template <class PreconditionerTypeA, class PreconditionerTypeMp>
1087 * class BlockSchurPreconditioner : public Subscriptor
1090 * BlockSchurPreconditioner(
1091 * const TrilinosWrappers::BlockSparseMatrix &S,
1092 * const InverseMatrix<TrilinosWrappers::SparseMatrix,
1093 * PreconditionerTypeMp> &Mpinv,
1094 * const PreconditionerTypeA & Apreconditioner);
1096 * void vmult(TrilinosWrappers::MPI::BlockVector & dst,
1097 * const TrilinosWrappers::MPI::BlockVector &src) const;
1100 * const SmartPointer<const TrilinosWrappers::BlockSparseMatrix>
1102 * const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
1103 * PreconditionerTypeMp>>
1105 * const PreconditionerTypeA &a_preconditioner;
1107 * mutable TrilinosWrappers::MPI::Vector tmp;
1112 * template <class PreconditionerTypeA, class PreconditionerTypeMp>
1113 * BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::
1114 * BlockSchurPreconditioner(
1115 * const TrilinosWrappers::BlockSparseMatrix &S,
1116 * const InverseMatrix<TrilinosWrappers::SparseMatrix,
1117 * PreconditionerTypeMp> &Mpinv,
1118 * const PreconditionerTypeA & Apreconditioner)
1119 * : darcy_matrix(&S)
1120 * , m_inverse(&Mpinv)
1121 * , a_preconditioner(Apreconditioner)
1122 * , tmp(complete_index_set(darcy_matrix->block(1, 1).m()))
1126 * template <class PreconditionerTypeA, class PreconditionerTypeMp>
1128 * BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::vmult(
1129 * TrilinosWrappers::MPI::BlockVector & dst,
1130 * const TrilinosWrappers::MPI::BlockVector &src) const
1132 * a_preconditioner.vmult(dst.block(0), src.block(0));
1133 * darcy_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
1135 * m_inverse->vmult(dst.block(1), tmp);
1137 * } // namespace LinearSolvers
1143 * <a name="TheTwoPhaseFlowProblemclass"></a>
1144 * <h3>The TwoPhaseFlowProblem class</h3>
1148 * The definition of the class that defines the top-level logic of solving
1149 * the time-dependent advection-dominated two-phase flow problem (or
1150 * Buckley-Leverett problem [Buckley 1942]) is mainly based on tutorial
1151 * programs @ref step_21 "step-21" and @ref step_33 "step-33", and in particular on @ref step_31 "step-31" where we have
1152 * used basically the same general structure as done here. As in @ref step_31 "step-31",
1153 * the key routines to look for in the implementation below are the
1154 * <code>run()</code> and <code>solve()</code> functions.
1158 * The main difference to @ref step_31 "step-31" is that, since adaptive operator splitting
1159 * is considered, we need a couple more member variables to hold the last
1160 * two computed Darcy (velocity/pressure) solutions in addition to the
1161 * current one (which is either computed directly, or extrapolated from the
1162 * previous two), and we need to remember the last two times we computed the
1163 * Darcy solution. We also need a helper function that figures out whether
1164 * we do indeed need to recompute the Darcy solution.
1168 * Unlike @ref step_31 "step-31", this step uses one more AffineConstraints object called
1169 * darcy_preconditioner_constraints. This constraint object is used only for
1170 * assembling the matrix for the Darcy preconditioner and includes hanging
1171 * node constraints as well as Dirichlet boundary value constraints for the
1172 * pressure variable. We need this because we are building a Laplace matrix
1173 * for the pressure as an approximation of the Schur complement) which is
1174 * only positive definite if boundary conditions are applied.
1178 * The collection of member functions and variables thus declared in this
1179 * class is then rather similar to those in @ref step_31 "step-31":
1182 * template <int dim>
1183 * class TwoPhaseFlowProblem
1186 * TwoPhaseFlowProblem(const unsigned int degree);
1190 * void setup_dofs();
1191 * void assemble_darcy_preconditioner();
1192 * void build_darcy_preconditioner();
1193 * void assemble_darcy_system();
1194 * void assemble_saturation_system();
1195 * void assemble_saturation_matrix();
1196 * void assemble_saturation_rhs();
1197 * void assemble_saturation_rhs_cell_term(
1198 * const FEValues<dim> & saturation_fe_values,
1199 * const FEValues<dim> & darcy_fe_values,
1200 * const double global_max_u_F_prime,
1201 * const double global_S_variation,
1202 * const std::vector<types::global_dof_index> &local_dof_indices);
1203 * void assemble_saturation_rhs_boundary_term(
1204 * const FEFaceValues<dim> & saturation_fe_face_values,
1205 * const FEFaceValues<dim> & darcy_fe_face_values,
1206 * const std::vector<types::global_dof_index> &local_dof_indices);
1208 * void refine_mesh(const unsigned int min_grid_level,
1209 * const unsigned int max_grid_level);
1210 * void output_results() const;
1214 * We follow with a number of helper functions that are used in a variety
1215 * of places throughout the program:
1218 * double get_max_u_F_prime() const;
1219 * std::pair<double, double> get_extrapolated_saturation_range() const;
1220 * bool determine_whether_to_solve_for_pressure_and_velocity() const;
1221 * void project_back_saturation();
1222 * double compute_viscosity(
1223 * const std::vector<double> & old_saturation,
1224 * const std::vector<double> & old_old_saturation,
1225 * const std::vector<Tensor<1, dim>> &old_saturation_grads,
1226 * const std::vector<Tensor<1, dim>> &old_old_saturation_grads,
1227 * const std::vector<Vector<double>> &present_darcy_values,
1228 * const double global_max_u_F_prime,
1229 * const double global_S_variation,
1230 * const double cell_diameter) const;
1235 * This all is followed by the member variables, most of which are similar
1236 * to the ones in @ref step_31 "step-31", with the exception of the ones that pertain to
1237 * the macro time stepping for the velocity/pressure system:
1240 * Triangulation<dim> triangulation;
1241 * double global_Omega_diameter;
1243 * const unsigned int degree;
1245 * const unsigned int darcy_degree;
1246 * FESystem<dim> darcy_fe;
1247 * DoFHandler<dim> darcy_dof_handler;
1248 * AffineConstraints<double> darcy_constraints;
1250 * AffineConstraints<double> darcy_preconditioner_constraints;
1252 * TrilinosWrappers::BlockSparseMatrix darcy_matrix;
1253 * TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix;
1255 * TrilinosWrappers::MPI::BlockVector darcy_solution;
1256 * TrilinosWrappers::MPI::BlockVector darcy_rhs;
1258 * TrilinosWrappers::MPI::BlockVector last_computed_darcy_solution;
1259 * TrilinosWrappers::MPI::BlockVector second_last_computed_darcy_solution;
1262 * const unsigned int saturation_degree;
1263 * FE_Q<dim> saturation_fe;
1264 * DoFHandler<dim> saturation_dof_handler;
1265 * AffineConstraints<double> saturation_constraints;
1267 * TrilinosWrappers::SparseMatrix saturation_matrix;
1270 * TrilinosWrappers::MPI::Vector saturation_solution;
1271 * TrilinosWrappers::MPI::Vector old_saturation_solution;
1272 * TrilinosWrappers::MPI::Vector old_old_saturation_solution;
1273 * TrilinosWrappers::MPI::Vector saturation_rhs;
1275 * TrilinosWrappers::MPI::Vector
1276 * saturation_matching_last_computed_darcy_solution;
1278 * const double saturation_refinement_threshold;
1281 * const double end_time;
1283 * double current_macro_time_step;
1284 * double old_macro_time_step;
1287 * double old_time_step;
1288 * unsigned int timestep_number;
1290 * const double viscosity;
1291 * const double porosity;
1292 * const double AOS_threshold;
1294 * std::shared_ptr<TrilinosWrappers::PreconditionIC> Amg_preconditioner;
1295 * std::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner;
1297 * bool rebuild_saturation_matrix;
1301 * At the very end we declare a variable that denotes the material
1302 * model. Compared to @ref step_21 "step-21", we do this here as a member variable since
1303 * we will want to use it in a variety of places and so having a central
1304 * place where such a variable is declared will make it simpler to replace
1305 * one class by another (e.g. replace RandomMedium::KInverse by
1306 * SingleCurvingCrack::KInverse).
1309 * const RandomMedium::KInverse<dim> k_inverse;
1316 * <a name="TwoPhaseFlowProblemdimTwoPhaseFlowProblem"></a>
1317 * <h3>TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem</h3>
1321 * The constructor of this class is an extension of the constructors in
1322 * @ref step_21 "step-21" and @ref step_31 "step-31". We need to add the various variables that concern
1323 * the saturation. As discussed in the introduction, we are going to use
1324 * @f$Q_2 \times Q_1@f$ (Taylor-Hood) elements again for the Darcy system, an
1325 * element combination that fulfills the Ladyzhenskaya-Babuska-Brezzi (LBB)
1326 * conditions [Brezzi and Fortin 1991, Chen 2005], and @f$Q_1@f$ elements for
1327 * the saturation. However, by using variables that store the polynomial
1328 * degree of the Darcy and temperature finite elements, it is easy to
1329 * consistently modify the degree of the elements as well as all quadrature
1330 * formulas used on them downstream. Moreover, we initialize the time
1331 * stepping variables related to operator splitting as well as the option
1332 * for matrix assembly and preconditioning:
1335 * template <int dim>
1336 * TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem(const unsigned int degree)
1337 * : triangulation(Triangulation<dim>::maximum_smoothing)
1338 * , global_Omega_diameter(std::numeric_limits<double>::quiet_NaN())
1340 * , darcy_degree(degree)
1341 * , darcy_fe(FE_Q<dim>(darcy_degree + 1), dim, FE_Q<dim>(darcy_degree), 1)
1342 * , darcy_dof_handler(triangulation)
1345 * saturation_degree(degree + 1)
1346 * , saturation_fe(saturation_degree)
1347 * , saturation_dof_handler(triangulation)
1350 * saturation_refinement_threshold(0.5)
1357 * current_macro_time_step(0)
1358 * , old_macro_time_step(0)
1362 * , old_time_step(0)
1363 * , timestep_number(0)
1366 * , AOS_threshold(3.0)
1369 * rebuild_saturation_matrix(true)
1376 * <a name="TwoPhaseFlowProblemdimsetup_dofs"></a>
1377 * <h3>TwoPhaseFlowProblem<dim>::setup_dofs</h3>
1381 * This is the function that sets up the DoFHandler objects we have here
1382 * (one for the Darcy part and one for the saturation part) as well as set
1383 * to the right sizes the various objects required for the linear algebra in
1384 * this program. Its basic operations are similar to what @ref step_31 "step-31" did.
1388 * The body of the function first enumerates all degrees of freedom for the
1389 * Darcy and saturation systems. For the Darcy part, degrees of freedom are
1390 * then sorted to ensure that velocities precede pressure DoFs so that we
1391 * can partition the Darcy matrix into a @f$2 \times 2@f$ matrix.
1395 * Then, we need to incorporate hanging node constraints and Dirichlet
1396 * boundary value constraints into darcy_preconditioner_constraints. The
1397 * boundary condition constraints are only set on the pressure component
1398 * since the Schur complement preconditioner that corresponds to the porous
1399 * media flow operator in non-mixed form, @f$-\nabla \cdot [\mathbf K
1400 * \lambda_t(S)]\nabla@f$, acts only on the pressure variable. Therefore, we
1401 * use a component_mask that filters out the velocity component, so that the
1402 * condensation is performed on pressure degrees of freedom only.
1406 * After having done so, we count the number of degrees of freedom in the
1407 * various blocks. This information is then used to create the sparsity
1408 * pattern for the Darcy and saturation system matrices as well as the
1409 * preconditioner matrix from which we build the Darcy preconditioner. As in
1410 * @ref step_31 "step-31", we choose to create the pattern using the blocked version of
1411 * DynamicSparsityPattern. So, for this, we follow the same way as @ref step_31 "step-31"
1412 * did and we don't have to repeat descriptions again
for the rest of the
1416 * template <int dim>
1417 *
void TwoPhaseFlowProblem<dim>::setup_dofs()
1419 * std::vector<unsigned int> darcy_block_component(dim + 1, 0);
1420 * darcy_block_component[dim] = 1;
1422 * darcy_dof_handler.distribute_dofs(darcy_fe);
1426 * darcy_constraints.clear();
1428 * darcy_constraints);
1429 * darcy_constraints.close();
1432 * saturation_dof_handler.distribute_dofs(saturation_fe);
1434 * saturation_constraints.clear();
1436 * saturation_constraints);
1437 * saturation_constraints.close();
1440 * darcy_preconditioner_constraints.clear();
1445 * darcy_preconditioner_constraints);
1447 * darcy_preconditioner_constraints,
1448 * darcy_fe.component_mask(
1451 * darcy_preconditioner_constraints.close();
1455 *
const std::vector<types::global_dof_index> darcy_dofs_per_block =
1457 * darcy_block_component);
1458 *
const unsigned int n_u = darcy_dofs_per_block[0],
1459 * n_p = darcy_dofs_per_block[1],
1460 * n_s = saturation_dof_handler.n_dofs();
1462 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
1463 * <<
" (on " <<
triangulation.n_levels() <<
" levels)" << std::endl
1464 * <<
"Number of degrees of freedom: " << n_u + n_p + n_s <<
" ("
1465 * << n_u <<
'+' << n_p <<
'+' << n_s <<
')' << std::endl
1469 * darcy_matrix.clear();
1473 * dsp.block(0, 0).reinit(n_u, n_u);
1474 * dsp.block(0, 1).reinit(n_u, n_p);
1475 * dsp.block(1, 0).reinit(n_p, n_u);
1476 * dsp.block(1, 1).reinit(n_p, n_p);
1478 * dsp.collect_sizes();
1482 *
for (
unsigned int c = 0; c < dim + 1; ++c)
1483 *
for (
unsigned int d = 0;
d < dim + 1; ++
d)
1484 *
if (!((c == dim) && (
d == dim)))
1491 * darcy_dof_handler, coupling, dsp, darcy_constraints,
false);
1493 * darcy_matrix.reinit(dsp);
1497 * Amg_preconditioner.reset();
1498 * Mp_preconditioner.reset();
1499 * darcy_preconditioner_matrix.clear();
1503 * dsp.block(0, 0).reinit(n_u, n_u);
1504 * dsp.block(0, 1).reinit(n_u, n_p);
1505 * dsp.block(1, 0).reinit(n_p, n_u);
1506 * dsp.block(1, 1).reinit(n_p, n_p);
1508 * dsp.collect_sizes();
1511 *
for (
unsigned int c = 0; c < dim + 1; ++c)
1512 *
for (
unsigned int d = 0;
d < dim + 1; ++
d)
1519 * darcy_dof_handler, coupling, dsp, darcy_constraints,
false);
1521 * darcy_preconditioner_matrix.reinit(dsp);
1526 * saturation_matrix.clear();
1532 * saturation_constraints,
1536 * saturation_matrix.reinit(dsp);
1539 * std::vector<IndexSet> darcy_partitioning(2);
1542 * darcy_solution.reinit(darcy_partitioning, MPI_COMM_WORLD);
1543 * darcy_solution.collect_sizes();
1545 * last_computed_darcy_solution.reinit(darcy_partitioning, MPI_COMM_WORLD);
1546 * last_computed_darcy_solution.collect_sizes();
1548 * second_last_computed_darcy_solution.reinit(darcy_partitioning,
1550 * second_last_computed_darcy_solution.collect_sizes();
1552 * darcy_rhs.reinit(darcy_partitioning, MPI_COMM_WORLD);
1553 * darcy_rhs.collect_sizes();
1556 * saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1557 * old_saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1558 * old_old_saturation_solution.reinit(saturation_partitioning, MPI_COMM_WORLD);
1560 * saturation_matching_last_computed_darcy_solution.reinit(
1561 * saturation_partitioning, MPI_COMM_WORLD);
1563 * saturation_rhs.reinit(saturation_partitioning, MPI_COMM_WORLD);
1570 * <a name=
"Assemblingmatricesandpreconditioners"></a>
1571 * <h3>Assembling matrices and preconditioners</h3>
1575 * The next few
functions are devoted to setting up the various system and
1576 * preconditioner matrices and right hand sides that we have to deal with in
1582 * <a name=
"TwoPhaseFlowProblemdimassemble_darcy_preconditioner"></a>
1583 * <h4>TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner</h4>
1587 * This
function assembles the
matrix we use
for preconditioning the Darcy
1588 * system. What we need are a vector mass
matrix weighted by
1589 * @f$\left(\mathbf{K} \lambda_t\right)^{-1}@f$ on the velocity components and a
1590 * mass
matrix weighted by @f$\left(\mathbf{K} \lambda_t\right)@f$ on the
1591 * pressure component. We start by generating a quadrature
object of
1592 * appropriate order, the
FEValues object that can give values and gradients
1593 * at the quadrature points (together with quadrature weights). Next we
1594 * create data structures
for the cell
matrix and the relation between local
1595 * and global DoFs. The vectors phi_u and grad_phi_p are going to hold the
1596 * values of the basis
functions in order to faster build up the local
1597 * matrices, as was already done in @ref step_22
"step-22". Before we start the
loop over
1598 * all active cells, we have to specify which components are pressure and
1599 * which are velocity.
1603 * The creation of the local
matrix is rather simple. There are only a term
1604 * weighted by @f$\left(\mathbf{K} \lambda_t\right)^{-1}@f$ (on the velocity)
1605 * and a Laplace
matrix weighted by @f$\left(\mathbf{K} \lambda_t\right)@f$ to
1606 * be generated, so the creation of the local
matrix is done in essentially
1607 * two lines. Since the material model
functions at the top of
this file
1608 * only provide the inverses of the permeability and mobility, we have to
1609 * compute @f$\mathbf K@f$ and @f$\lambda_t@f$ by hand from the given values, once
1610 * per quadrature
point.
1614 * Once the local
matrix is ready (
loop over rows and columns in the local
1615 *
matrix on each quadrature
point), we get the local DoF indices and write
1616 * the local information into the global
matrix. We
do this by directly
1617 * applying the constraints (i.e. darcy_preconditioner_constraints) that
1618 * takes care of hanging node and
zero Dirichlet boundary condition
1619 * constraints. By doing so, we don
't have to do that afterwards, and we
1622 * modify
matrix and vector entries and so are difficult to write
for the
1623 * Trilinos classes where we don
't immediately have access to individual
1627 * template <int dim>
1628 * void TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner()
1630 * std::cout << " Rebuilding darcy preconditioner..." << std::endl;
1632 * darcy_preconditioner_matrix = 0;
1634 * const QGauss<dim> quadrature_formula(darcy_degree + 2);
1635 * FEValues<dim> darcy_fe_values(darcy_fe,
1636 * quadrature_formula,
1637 * update_JxW_values | update_values |
1638 * update_gradients |
1639 * update_quadrature_points);
1640 * FEValues<dim> saturation_fe_values(saturation_fe,
1641 * quadrature_formula,
1644 * const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell;
1645 * const unsigned int n_q_points = quadrature_formula.size();
1647 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1649 * std::vector<double> old_saturation_values(n_q_points);
1651 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1652 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1654 * std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
1655 * std::vector<Tensor<1, dim>> grad_phi_p(dofs_per_cell);
1657 * const FEValuesExtractors::Vector velocities(0);
1658 * const FEValuesExtractors::Scalar pressure(dim);
1660 * auto cell = darcy_dof_handler.begin_active();
1661 * const auto endc = darcy_dof_handler.end();
1662 * auto saturation_cell = saturation_dof_handler.begin_active();
1664 * for (; cell != endc; ++cell, ++saturation_cell)
1666 * darcy_fe_values.reinit(cell);
1667 * saturation_fe_values.reinit(saturation_cell);
1671 * saturation_fe_values.get_function_values(old_saturation_solution,
1672 * old_saturation_values);
1674 * k_inverse.value_list(darcy_fe_values.get_quadrature_points(),
1675 * k_inverse_values);
1677 * for (unsigned int q = 0; q < n_q_points; ++q)
1679 * const double old_s = old_saturation_values[q];
1681 * const double inverse_mobility = mobility_inverse(old_s, viscosity);
1682 * const double mobility = 1.0 / inverse_mobility;
1683 * const Tensor<2, dim> permeability = invert(k_inverse_values[q]);
1685 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
1687 * phi_u[k] = darcy_fe_values[velocities].value(k, q);
1688 * grad_phi_p[k] = darcy_fe_values[pressure].gradient(k, q);
1691 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1692 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1694 * local_matrix(i, j) +=
1695 * (k_inverse_values[q] * inverse_mobility * phi_u[i] *
1697 * permeability * mobility * grad_phi_p[i] * grad_phi_p[j]) *
1698 * darcy_fe_values.JxW(q);
1702 * cell->get_dof_indices(local_dof_indices);
1703 * darcy_preconditioner_constraints.distribute_local_to_global(
1704 * local_matrix, local_dof_indices, darcy_preconditioner_matrix);
1712 * <a name="TwoPhaseFlowProblemdimbuild_darcy_preconditioner"></a>
1713 * <h4>TwoPhaseFlowProblem<dim>::build_darcy_preconditioner</h4>
1717 * After calling the above functions to assemble the preconditioner matrix,
1718 * this function generates the inner preconditioners that are going to be
1719 * used for the Schur complement block preconditioner. The preconditioners
1720 * need to be regenerated at every saturation time step since they depend on
1721 * the saturation @f$S@f$ that varies with time.
1725 * In here, we set up the preconditioner for the velocity-velocity matrix
1726 * @f$\mathbf{M}^{\mathbf{u}}@f$ and the Schur complement @f$\mathbf{S}@f$. As
1727 * explained in the introduction, we are going to use an IC preconditioner
1728 * based on the vector matrix @f$\mathbf{M}^{\mathbf{u}}@f$ and another based on
1729 * the scalar Laplace matrix @f$\tilde{\mathbf{S}}^p@f$ (which is spectrally
1730 * close to the Schur complement of the Darcy matrix). Usually, the
1731 * TrilinosWrappers::PreconditionIC class can be seen as a good black-box
1732 * preconditioner which does not need any special knowledge of the matrix
1733 * structure and/or the operator that's behind it.
1736 *
template <
int dim>
1737 *
void TwoPhaseFlowProblem<dim>::build_darcy_preconditioner()
1739 * assemble_darcy_preconditioner();
1741 * Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionIC>();
1742 * Amg_preconditioner->initialize(darcy_preconditioner_matrix.block(0, 0));
1744 * Mp_preconditioner = std::make_shared<TrilinosWrappers::PreconditionIC>();
1745 * Mp_preconditioner->initialize(darcy_preconditioner_matrix.block(1, 1));
1752 * <a name=
"TwoPhaseFlowProblemdimassemble_darcy_system"></a>
1753 * <h4>TwoPhaseFlowProblem<dim>::assemble_darcy_system</h4>
1757 * This is the
function that assembles the linear system
for the Darcy
1762 * Regarding the technical details of implementation, the procedures are
1763 * similar to those in @ref step_22
"step-22" and @ref step_31
"step-31". We reset
matrix and vector,
1764 * create a quadrature formula on the cells, and then create the respective
1769 * There is
one thing that needs to be commented: since we have a separate
1770 * finite element and
DoFHandler for the saturation, we need to generate a
1771 *
second FEValues object for the proper evaluation of the saturation
1772 * solution. This isn
't too complicated to realize here: just use the
1773 * saturation structures and set an update flag for the basis function
1774 * values which we need for evaluation of the saturation solution. The only
1775 * important part to remember here is that the same quadrature formula is
1776 * used for both FEValues objects to ensure that we get matching information
1777 * when we loop over the quadrature points of the two objects.
1781 * The declarations proceed with some shortcuts for array sizes, the
1782 * creation of the local matrix, right hand side as well as the vector for
1783 * the indices of the local dofs compared to the global system.
1786 * template <int dim>
1787 * void TwoPhaseFlowProblem<dim>::assemble_darcy_system()
1792 * QGauss<dim> quadrature_formula(darcy_degree + 2);
1793 * QGauss<dim - 1> face_quadrature_formula(darcy_degree + 2);
1795 * FEValues<dim> darcy_fe_values(darcy_fe,
1796 * quadrature_formula,
1797 * update_values | update_gradients |
1798 * update_quadrature_points |
1799 * update_JxW_values);
1801 * FEValues<dim> saturation_fe_values(saturation_fe,
1802 * quadrature_formula,
1805 * FEFaceValues<dim> darcy_fe_face_values(darcy_fe,
1806 * face_quadrature_formula,
1808 * update_normal_vectors |
1809 * update_quadrature_points |
1810 * update_JxW_values);
1812 * const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell;
1814 * const unsigned int n_q_points = quadrature_formula.size();
1815 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1817 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1818 * Vector<double> local_rhs(dofs_per_cell);
1820 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1822 * const Functions::ZeroFunction<dim> pressure_right_hand_side;
1823 * const PressureBoundaryValues<dim> pressure_boundary_values;
1825 * std::vector<double> pressure_rhs_values(n_q_points);
1826 * std::vector<double> boundary_values(n_face_q_points);
1827 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
1831 * Next we need a vector that will contain the values of the saturation
1832 * solution at the previous time level at the quadrature points to
1833 * assemble the saturation dependent coefficients in the Darcy equations.
1837 * The set of vectors we create next hold the evaluations of the basis
1838 * functions as well as their gradients that will be used for creating the
1839 * matrices. Putting these into their own arrays rather than asking the
1840 * FEValues object for this information each time it is needed is an
1841 * optimization to accelerate the assembly process, see @ref step_22 "step-22" for
1846 * The last two declarations are used to extract the individual blocks
1847 * (velocity, pressure, saturation) from the total FE system.
1850 * std::vector<double> old_saturation_values(n_q_points);
1852 * std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
1853 * std::vector<double> div_phi_u(dofs_per_cell);
1854 * std::vector<double> phi_p(dofs_per_cell);
1856 * const FEValuesExtractors::Vector velocities(0);
1857 * const FEValuesExtractors::Scalar pressure(dim);
1861 * Now start the loop over all cells in the problem. We are working on two
1862 * different DoFHandlers for this assembly routine, so we must have two
1863 * different cell iterators for the two objects in use. This might seem a
1864 * bit peculiar, but since both the Darcy system and the saturation system
1865 * use the same grid we can assume that the two iterators run in sync over
1866 * the cells of the two DoFHandler objects.
1870 * The first statements within the loop are again all very familiar, doing
1871 * the update of the finite element data as specified by the update flags,
1872 * zeroing out the local arrays and getting the values of the old solution
1873 * at the quadrature points. At this point we also have to get the values
1874 * of the saturation function of the previous time step at the quadrature
1875 * points. To this end, we can use the FEValues::get_function_values
1876 * (previously already used in @ref step_9 "step-9", @ref step_14 "step-14" and @ref step_15 "step-15"), a function
1877 * that takes a solution vector and returns a list of function values at
1878 * the quadrature points of the present cell. In fact, it returns the
1879 * complete vector-valued solution at each quadrature point, i.e. not only
1880 * the saturation but also the velocities and pressure.
1884 * Then we are ready to loop over the quadrature points on the cell to do
1885 * the integration. The formula for this follows in a straightforward way
1886 * from what has been discussed in the introduction.
1890 * Once this is done, we start the loop over the rows and columns of the
1891 * local matrix and feed the matrix with the relevant products.
1895 * The last step in the loop over all cells is to enter the local
1896 * contributions into the global matrix and vector structures to the
1897 * positions specified in local_dof_indices. Again, we let the
1898 * AffineConstraints class do the insertion of the cell matrix
1899 * elements to the global matrix, which already condenses the hanging node
1903 * auto cell = darcy_dof_handler.begin_active();
1904 * const auto endc = darcy_dof_handler.end();
1905 * auto saturation_cell = saturation_dof_handler.begin_active();
1907 * for (; cell != endc; ++cell, ++saturation_cell)
1909 * darcy_fe_values.reinit(cell);
1910 * saturation_fe_values.reinit(saturation_cell);
1915 * saturation_fe_values.get_function_values(old_saturation_solution,
1916 * old_saturation_values);
1918 * pressure_right_hand_side.value_list(
1919 * darcy_fe_values.get_quadrature_points(), pressure_rhs_values);
1920 * k_inverse.value_list(darcy_fe_values.get_quadrature_points(),
1921 * k_inverse_values);
1923 * for (unsigned int q = 0; q < n_q_points; ++q)
1925 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
1927 * phi_u[k] = darcy_fe_values[velocities].value(k, q);
1928 * div_phi_u[k] = darcy_fe_values[velocities].divergence(k, q);
1929 * phi_p[k] = darcy_fe_values[pressure].value(k, q);
1931 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1933 * const double old_s = old_saturation_values[q];
1934 * for (unsigned int j = 0; j <= i; ++j)
1936 * local_matrix(i, j) +=
1937 * (phi_u[i] * k_inverse_values[q] *
1938 * mobility_inverse(old_s, viscosity) * phi_u[j] -
1939 * div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
1940 * darcy_fe_values.JxW(q);
1944 * (-phi_p[i] * pressure_rhs_values[q]) * darcy_fe_values.JxW(q);
1948 * for (const auto &face : cell->face_iterators())
1949 * if (face->at_boundary())
1951 * darcy_fe_face_values.reinit(cell, face);
1953 * pressure_boundary_values.value_list(
1954 * darcy_fe_face_values.get_quadrature_points(), boundary_values);
1956 * for (unsigned int q = 0; q < n_face_q_points; ++q)
1957 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1959 * const Tensor<1, dim> phi_i_u =
1960 * darcy_fe_face_values[velocities].value(i, q);
1963 * -(phi_i_u * darcy_fe_face_values.normal_vector(q) *
1964 * boundary_values[q] * darcy_fe_face_values.JxW(q));
1968 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1969 * for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
1970 * local_matrix(i, j) = local_matrix(j, i);
1972 * cell->get_dof_indices(local_dof_indices);
1974 * darcy_constraints.distribute_local_to_global(
1975 * local_matrix, local_rhs, local_dof_indices, darcy_matrix, darcy_rhs);
1983 * <a name="TwoPhaseFlowProblemdimassemble_saturation_system"></a>
1984 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_system</h4>
1988 * This function is to assemble the linear system for the saturation
1989 * transport equation. It calls, if necessary, two other member functions:
1990 * assemble_saturation_matrix() and assemble_saturation_rhs(). The former
1991 * function then assembles the saturation matrix that only needs to be
1992 * changed occasionally. On the other hand, the latter function that
1993 * assembles the right hand side must be called at every saturation time
1997 * template <int dim>
1998 * void TwoPhaseFlowProblem<dim>::assemble_saturation_system()
2000 * if (rebuild_saturation_matrix == true)
2002 * saturation_matrix = 0;
2003 * assemble_saturation_matrix();
2006 * saturation_rhs = 0;
2007 * assemble_saturation_rhs();
2015 * <a name="TwoPhaseFlowProblemdimassemble_saturation_matrix"></a>
2016 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_matrix</h4>
2020 * This function is easily understood since it only forms a simple mass
2021 * matrix for the left hand side of the saturation linear system by basis
2022 * functions phi_i_s and phi_j_s only. Finally, as usual, we enter the local
2023 * contribution into the global matrix by specifying the position in
2024 * local_dof_indices. This is done by letting the AffineConstraints class do
2025 * the insertion of the cell matrix elements to the global matrix, which
2026 * already condenses the hanging node constraints.
2029 * template <int dim>
2030 * void TwoPhaseFlowProblem<dim>::assemble_saturation_matrix()
2032 * QGauss<dim> quadrature_formula(saturation_degree + 2);
2034 * FEValues<dim> saturation_fe_values(saturation_fe,
2035 * quadrature_formula,
2036 * update_values | update_JxW_values);
2038 * const unsigned int dofs_per_cell = saturation_fe.dofs_per_cell;
2040 * const unsigned int n_q_points = quadrature_formula.size();
2042 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2043 * Vector<double> local_rhs(dofs_per_cell);
2045 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2047 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2049 * saturation_fe_values.reinit(cell);
2053 * for (unsigned int q = 0; q < n_q_points; ++q)
2054 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2056 * const double phi_i_s = saturation_fe_values.shape_value(i, q);
2057 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2059 * const double phi_j_s = saturation_fe_values.shape_value(j, q);
2060 * local_matrix(i, j) +=
2061 * porosity * phi_i_s * phi_j_s * saturation_fe_values.JxW(q);
2064 * cell->get_dof_indices(local_dof_indices);
2066 * saturation_constraints.distribute_local_to_global(local_matrix,
2067 * local_dof_indices,
2068 * saturation_matrix);
2077 * <a name="TwoPhaseFlowProblemdimassemble_saturation_rhs"></a>
2078 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs</h4>
2082 * This function is to assemble the right hand side of the saturation
2083 * transport equation. Before going about it, we have to create two FEValues
2084 * objects for the Darcy and saturation systems respectively and, in
2085 * addition, two FEFaceValues objects for the two systems because we have a
2086 * boundary integral term in the weak form of saturation equation. For the
2087 * FEFaceValues object of the saturation system, we also require normal
2088 * vectors, which we request using the update_normal_vectors flag.
2092 * Next, before looping over all the cells, we have to compute some
2093 * parameters (e.g. global_u_infty, global_S_variation, and
2094 * global_Omega_diameter) that the artificial viscosity @f$\nu@f$ needs. This is
2095 * largely the same as was done in @ref step_31 "step-31", so you may see there for more
2100 * The real works starts with the loop over all the saturation and Darcy
2101 * cells to put the local contributions into the global vector. In this
2102 * loop, in order to simplify the implementation, we split some of the work
2103 * into two helper functions: assemble_saturation_rhs_cell_term and
2104 * assemble_saturation_rhs_boundary_term. We note that we insert cell or
2105 * boundary contributions into the global vector in the two functions rather
2106 * than in this present function.
2109 * template <int dim>
2110 * void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs()
2112 * QGauss<dim> quadrature_formula(saturation_degree + 2);
2113 * QGauss<dim - 1> face_quadrature_formula(saturation_degree + 2);
2115 * FEValues<dim> saturation_fe_values(saturation_fe,
2116 * quadrature_formula,
2117 * update_values | update_gradients |
2118 * update_quadrature_points |
2119 * update_JxW_values);
2120 * FEValues<dim> darcy_fe_values(darcy_fe, quadrature_formula, update_values);
2121 * FEFaceValues<dim> saturation_fe_face_values(saturation_fe,
2122 * face_quadrature_formula,
2124 * update_normal_vectors |
2125 * update_quadrature_points |
2126 * update_JxW_values);
2127 * FEFaceValues<dim> darcy_fe_face_values(darcy_fe,
2128 * face_quadrature_formula,
2130 * FEFaceValues<dim> saturation_fe_face_values_neighbor(
2131 * saturation_fe, face_quadrature_formula, update_values);
2133 * const unsigned int dofs_per_cell =
2134 * saturation_dof_handler.get_fe().dofs_per_cell;
2135 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2137 * const double global_max_u_F_prime = get_max_u_F_prime();
2138 * const std::pair<double, double> global_S_range =
2139 * get_extrapolated_saturation_range();
2140 * const double global_S_variation =
2141 * global_S_range.second - global_S_range.first;
2143 * auto cell = saturation_dof_handler.begin_active();
2144 * const auto endc = saturation_dof_handler.end();
2145 * auto darcy_cell = darcy_dof_handler.begin_active();
2146 * for (; cell != endc; ++cell, ++darcy_cell)
2148 * saturation_fe_values.reinit(cell);
2149 * darcy_fe_values.reinit(darcy_cell);
2151 * cell->get_dof_indices(local_dof_indices);
2153 * assemble_saturation_rhs_cell_term(saturation_fe_values,
2155 * global_max_u_F_prime,
2156 * global_S_variation,
2157 * local_dof_indices);
2159 * for (const auto &face : cell->face_iterators())
2160 * if (face->at_boundary())
2162 * darcy_fe_face_values.reinit(darcy_cell, face);
2163 * saturation_fe_face_values.reinit(cell, face);
2164 * assemble_saturation_rhs_boundary_term(saturation_fe_face_values,
2165 * darcy_fe_face_values,
2166 * local_dof_indices);
2176 * <a name="TwoPhaseFlowProblemdimassemble_saturation_rhs_cell_term"></a>
2177 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_cell_term</h4>
2181 * This function takes care of integrating the cell terms of the right hand
2182 * side of the saturation equation, and then assembling it into the global
2183 * right hand side vector. Given the discussion in the introduction, the
2184 * form of these contributions is clear. The only tricky part is getting the
2185 * artificial viscosity and all that is necessary to compute it. The first
2186 * half of the function is devoted to this task.
2190 * The last part of the function is copying the local contributions into the
2191 * global vector with position specified in local_dof_indices.
2194 * template <int dim>
2195 * void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_cell_term(
2196 * const FEValues<dim> & saturation_fe_values,
2197 * const FEValues<dim> & darcy_fe_values,
2198 * const double global_max_u_F_prime,
2199 * const double global_S_variation,
2200 * const std::vector<types::global_dof_index> &local_dof_indices)
2202 * const unsigned int dofs_per_cell = saturation_fe_values.dofs_per_cell;
2203 * const unsigned int n_q_points = saturation_fe_values.n_quadrature_points;
2205 * std::vector<double> old_saturation_solution_values(n_q_points);
2206 * std::vector<double> old_old_saturation_solution_values(n_q_points);
2207 * std::vector<Tensor<1, dim>> old_grad_saturation_solution_values(n_q_points);
2208 * std::vector<Tensor<1, dim>> old_old_grad_saturation_solution_values(
2210 * std::vector<Vector<double>> present_darcy_solution_values(
2211 * n_q_points, Vector<double>(dim + 1));
2213 * saturation_fe_values.get_function_values(old_saturation_solution,
2214 * old_saturation_solution_values);
2215 * saturation_fe_values.get_function_values(
2216 * old_old_saturation_solution, old_old_saturation_solution_values);
2217 * saturation_fe_values.get_function_gradients(
2218 * old_saturation_solution, old_grad_saturation_solution_values);
2219 * saturation_fe_values.get_function_gradients(
2220 * old_old_saturation_solution, old_old_grad_saturation_solution_values);
2221 * darcy_fe_values.get_function_values(darcy_solution,
2222 * present_darcy_solution_values);
2225 * compute_viscosity(old_saturation_solution_values,
2226 * old_old_saturation_solution_values,
2227 * old_grad_saturation_solution_values,
2228 * old_old_grad_saturation_solution_values,
2229 * present_darcy_solution_values,
2230 * global_max_u_F_prime,
2231 * global_S_variation,
2232 * saturation_fe_values.get_cell()->diameter());
2234 * Vector<double> local_rhs(dofs_per_cell);
2236 * for (unsigned int q = 0; q < n_q_points; ++q)
2237 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2239 * const double old_s = old_saturation_solution_values[q];
2240 * Tensor<1, dim> present_u;
2241 * for (unsigned int d = 0; d < dim; ++d)
2242 * present_u[d] = present_darcy_solution_values[q](d);
2244 * const double phi_i_s = saturation_fe_values.shape_value(i, q);
2245 * const Tensor<1, dim> grad_phi_i_s =
2246 * saturation_fe_values.shape_grad(i, q);
2249 * (time_step * fractional_flow(old_s, viscosity) * present_u *
2251 * time_step * nu * old_grad_saturation_solution_values[q] *
2253 * porosity * old_s * phi_i_s) *
2254 * saturation_fe_values.JxW(q);
2257 * saturation_constraints.distribute_local_to_global(local_rhs,
2258 * local_dof_indices,
2266 * <a name="TwoPhaseFlowProblemdimassemble_saturation_rhs_boundary_term"></a>
2267 * <h4>TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_boundary_term</h4>
2271 * The next function is responsible for the boundary integral terms in the
2272 * right hand side form of the saturation equation. For these, we have to
2273 * compute the upwinding flux on the global boundary faces, i.e. we impose
2274 * Dirichlet boundary conditions weakly only on inflow parts of the global
2275 * boundary. As before, this has been described in @ref step_21 "step-21" so we refrain
2276 * from giving more descriptions about that.
2279 * template <int dim>
2280 * void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_boundary_term(
2281 * const FEFaceValues<dim> & saturation_fe_face_values,
2282 * const FEFaceValues<dim> & darcy_fe_face_values,
2283 * const std::vector<types::global_dof_index> &local_dof_indices)
2285 * const unsigned int dofs_per_cell = saturation_fe_face_values.dofs_per_cell;
2286 * const unsigned int n_face_q_points =
2287 * saturation_fe_face_values.n_quadrature_points;
2289 * Vector<double> local_rhs(dofs_per_cell);
2291 * std::vector<double> old_saturation_solution_values_face(n_face_q_points);
2292 * std::vector<Vector<double>> present_darcy_solution_values_face(
2293 * n_face_q_points, Vector<double>(dim + 1));
2294 * std::vector<double> neighbor_saturation(n_face_q_points);
2296 * saturation_fe_face_values.get_function_values(
2297 * old_saturation_solution, old_saturation_solution_values_face);
2298 * darcy_fe_face_values.get_function_values(
2299 * darcy_solution, present_darcy_solution_values_face);
2301 * SaturationBoundaryValues<dim> saturation_boundary_values;
2302 * saturation_boundary_values.value_list(
2303 * saturation_fe_face_values.get_quadrature_points(), neighbor_saturation);
2305 * for (unsigned int q = 0; q < n_face_q_points; ++q)
2307 * Tensor<1, dim> present_u_face;
2308 * for (unsigned int d = 0; d < dim; ++d)
2309 * present_u_face[d] = present_darcy_solution_values_face[q](d);
2311 * const double normal_flux =
2312 * present_u_face * saturation_fe_face_values.normal_vector(q);
2314 * const bool is_outflow_q_point = (normal_flux >= 0);
2316 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2318 * time_step * normal_flux *
2319 * fractional_flow((is_outflow_q_point == true ?
2320 * old_saturation_solution_values_face[q] :
2321 * neighbor_saturation[q]),
2323 * saturation_fe_face_values.shape_value(i, q) *
2324 * saturation_fe_face_values.JxW(q);
2326 * saturation_constraints.distribute_local_to_global(local_rhs,
2327 * local_dof_indices,
2335 * <a name="TwoPhaseFlowProblemdimsolve"></a>
2336 * <h3>TwoPhaseFlowProblem<dim>::solve</h3>
2340 * This function implements the operator splitting algorithm, i.e. in each
2341 * time step it either re-computes the solution of the Darcy system or
2342 * extrapolates velocity/pressure from previous time steps, then determines
2343 * the size of the time step, and then updates the saturation variable. The
2344 * implementation largely follows similar code in @ref step_31 "step-31". It is, next to
2345 * the run() function, the central one in this program.
2349 * At the beginning of the function, we ask whether to solve the
2350 * pressure-velocity part by evaluating the a posteriori criterion (see the
2351 * following function). If necessary, we will solve the pressure-velocity
2352 * part using the GMRES solver with the Schur complement block
2353 * preconditioner as is described in the introduction.
2356 * template <int dim>
2357 * void TwoPhaseFlowProblem<dim>::solve()
2359 * const bool solve_for_pressure_and_velocity =
2360 * determine_whether_to_solve_for_pressure_and_velocity();
2362 * if (solve_for_pressure_and_velocity == true)
2364 * std::cout << " Solving Darcy (pressure-velocity) system..."
2367 * assemble_darcy_system();
2368 * build_darcy_preconditioner();
2371 * const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
2372 * TrilinosWrappers::PreconditionIC>
2373 * mp_inverse(darcy_preconditioner_matrix.block(1, 1),
2374 * *Mp_preconditioner);
2376 * const LinearSolvers::BlockSchurPreconditioner<
2377 * TrilinosWrappers::PreconditionIC,
2378 * TrilinosWrappers::PreconditionIC>
2379 * preconditioner(darcy_matrix, mp_inverse, *Amg_preconditioner);
2381 * SolverControl solver_control(darcy_matrix.m(),
2382 * 1e-16 * darcy_rhs.l2_norm());
2384 * SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
2386 * SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(
2389 * for (unsigned int i = 0; i < darcy_solution.size(); ++i)
2390 * if (darcy_constraints.is_constrained(i))
2391 * darcy_solution(i) = 0;
2393 * gmres.solve(darcy_matrix, darcy_solution, darcy_rhs, preconditioner);
2395 * darcy_constraints.distribute(darcy_solution);
2397 * std::cout << " ..." << solver_control.last_step()
2398 * << " GMRES iterations." << std::endl;
2402 * second_last_computed_darcy_solution = last_computed_darcy_solution;
2403 * last_computed_darcy_solution = darcy_solution;
2405 * saturation_matching_last_computed_darcy_solution =
2406 * saturation_solution;
2411 * On the other hand, if we have decided that we don't want to compute the
2412 * solution of the Darcy system
for the current time step, then we need to
2413 * simply
extrapolate the previous two Darcy solutions to the same time as
2414 * we would have computed the velocity/pressure at. We
do a simple linear
2415 * extrapolation, i.e. given the current length @f$dt@f$ of the macro time
2416 * step from the time when we last computed the Darcy solution to now
2417 * (given by <code>current_macro_time_step</code>), and @f$DT@f$ the length of
2418 * the last macro time step (given by <code>old_macro_time_step</code>),
2419 * then we get @f$u^\ast = u_p + dt \frac{u_p-u_{pp}}{DT} = (1+dt/DT)u_p -
2420 * dt/DT u_{pp}@f$, where @f$u_p@f$ and @f$u_{pp}@f$ are the last two computed Darcy
2421 * solutions. We can implement
this formula
using just two lines of code.
2425 * Note that the algorithm here only works
if we have at least two
2426 * previously computed Darcy solutions from which we can
extrapolate to
2427 * the current time, and
this is ensured by requiring re-computation of
2428 * the Darcy solution
for the
first 2 time steps.
2433 * darcy_solution = last_computed_darcy_solution;
2434 * darcy_solution.sadd(1 + current_macro_time_step / old_macro_time_step,
2435 * -current_macro_time_step / old_macro_time_step,
2436 * second_last_computed_darcy_solution);
2442 * With the so computed velocity vector, compute the optimal time step
2443 * based on the CFL criterion discussed in the introduction...
2447 * old_time_step = time_step;
2449 *
const double max_u_F_prime = get_max_u_F_prime();
2450 *
if (max_u_F_prime > 0)
2452 * saturation_degree / max_u_F_prime / 50;
2454 * time_step = end_time - time;
2461 * ...and then also update the length of the macro time steps we use
while
2462 * we
're dealing with time step sizes. In particular, this involves: (i)
2463 * If we have just recomputed the Darcy solution, then the length of the
2464 * previous macro time step is now fixed and the length of the current
2465 * macro time step is, up to now, simply the length of the current (micro)
2466 * time step. (ii) If we have not recomputed the Darcy solution, then the
2467 * length of the current macro time step has just grown by
2468 * <code>time_step</code>.
2471 * if (solve_for_pressure_and_velocity == true)
2473 * old_macro_time_step = current_macro_time_step;
2474 * current_macro_time_step = time_step;
2477 * current_macro_time_step += time_step;
2481 * The last step in this function is to recompute the saturation solution
2482 * based on the velocity field we've just obtained. This naturally happens
2483 * in every time step, and we don
't skip any of these computations. At the
2484 * end of computing the saturation, we project back into the allowed
2485 * interval @f$[0,1]@f$ to make sure our solution remains physical.
2489 * std::cout << " Solving saturation transport equation..." << std::endl;
2491 * assemble_saturation_system();
2493 * SolverControl solver_control(saturation_matrix.m(),
2494 * 1e-16 * saturation_rhs.l2_norm());
2495 * SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
2497 * TrilinosWrappers::PreconditionIC preconditioner;
2498 * preconditioner.initialize(saturation_matrix);
2500 * cg.solve(saturation_matrix,
2501 * saturation_solution,
2505 * saturation_constraints.distribute(saturation_solution);
2506 * project_back_saturation();
2508 * std::cout << " ..." << solver_control.last_step()
2509 * << " CG iterations." << std::endl;
2517 * <a name="TwoPhaseFlowProblemdimrefine_mesh"></a>
2518 * <h3>TwoPhaseFlowProblem<dim>::refine_mesh</h3>
2522 * The next function does the refinement and coarsening of the mesh. It does
2523 * its work in three blocks: (i) Compute refinement indicators by looking at
2524 * the gradient of a solution vector extrapolated linearly from the previous
2525 * two using the respective sizes of the time step (or taking the only
2526 * solution we have if this is the first time step). (ii) Flagging those
2527 * cells for refinement and coarsening where the gradient is larger or
2528 * smaller than a certain threshold, preserving minimal and maximal levels
2529 * of mesh refinement. (iii) Transferring the solution from the old to the
2530 * new mesh. None of this is particularly difficult.
2533 * template <int dim>
2534 * void TwoPhaseFlowProblem<dim>::refine_mesh(const unsigned int min_grid_level,
2535 * const unsigned int max_grid_level)
2537 * Vector<double> refinement_indicators(triangulation.n_active_cells());
2539 * const QMidpoint<dim> quadrature_formula;
2540 * FEValues<dim> fe_values(saturation_fe,
2541 * quadrature_formula,
2542 * update_gradients);
2543 * std::vector<Tensor<1, dim>> grad_saturation(1);
2545 * TrilinosWrappers::MPI::Vector extrapolated_saturation_solution(
2546 * saturation_solution);
2547 * if (timestep_number != 0)
2548 * extrapolated_saturation_solution.sadd((1. + time_step / old_time_step),
2549 * time_step / old_time_step,
2550 * old_saturation_solution);
2552 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2554 * const unsigned int cell_no = cell->active_cell_index();
2555 * fe_values.reinit(cell);
2556 * fe_values.get_function_gradients(extrapolated_saturation_solution,
2559 * refinement_indicators(cell_no) = grad_saturation[0].norm();
2564 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2566 * const unsigned int cell_no = cell->active_cell_index();
2567 * cell->clear_coarsen_flag();
2568 * cell->clear_refine_flag();
2570 * if ((static_cast<unsigned int>(cell->level()) < max_grid_level) &&
2571 * (std::fabs(refinement_indicators(cell_no)) >
2572 * saturation_refinement_threshold))
2573 * cell->set_refine_flag();
2574 * else if ((static_cast<unsigned int>(cell->level()) >
2575 * min_grid_level) &&
2576 * (std::fabs(refinement_indicators(cell_no)) <
2577 * 0.5 * saturation_refinement_threshold))
2578 * cell->set_coarsen_flag();
2582 * triangulation.prepare_coarsening_and_refinement();
2585 * std::vector<TrilinosWrappers::MPI::Vector> x_saturation(3);
2586 * x_saturation[0] = saturation_solution;
2587 * x_saturation[1] = old_saturation_solution;
2588 * x_saturation[2] = saturation_matching_last_computed_darcy_solution;
2590 * std::vector<TrilinosWrappers::MPI::BlockVector> x_darcy(2);
2591 * x_darcy[0] = last_computed_darcy_solution;
2592 * x_darcy[1] = second_last_computed_darcy_solution;
2594 * SolutionTransfer<dim, TrilinosWrappers::MPI::Vector> saturation_soltrans(
2595 * saturation_dof_handler);
2597 * SolutionTransfer<dim, TrilinosWrappers::MPI::BlockVector> darcy_soltrans(
2598 * darcy_dof_handler);
2601 * triangulation.prepare_coarsening_and_refinement();
2602 * saturation_soltrans.prepare_for_coarsening_and_refinement(x_saturation);
2604 * darcy_soltrans.prepare_for_coarsening_and_refinement(x_darcy);
2606 * triangulation.execute_coarsening_and_refinement();
2609 * std::vector<TrilinosWrappers::MPI::Vector> tmp_saturation(3);
2610 * tmp_saturation[0].reinit(saturation_solution);
2611 * tmp_saturation[1].reinit(saturation_solution);
2612 * tmp_saturation[2].reinit(saturation_solution);
2613 * saturation_soltrans.interpolate(x_saturation, tmp_saturation);
2615 * saturation_solution = tmp_saturation[0];
2616 * old_saturation_solution = tmp_saturation[1];
2617 * saturation_matching_last_computed_darcy_solution = tmp_saturation[2];
2619 * saturation_constraints.distribute(saturation_solution);
2620 * saturation_constraints.distribute(old_saturation_solution);
2621 * saturation_constraints.distribute(
2622 * saturation_matching_last_computed_darcy_solution);
2624 * std::vector<TrilinosWrappers::MPI::BlockVector> tmp_darcy(2);
2625 * tmp_darcy[0].reinit(darcy_solution);
2626 * tmp_darcy[1].reinit(darcy_solution);
2627 * darcy_soltrans.interpolate(x_darcy, tmp_darcy);
2629 * last_computed_darcy_solution = tmp_darcy[0];
2630 * second_last_computed_darcy_solution = tmp_darcy[1];
2632 * darcy_constraints.distribute(last_computed_darcy_solution);
2633 * darcy_constraints.distribute(second_last_computed_darcy_solution);
2635 * rebuild_saturation_matrix = true;
2644 * <a name="TwoPhaseFlowProblemdimoutput_results"></a>
2645 * <h3>TwoPhaseFlowProblem<dim>::output_results</h3>
2649 * This function generates graphical output. It is in essence a copy of the
2650 * implementation in @ref step_31 "step-31".
2653 * template <int dim>
2654 * void TwoPhaseFlowProblem<dim>::output_results() const
2656 * const FESystem<dim> joint_fe(darcy_fe, 1, saturation_fe, 1);
2657 * DoFHandler<dim> joint_dof_handler(triangulation);
2658 * joint_dof_handler.distribute_dofs(joint_fe);
2659 * Assert(joint_dof_handler.n_dofs() ==
2660 * darcy_dof_handler.n_dofs() + saturation_dof_handler.n_dofs(),
2661 * ExcInternalError());
2663 * Vector<double> joint_solution(joint_dof_handler.n_dofs());
2666 * std::vector<types::global_dof_index> local_joint_dof_indices(
2667 * joint_fe.dofs_per_cell);
2668 * std::vector<types::global_dof_index> local_darcy_dof_indices(
2669 * darcy_fe.dofs_per_cell);
2670 * std::vector<types::global_dof_index> local_saturation_dof_indices(
2671 * saturation_fe.dofs_per_cell);
2673 * auto joint_cell = joint_dof_handler.begin_active();
2674 * const auto joint_endc = joint_dof_handler.end();
2675 * auto darcy_cell = darcy_dof_handler.begin_active();
2676 * auto saturation_cell = saturation_dof_handler.begin_active();
2678 * for (; joint_cell != joint_endc;
2679 * ++joint_cell, ++darcy_cell, ++saturation_cell)
2681 * joint_cell->get_dof_indices(local_joint_dof_indices);
2682 * darcy_cell->get_dof_indices(local_darcy_dof_indices);
2683 * saturation_cell->get_dof_indices(local_saturation_dof_indices);
2685 * for (unsigned int i = 0; i < joint_fe.dofs_per_cell; ++i)
2686 * if (joint_fe.system_to_base_index(i).first.first == 0)
2688 * Assert(joint_fe.system_to_base_index(i).second <
2689 * local_darcy_dof_indices.size(),
2690 * ExcInternalError());
2691 * joint_solution(local_joint_dof_indices[i]) = darcy_solution(
2692 * local_darcy_dof_indices[joint_fe.system_to_base_index(i)
2697 * Assert(joint_fe.system_to_base_index(i).first.first == 1,
2698 * ExcInternalError());
2699 * Assert(joint_fe.system_to_base_index(i).second <
2700 * local_darcy_dof_indices.size(),
2701 * ExcInternalError());
2702 * joint_solution(local_joint_dof_indices[i]) =
2703 * saturation_solution(
2704 * local_saturation_dof_indices
2705 * [joint_fe.system_to_base_index(i).second]);
2709 * std::vector<std::string> joint_solution_names(dim, "velocity");
2710 * joint_solution_names.emplace_back("pressure");
2711 * joint_solution_names.emplace_back("saturation");
2713 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2714 * data_component_interpretation(
2715 * dim, DataComponentInterpretation::component_is_part_of_vector);
2716 * data_component_interpretation.push_back(
2717 * DataComponentInterpretation::component_is_scalar);
2718 * data_component_interpretation.push_back(
2719 * DataComponentInterpretation::component_is_scalar);
2721 * DataOut<dim> data_out;
2723 * data_out.attach_dof_handler(joint_dof_handler);
2724 * data_out.add_data_vector(joint_solution,
2725 * joint_solution_names,
2726 * DataOut<dim>::type_dof_data,
2727 * data_component_interpretation);
2729 * data_out.build_patches();
2731 * std::string filename =
2732 * "solution-" + Utilities::int_to_string(timestep_number, 5) + ".vtu";
2733 * std::ofstream output(filename);
2734 * data_out.write_vtu(output);
2742 * <a name="Toolfunctions"></a>
2743 * <h3>Tool functions</h3>
2748 * <a name="TwoPhaseFlowProblemdimdetermine_whether_to_solve_for_pressure_and_velocity"></a>
2749 * <h4>TwoPhaseFlowProblem<dim>::determine_whether_to_solve_for_pressure_and_velocity</h4>
2753 * This function implements the a posteriori criterion for adaptive operator
2754 * splitting. The function is relatively straightforward given the way we
2755 * have implemented other functions above and given the formula for the
2756 * criterion derived in the paper.
2760 * If one decides that one wants the original IMPES method in which the
2761 * Darcy equation is solved in every time step, then this can be achieved by
2762 * setting the threshold value <code>AOS_threshold</code> (with a default of
2763 * @f$5.0@f$) to zero, thereby forcing the function to always return true.
2767 * Finally, note that the function returns true unconditionally for the
2768 * first two time steps to ensure that we have always solved the Darcy
2769 * system at least twice when skipping its solution, thereby allowing us to
2770 * extrapolate the velocity from the last two solutions in
2771 * <code>solve()</code>.
2774 * template <int dim>
2775 * bool TwoPhaseFlowProblem<
2776 * dim>::determine_whether_to_solve_for_pressure_and_velocity() const
2778 * if (timestep_number <= 2)
2781 * const QGauss<dim> quadrature_formula(saturation_degree + 2);
2782 * const unsigned int n_q_points = quadrature_formula.size();
2784 * FEValues<dim> fe_values(saturation_fe,
2785 * quadrature_formula,
2786 * update_values | update_quadrature_points);
2788 * std::vector<double> old_saturation_after_solving_pressure(n_q_points);
2789 * std::vector<double> present_saturation(n_q_points);
2791 * std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
2793 * double max_global_aop_indicator = 0.0;
2795 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2797 * double max_local_mobility_reciprocal_difference = 0.0;
2798 * double max_local_permeability_inverse_l1_norm = 0.0;
2800 * fe_values.reinit(cell);
2801 * fe_values.get_function_values(
2802 * saturation_matching_last_computed_darcy_solution,
2803 * old_saturation_after_solving_pressure);
2804 * fe_values.get_function_values(saturation_solution, present_saturation);
2806 * k_inverse.value_list(fe_values.get_quadrature_points(),
2807 * k_inverse_values);
2809 * for (unsigned int q = 0; q < n_q_points; ++q)
2811 * const double mobility_reciprocal_difference = std::fabs(
2812 * mobility_inverse(present_saturation[q], viscosity) -
2813 * mobility_inverse(old_saturation_after_solving_pressure[q],
2816 * max_local_mobility_reciprocal_difference =
2817 * std::max(max_local_mobility_reciprocal_difference,
2818 * mobility_reciprocal_difference);
2820 * max_local_permeability_inverse_l1_norm =
2821 * std::max(max_local_permeability_inverse_l1_norm,
2822 * l1_norm(k_inverse_values[q]));
2825 * max_global_aop_indicator =
2826 * std::max(max_global_aop_indicator,
2827 * (max_local_mobility_reciprocal_difference *
2828 * max_local_permeability_inverse_l1_norm));
2831 * return (max_global_aop_indicator > AOS_threshold);
2839 * <a name="TwoPhaseFlowProblemdimproject_back_saturation"></a>
2840 * <h4>TwoPhaseFlowProblem<dim>::project_back_saturation</h4>
2844 * The next function simply makes sure that the saturation values always
2845 * remain within the physically reasonable range of @f$[0,1]@f$. While the
2846 * continuous equations guarantee that this is so, the discrete equations
2847 * don't. However,
if we allow the discrete solution to
escape this range we
2848 * get into trouble because terms like @f$F(S)@f$ and @f$F
'(S)@f$ will produce
2849 * unreasonable results (e.g. @f$F'(S)<0@f$
for @f$S<0@f$, which would imply that
2850 * the wetting fluid phase flows <i>against</i> the direction of the bulk
2851 * fluid velocity)). Consequently, at the
end of each time step, we simply
2852 *
project the saturation field back into the physically reasonable region.
2855 * template <int dim>
2856 *
void TwoPhaseFlowProblem<dim>::project_back_saturation()
2858 *
for (
unsigned int i = 0; i < saturation_solution.size(); ++i)
2859 *
if (saturation_solution(i) < 0.2)
2860 * saturation_solution(i) = 0.2;
2861 *
else if (saturation_solution(i) > 1)
2862 * saturation_solution(i) = 1;
2870 * <a name=
"TwoPhaseFlowProblemdimget_max_u_F_prime"></a>
2871 * <h4>TwoPhaseFlowProblem<dim>::get_max_u_F_prime</h4>
2875 * Another simpler helper
function: Compute the maximum of the total
2876 * velocity times the derivative of the fraction flow
function, i.e.,
2877 * compute @f$\|\mathbf{u}
F'(S)\|_{L_\infty(\Omega)}@f$. This term is used in
2878 * both the computation of the time step as well as in normalizing the
2879 * entropy-residual term in the artificial viscosity.
2882 * template <int dim>
2883 * double TwoPhaseFlowProblem<dim>::get_max_u_F_prime() const
2885 * const QGauss<dim> quadrature_formula(darcy_degree + 2);
2886 * const unsigned int n_q_points = quadrature_formula.size();
2888 * FEValues<dim> darcy_fe_values(darcy_fe, quadrature_formula, update_values);
2889 * FEValues<dim> saturation_fe_values(saturation_fe,
2890 * quadrature_formula,
2893 * std::vector<Vector<double>> darcy_solution_values(n_q_points,
2894 * Vector<double>(dim + 1));
2895 * std::vector<double> saturation_values(n_q_points);
2897 * double max_velocity_times_dF_dS = 0;
2899 * auto cell = darcy_dof_handler.begin_active();
2900 * const auto endc = darcy_dof_handler.end();
2901 * auto saturation_cell = saturation_dof_handler.begin_active();
2902 * for (; cell != endc; ++cell, ++saturation_cell)
2904 * darcy_fe_values.reinit(cell);
2905 * saturation_fe_values.reinit(saturation_cell);
2907 * darcy_fe_values.get_function_values(darcy_solution,
2908 * darcy_solution_values);
2909 * saturation_fe_values.get_function_values(old_saturation_solution,
2910 * saturation_values);
2912 * for (unsigned int q = 0; q < n_q_points; ++q)
2914 * Tensor<1, dim> velocity;
2915 * for (unsigned int i = 0; i < dim; ++i)
2916 * velocity[i] = darcy_solution_values[q](i);
2918 * const double dF_dS =
2919 * fractional_flow_derivative(saturation_values[q], viscosity);
2921 * max_velocity_times_dF_dS =
2922 * std::max(max_velocity_times_dF_dS, velocity.norm() * dF_dS);
2926 * return max_velocity_times_dF_dS;
2933 * <a name="TwoPhaseFlowProblemdimget_extrapolated_saturation_range"></a>
2934 * <h4>TwoPhaseFlowProblem<dim>::get_extrapolated_saturation_range</h4>
2938 * For computing the stabilization term, we need to know the range of the
2939 * saturation variable. Unlike in @ref step_31 "step-31", this range is trivially bounded
2940 * by the interval @f$[0,1]@f$ but we can do a bit better by looping over a
2941 * collection of quadrature points and seeing what the values are there. If
2942 * we can, i.e., if there are at least two timesteps around, we can even
2943 * take the values extrapolated to the next time step.
2947 * As before, the function is taken with minimal modifications from @ref step_31 "step-31".
2950 * template <int dim>
2951 * std::pair<double, double>
2952 * TwoPhaseFlowProblem<dim>::get_extrapolated_saturation_range() const
2954 * const QGauss<dim> quadrature_formula(saturation_degree + 2);
2955 * const unsigned int n_q_points = quadrature_formula.size();
2957 * FEValues<dim> fe_values(saturation_fe, quadrature_formula, update_values);
2958 * std::vector<double> old_saturation_values(n_q_points);
2959 * std::vector<double> old_old_saturation_values(n_q_points);
2961 * if (timestep_number != 0)
2963 * double min_saturation = std::numeric_limits<double>::max(),
2964 * max_saturation = -std::numeric_limits<double>::max();
2966 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2968 * fe_values.reinit(cell);
2969 * fe_values.get_function_values(old_saturation_solution,
2970 * old_saturation_values);
2971 * fe_values.get_function_values(old_old_saturation_solution,
2972 * old_old_saturation_values);
2974 * for (unsigned int q = 0; q < n_q_points; ++q)
2976 * const double saturation =
2977 * (1. + time_step / old_time_step) * old_saturation_values[q] -
2978 * time_step / old_time_step * old_old_saturation_values[q];
2980 * min_saturation = std::min(min_saturation, saturation);
2981 * max_saturation = std::max(max_saturation, saturation);
2985 * return std::make_pair(min_saturation, max_saturation);
2989 * double min_saturation = std::numeric_limits<double>::max(),
2990 * max_saturation = -std::numeric_limits<double>::max();
2992 * for (const auto &cell : saturation_dof_handler.active_cell_iterators())
2994 * fe_values.reinit(cell);
2995 * fe_values.get_function_values(old_saturation_solution,
2996 * old_saturation_values);
2998 * for (unsigned int q = 0; q < n_q_points; ++q)
3000 * const double saturation = old_saturation_values[q];
3002 * min_saturation = std::min(min_saturation, saturation);
3003 * max_saturation = std::max(max_saturation, saturation);
3007 * return std::make_pair(min_saturation, max_saturation);
3016 * <a name="TwoPhaseFlowProblemdimcompute_viscosity"></a>
3017 * <h4>TwoPhaseFlowProblem<dim>::compute_viscosity</h4>
3021 * The final tool function is used to compute the artificial viscosity on a
3022 * given cell. This isn't particularly complicated
if you have the formula
3023 *
for it in front of you, and looking at the implementation in @ref step_31
"step-31". The
3024 * major difference to that tutorial program is that the velocity here is
3025 * not simply @f$\mathbf u@f$ but @f$\mathbf u
F'(S)@f$ and some of the formulas
3026 * need to be adjusted accordingly.
3029 * template <int dim>
3030 * double TwoPhaseFlowProblem<dim>::compute_viscosity(
3031 * const std::vector<double> & old_saturation,
3032 * const std::vector<double> & old_old_saturation,
3033 * const std::vector<Tensor<1, dim>> &old_saturation_grads,
3034 * const std::vector<Tensor<1, dim>> &old_old_saturation_grads,
3035 * const std::vector<Vector<double>> &present_darcy_values,
3036 * const double global_max_u_F_prime,
3037 * const double global_S_variation,
3038 * const double cell_diameter) const
3040 * const double beta = .4 * dim;
3041 * const double alpha = 1;
3043 * if (global_max_u_F_prime == 0)
3044 * return 5e-3 * cell_diameter;
3046 * const unsigned int n_q_points = old_saturation.size();
3048 * double max_residual = 0;
3049 * double max_velocity_times_dF_dS = 0;
3051 * const bool use_dF_dS = true;
3053 * for (unsigned int q = 0; q < n_q_points; ++q)
3056 * for (unsigned int d = 0; d < dim; ++d)
3057 * u[d] = present_darcy_values[q](d);
3059 * const double dS_dt = porosity *
3060 * (old_saturation[q] - old_old_saturation[q]) /
3063 * const double dF_dS = fractional_flow_derivative(
3064 * (old_saturation[q] + old_old_saturation[q]) / 2.0, viscosity);
3066 * const double u_grad_S =
3067 * u * dF_dS * (old_saturation_grads[q] + old_old_saturation_grads[q]) /
3070 * const double residual =
3071 * std::abs((dS_dt + u_grad_S) *
3072 * std::pow((old_saturation[q] + old_old_saturation[q]) / 2,
3075 * max_residual = std::max(residual, max_residual);
3076 * max_velocity_times_dF_dS =
3077 * std::max(std::sqrt(u * u) * (use_dF_dS ? std::max(dF_dS, 1.) : 1),
3078 * max_velocity_times_dF_dS);
3081 * const double c_R = 1.0;
3082 * const double global_scaling = c_R * porosity *
3083 * (global_max_u_F_prime)*global_S_variation /
3084 * std::pow(global_Omega_diameter, alpha - 2.);
3087 * (max_velocity_times_dF_dS)*std::min(cell_diameter,
3088 * std::pow(cell_diameter, alpha) *
3097 * <a name="TwoPhaseFlowProblemdimrun"></a>
3098 * <h3>TwoPhaseFlowProblem<dim>::run</h3>
3102 * This function is, besides <code>solve()</code>, the primary function of
3103 * this program as it controls the time iteration as well as when the
3104 * solution is written into output files and when to do mesh refinement.
3108 * With the exception of the startup code that loops back to the beginning
3109 * of the function through the <code>goto start_time_iteration</code> label,
3110 * everything should be relatively straightforward. In any case, it mimics
3111 * the corresponding function in @ref step_31 "step-31".
3114 * template <int dim>
3115 * void TwoPhaseFlowProblem<dim>::run()
3117 * const unsigned int initial_refinement = (dim == 2 ? 5 : 2);
3118 * const unsigned int n_pre_refinement_steps = (dim == 2 ? 3 : 2);
3121 * GridGenerator::hyper_cube(triangulation, 0, 1);
3122 * triangulation.refine_global(initial_refinement);
3123 * global_Omega_diameter = GridTools::diameter(triangulation);
3127 * unsigned int pre_refinement_step = 0;
3129 * start_time_iteration:
3131 * VectorTools::project(saturation_dof_handler,
3132 * saturation_constraints,
3133 * QGauss<dim>(saturation_degree + 2),
3134 * SaturationInitialValues<dim>(),
3135 * old_saturation_solution);
3137 * time_step = old_time_step = 0;
3138 * current_macro_time_step = old_macro_time_step = 0;
3144 * std::cout << "Timestep " << timestep_number << ": t=" << time
3145 * << ", dt=" << time_step << std::endl;
3149 * std::cout << std::endl;
3151 * if (timestep_number % 200 == 0)
3154 * if (timestep_number % 25 == 0)
3155 * refine_mesh(initial_refinement,
3156 * initial_refinement + n_pre_refinement_steps);
3158 * if ((timestep_number == 0) &&
3159 * (pre_refinement_step < n_pre_refinement_steps))
3161 * ++pre_refinement_step;
3162 * goto start_time_iteration;
3165 * time += time_step;
3166 * ++timestep_number;
3168 * old_old_saturation_solution = old_saturation_solution;
3169 * old_saturation_solution = saturation_solution;
3171 * while (time <= end_time);
3173 * } // namespace Step43
3180 * <a name="Thecodemaincodefunction"></a>
3181 * <h3>The <code>main()</code> function</h3>
3185 * The main function looks almost the same as in all other programs. The need
3186 * to initialize the MPI subsystem for a program that uses Trilinos -- even
3187 * for programs that do not actually run in parallel -- is explained in
3188 * @ref step_31 "step-31".
3191 * int main(int argc, char *argv[])
3195 * using namespace dealii;
3196 * using namespace Step43;
3198 * Utilities::MPI::MPI_InitFinalize mpi_initialization(
3199 * argc, argv, numbers::invalid_unsigned_int);
3203 * This program can only be run in serial. Otherwise, throw an exception.
3206 * AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
3208 * "This program can only be run in serial, use ./step-43"));
3210 * TwoPhaseFlowProblem<2> two_phase_flow_problem(1);
3211 * two_phase_flow_problem.run();
3213 * catch (std::exception &exc)
3215 * std::cerr << std::endl
3217 * << "----------------------------------------------------"
3219 * std::cerr << "Exception on processing: " << std::endl
3220 * << exc.what() << std::endl
3221 * << "Aborting!" << std::endl
3222 * << "----------------------------------------------------"
3229 * std::cerr << std::endl
3231 * << "----------------------------------------------------"
3233 * std::cerr << "Unknown exception!" << std::endl
3234 * << "Aborting!" << std::endl
3235 * << "----------------------------------------------------"
3243 <a name="Results"></a><h1>Results</h1>
3247 The output of this program is not really much different from that of
3248 @ref step_21 "step-21": it solves the same problem, after all. Of more importance are
3249 quantitative metrics such as the accuracy of the solution as well as
3250 the time needed to compute it. These are documented in detail in the
3251 two publications listed at the top of this page and we won't repeat
3254 That said, no tutorial program is complete without a couple of good
3255 pictures, so here is some output of a
run in 3
d:
3257 <table align=
"center" class=
"tutorial" cellspacing=
"3" cellpadding=
"3">
3260 <img src=
"https://www.dealii.org/images/steps/developer/step-43.3d.velocity.png" alt=
"">
3262 Velocity vectors of flow through the porous medium with
random
3263 permeability model. Streaming paths of high permeability and resulting
3264 high velocity are clearly visible.
3268 <img src=
"https://www.dealii.org/images/steps/developer/step-43.3d.streamlines.png" alt=
"">
3270 Streamlines colored by the saturation along the streamline path. Blue
3271 streamlines indicate low saturations, i.e., the flow along these
3272 streamlines must be slow or
else more fluid would have been
3273 transported along them. On the other hand, green paths indicate high
3274 velocities since the fluid front has already reached further into the
3281 <img src=
"https://www.dealii.org/images/steps/developer/step-43.3d.saturation.png" alt=
"">
3283 Streamlines with a
volume rendering of the saturation, showing how far
3284 the fluid front has advanced at
this time.
3288 <img src=
"https://www.dealii.org/images/steps/developer/step-43.3d.mesh.png" alt=
"">
3290 Surface of the mesh showing the adaptive refinement along the front.
3297 <a name=
"extensions"></a>
3298 <a name=
"Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
3301 The primary objection
one may have to
this program is that it is still too
3302 slow: 3
d computations on reasonably fine meshes are simply too expensive to be
3303 done routinely and with reasonably quick turn-around. This is similar to the
3304 situation we were in when we wrote @ref step_31
"step-31", from which
this program has taken
3305 much inspiration. The solution is similar as it was there as well: We need to
3306 parallelize the program in a way similar to how we derived @ref step_32
"step-32" out of
3307 @ref step_31
"step-31". In fact, all of the techniques used in @ref step_32
"step-32" would be transferable
3308 to
this program as well, making the program
run on dozens or hundreds of
3309 processors immediately.
3311 A different direction is to make the program more relevant to many other
3312 porous media applications. Specifically,
one avenue is to go to the primary
3313 user of porous media flow simulators, namely the oil industry. There,
3314 applications in
this area are dominated by multiphase flow (i.e., more than
3315 the two phases we have here), and the reactions they may have with each other
3316 (or any other way phases may exchange mass, such as through dissolution in and
3317 bubbling out of gas from the oil phase). Furthermore, the presence of gas
3318 often leads to compressibility effects of the fluid. Jointly, these effects
3319 are typically formulated in the widely-used
"black oil model". True reactions
3320 between multiple phases also play a role in oil reservoir modeling when
3321 considering controlled burns of oil in the reservoir to
raise pressure and
3322 temperature. These are much more complex problems, though, and left
for future
3325 Finally, from a mathematical perspective, we have derived the
3326 criterion
for re-computing the velocity/pressure solution at a given
3327 time step under the assumption that we want to compare the solution we
3328 would get at the current time step with that computed the last time we
3329 actually solved
this system. However, in the program, whenever we did
3330 not re-compute the solution, we didn
't just use the previously
3331 computed solution but instead extrapolated from the previous two times
3332 we solved the system. Consequently, the criterion was pessimistically
3333 stated: what we should really compare is the solution we would get at
3334 the current time step with the extrapolated one. Re-stating the
3335 theorem in this regard is left as an exercise.
3337 There are also other ways to extend the mathematical foundation of
3338 this program; for example, one may say that it isn't the velocity we
3339 care about, but in fact the saturation. Thus,
one may ask whether the
3340 criterion we use here to decide whether @f$\mathbf u@f$ needs to be
3341 recomputed is appropriate;
one may,
for example, suggest that it is
3342 also important to decide whether (and by how much) a wrong velocity
3343 field in fact affects the solution of the saturation equation. This
3344 would then naturally lead to a sensitivity analysis.
3346 From an algorithmic viewpoint, we have here used a criterion
for refinement
3347 that is often used in engineering, namely by looking at the gradient of
3348 the solution. However,
if you inspect the solution, you will find that
3349 it quickly leads to refinement almost everywhere, even in regions where it
3350 is clearly not necessary: frequently used therefore does not need to imply
3351 that it is a useful criterion to
begin with. On the other hand, replacing
3352 this criterion by a different and better
one should not be very difficult.
3354 should certainly be applicable to the current problem as well.
3357 <a name=
"PlainProg"></a>
3358 <h1> The plain program</h1>
3359 @include
"step-43.cc"