more.
<h3>The test cases.</h3>
-The implementation discussed here uses and extends
-parts of the step-21, step-31 and step-33 tutorial programs of this library. We
-use the implementation of the incomplete Cholesky decomposition provided by
-the <code>Trilinos</code> library (see [Trilinos 2004]).
-We show numerical results that illustrate the
-efficiency and accuracy of our combined methods in solving the two-phase flow equations
+@note
+The implementation discussed here uses and extends
+parts of the step-21, step-31 and step-33 tutorial programs of this
+library. In particular, if you want to understand how it works, please
+consult step-21 for a discussion of the mathematical problem, and
+step-31 from which most of the implementation is derived. We will not
+discuss aspects of the implementation that have already been discussed
+in step-31.
+
+We show numerical results for some two-phase flow equations
augmented by
appropriate initial and boundary conditions in conjunction with two different choices of the
-permeability model. In the problems considered, there is no internal source term ($q=0$).
+permeability model. In the problems considered, there is no internal
+source term ($q=0$). As mentioned above, quantitative numerical
+results are presented in [Chueh, Djilali
+and Bangerth 2011].
For simplicity, we choose $\Omega=[0,1]^d,d=2,3$, though all methods (as well
as our implementation) should work equally well on general unstructured meshes.
<br> <i>Matrix Computations</i>.
<br> 3rd Edition, Johns Hopkins, 1996.
-<li>
-MA Heroux, RA Bartlett, VE Howle, RJ Hoekstra, JJ Hu, TG Kolda, RB Lehoucq,
- KR Long, RP Pawlowski, ET Phipps, AG Salinger, HK Thornquist, RS Tuminaro,
- JM Willenbring, A Williams, and KS Stanley.
-<br> An overview of the trilinos project.
-<br> <i>ACM Trans. Math. Softw.</i>, 31(3):397-423, 2005.
-
<li>
SE Buckley and MC Leverett.
<br> Mechanism of fluid displacements in sands.
-/* Author: Chih-Che Chueh, University of Victoria, 2010 */
+/* Author: Chih-Che Chueh, University of Victoria, 2010 */
/* Wolfgang Bangerth, Texas A&M University, 2010 */
/* $Id$ */
// @sect3{Include files}
- // The first step, as always, is to include
- // the functionality of these well-known
- // deal.II library files and some C++ header
+ // The first step, as always, is to
+ // include the functionality of a
+ // number of deal.II and C++ header
// files.
//
- // In this program, we use a tensor-valued
- // coefficient. Since it may have a spatial
- // dependence, we consider it a tensor-valued
- // function. The following include file
- // provides the TensorFunction class that
- // offers such functionality:
- //
- // Then we need to include some header files
- // that provide vector, matrix, and
- // preconditioner classes that implement
- // interfaces to the respective Trilinos
- // classes, which has been used in
- // step-31. In particular, we will need
- // interfaces to the matrix and vector
- // classes based on Trilinos as well as
- // Trilinos preconditioners:
- //
- // At the end of this top-matter, we import
- // all deal.II names into the global
- // namespace:
+ // The list includes some header
+ // files that provide vector, matrix,
+ // and preconditioner classes that
+ // implement interfaces to the
+ // respective Trilinos classes; some
+ // more information on these may be
+ // found in step-31.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <fstream>
#include <sstream>
+
+ // At the end of this top-matter, we
+ // open a namespace for the current
+ // project into which all the
+ // following material will go, and
+ // then import all deal.II names into
+ // this namespace:
namespace Step43
{
using namespace dealii;
- // @sect3{The InverseMatrix class template}
-
- // This part is exactly the same as that used in step-31.
-
- // @sect3{Schur complement preconditioner}
-
- // This part for the Schur complement
- // preconditioner is almost the same as that
- // used in step-31. The only difference is
- // that the original variable name
- // stokes_matrix is replaced by another name
- // darcy_matrix to satisfy our problem.
+ // @sect3{Helper classes for solvers and preconditioners}
+
+ // In this first part we define a
+ // number of classes that we need
+ // in the construction of linear
+ // solvers and
+ // preconditioners. This part is
+ // essentially the same as that
+ // used in step-31. The only
+ // difference is that the original
+ // variable name stokes_matrix is
+ // replaced by another name
+ // darcy_matrix to match our
+ // problem.
namespace LinearSolvers
{
template <class Matrix, class Preconditioner>
// @sect3{The TwoPhaseFlowProblem class}
- // The definition of the class that defines
- // the top-level logic of solving the
- // time-dependent advection-dominated
- // two-phase flow problem (or
+ // The definition of the class that
+ // defines the top-level logic of
+ // solving the time-dependent
+ // advection-dominated two-phase
+ // flow problem (or
// Buckley-Leverett problem
- // [Buckley 1942]) is mainly based on
- // three tutorial programs (step-21, step-31,
- // step-33). The main difference is that,
- // since adaptive operator splitting is
- // considered, we need a bool-type variable
- // solve_pressure_velocity_part to tell us
- // when we need to solve the pressure and
- // velocity part, need another bool-type
- // variable
- // previous_solve_pressure_velocity_part to
- // determine if we have to cumulate
- // micro-time steps that we need them to do
- // extrapolation for the total velocity, and
- // some solution vectors
+ // [Buckley 1942]) is mainly based
+ // on tutorial programs step-21 and
+ // step-33, and in particular on
+ // step-31 where we have used
+ // basically the same general
+ // structure as done here. The main
+ // difference to step-31 is that,
+ // since adaptive operator
+ // splitting is considered, we need
+ // a bool-type variable
+ // <code>solve_for_pressure_and_velocity</code> to
+ // tell us whether we need to solve
+ // the pressure and velocity part,
+ // need another bool-type variable
+ // <code>previous_solve_for_pressure_and_velocity</code>
+ // to determine if we have to
+ // cumulate micro-time steps that
+ // we need them to do extrapolation
+ // for the total velocity, and some
+ // solution vectors
// (e.g. nth_darcy_solution_after_solving_pressure_part
// and
// n_minus_oneth_darcy_solution_after_solving_pressure_part)
- // to store some solutions in previous time
- // steps after the solution of the pressure
- // and velocity part.
+ // to store some solutions in
+ // previous time steps after the
+ // solution of the pressure and
+ // velocity part.
//
// The member functions within this class
// have been named so properly so that
double get_maximal_velocity_times_dF_dS () const;
std::pair<double,double> get_extrapolated_saturation_range () const;
void solve ();
- bool determine_whether_to_solve_pressure_velocity_part () const;
+ bool determine_whether_to_solve_for_pressure_and_velocity () const;
void compute_refinement_indicators (Vector<double> &indicator) const;
void refine_grid (const Vector<double> &indicator);
void project_back_saturation ();
TrilinosWrappers::Vector nth_saturation_solution_after_solving_pressure_part;
const unsigned int n_refinement_steps;
- bool solve_pressure_velocity_part;
- bool previous_solve_pressure_velocity_part;
+ bool solve_for_pressure_and_velocity;
+ bool previous_solve_for_pressure_and_velocity;
const double saturation_level;
const double saturation_refinement_threshold;
saturation_dof_handler (triangulation),
n_refinement_steps (4),
- solve_pressure_velocity_part (false),
- previous_solve_pressure_velocity_part (false),
+ solve_for_pressure_and_velocity (false),
+ previous_solve_for_pressure_and_velocity (false),
saturation_level (2),
saturation_refinement_threshold (0.5),
// compute the present time step.
//
// Next, we need to use two bool variables
- // solve_pressure_velocity_part and
- // previous_solve_pressure_velocity_part to
+ // <code>solve_for_pressure_and_velocity</code> and
+ // <code>previous_solve_for_pressure_and_velocity</code> to
// decide whether we stop or continue
// cumulating the micro time steps for linear
// extropolations in the next iteration. With
template <int dim>
void TwoPhaseFlowProblem<dim>::solve ()
{
- solve_pressure_velocity_part = determine_whether_to_solve_pressure_velocity_part ();
+ solve_for_pressure_and_velocity = determine_whether_to_solve_for_pressure_and_velocity ();
- if ( timestep_number <= 3 || solve_pressure_velocity_part == true )
+ if ( timestep_number <= 3 || solve_for_pressure_and_velocity == true )
{
std::cout << " Solving darcy system (pressure-velocity part)..." << std::endl;
// if we haven't computed the
// velocity before, then
// extrapolate now
- if ( !(timestep_number <= 3 || solve_pressure_velocity_part == true ))
+ if ( !(timestep_number <= 3 || solve_for_pressure_and_velocity == true ))
{
darcy_solution.block(0) = nth_darcy_solution_after_solving_pressure_part.block(0);
darcy_solution.block(0).sadd (2.0, -1.0, n_minus_oneth_darcy_solution_after_solving_pressure_part.block(0) );
}
if ( timestep_number <= 3 ||
- ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == true ) )
+ ( solve_for_pressure_and_velocity == true && previous_solve_for_pressure_and_velocity == true ) )
{
n_minus_oneth_time_step = time_step;
cumulative_nth_time_step = 0.0;
}
- else if ( solve_pressure_velocity_part == true && previous_solve_pressure_velocity_part == false )
+ else if ( solve_for_pressure_and_velocity == true && previous_solve_for_pressure_and_velocity == false )
{
n_minus_oneth_time_step = cumulative_nth_time_step;
cumulative_nth_time_step = 0.0;
cumulative_nth_time_step += time_step;
}
- previous_solve_pressure_velocity_part = solve_pressure_velocity_part;
+ previous_solve_for_pressure_and_velocity = solve_for_pressure_and_velocity;
std::cout << " Solving saturation transport equation..." << std::endl;
- // @sect3{TwoPhaseFlowProblem<dim>::determine_whether_to_solve_pressure_velocity_part}
+ // @sect3{TwoPhaseFlowProblem<dim>::determine_whether_to_solve_for_pressure_and_velocity}
// This function is to implement the a
// posteriori criterion for
// method.
template <int dim>
bool
- TwoPhaseFlowProblem<dim>::determine_whether_to_solve_pressure_velocity_part () const
+ TwoPhaseFlowProblem<dim>::determine_whether_to_solve_for_pressure_and_velocity () const
{
if (timestep_number <= 3)
return true;
template <int dim>
void TwoPhaseFlowProblem<dim>::output_results () const
{
- if ( solve_pressure_velocity_part == false )
+ if ( solve_for_pressure_and_velocity == false )
return;
const FESystem<dim> joint_fe (darcy_fe, 1,
output_results ();
- solve_pressure_velocity_part = false;
+ solve_for_pressure_and_velocity = false;
if ((timestep_number == 0) &&
(pre_refinement_step < saturation_level))