/* $Id$ */
/* */
-/* Copyright (C) 2006, 2007, 2008, 2009 by the deal.II authors */
+/* Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// @sect3{The <code>TwoPhaseFlowProblem</code> class}
-
+
// This is the main class of the program. It
// is close to the one of step-20, but with a
// few additional functions:
// several of the formulas in the nonlinear
// equations.
template <int dim>
-class TwoPhaseFlowProblem
+class TwoPhaseFlowProblem
{
public:
TwoPhaseFlowProblem (const unsigned int degree);
void run ();
-
+
private:
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void project_back_saturation ();
void output_results () const;
-
+
const unsigned int degree;
-
+
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
BlockSparseMatrix<double> system_matrix;
const unsigned int n_refinement_steps;
-
+
double time_step;
unsigned int timestep_number;
- double viscosity;
-
+ double viscosity;
+
BlockVector<double> solution;
BlockVector<double> old_solution;
BlockVector<double> system_rhs;
// is fully equipped to deal with anything
// else, if this is desired:
template <int dim>
-class PressureRightHandSide : public Function<dim>
+class PressureRightHandSide : public Function<dim>
{
public:
PressureRightHandSide () : Function<dim>(1) {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
PressureRightHandSide<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
+ const unsigned int /*component*/) const
{
return 0;
}
// mentioned in the introduction, we choose a
// linear pressure field:
template <int dim>
-class PressureBoundaryValues : public Function<dim>
+class PressureBoundaryValues : public Function<dim>
{
public:
PressureBoundaryValues () : Function<dim>(1) {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
PressureBoundaryValues<dim>::value (const Point<dim> &p,
- const unsigned int /*component*/) const
+ const unsigned int /*component*/) const
{
return 1-p[0];
}
// values. This is as explained in the
// introduction:
template <int dim>
-class SaturationBoundaryValues : public Function<dim>
+class SaturationBoundaryValues : public Function<dim>
{
public:
SaturationBoundaryValues () : Function<dim>(1) {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
SaturationBoundaryValues<dim>::value (const Point<dim> &p,
- const unsigned int /*component*/) const
+ const unsigned int /*component*/) const
{
if (p[0] == 0)
return 1;
// later go back and choose a different
// function for initial values.
template <int dim>
-class InitialValues : public Function<dim>
+class InitialValues : public Function<dim>
{
public:
InitialValues () : Function<dim>(dim+2) {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
- virtual void vector_value (const Point<dim> &p,
+ virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
InitialValues<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+ const unsigned int component) const
{
return ZeroFunction<dim>(dim+2).value (p, component);
}
template <int dim>
void
InitialValues<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
+ Vector<double> &values) const
{
ZeroFunction<dim>(dim+2).vector_value (p, values);
}
:
TensorFunction<2,dim> ()
{}
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const;
};
const double distance_to_flowline
= std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0]));
-
+
const double permeability = std::max(std::exp(-(distance_to_flowline*
distance_to_flowline)
/ (0.1 * 0.1)),
0.01);
-
+
for (unsigned int d=0; d<dim; ++d)
values[p][d][d] = 1./permeability;
}
:
TensorFunction<2,dim> ()
{}
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<2,dim> > &values) const;
(dim == 3 ?
100 :
throw ExcNotImplemented()));
-
+
std::vector<Point<dim> > centers_list (N);
for (unsigned int i=0; i<N; ++i)
for (unsigned int d=0; d<dim; ++d)
for (unsigned int i=0; i<centers.size(); ++i)
permeability += std::exp(-(points[p]-centers[i]).square()
/ (0.05 * 0.05));
-
+
const double normalized_permeability
= std::min (std::max(permeability, 0.01), 4.);
-
+
for (unsigned int d=0; d<dim; ++d)
values[p][d][d] = 1./normalized_permeability;
}
double f_saturation (const double S,
const double viscosity)
-{
+{
return S*S /( S * S +viscosity * (1-S) * (1-S));
}
{}
-
+
template <class Matrix>
void InverseMatrix<Matrix>::vmult (Vector<double> &dst,
const Vector<double> &src) const
SolverCG<> cg (solver_control);
dst = 0;
-
- cg.solve (*matrix, dst, src, PreconditionIdentity());
+
+ cg.solve (*matrix, dst, src, PreconditionIdentity());
}
private:
const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
const SmartPointer<const InverseMatrix<SparseMatrix<double> > > m_inverse;
-
+
mutable Vector<double> tmp1, tmp2;
};
private:
const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-
+
mutable Vector<double> tmp1, tmp2;
};
template <int dim>
void TwoPhaseFlowProblem<dim>::make_grid_and_dofs ()
{
- GridGenerator::hyper_cube (triangulation, 0, 1);
+ GridGenerator::hyper_cube (triangulation, 0, 1);
triangulation.refine_global (n_refinement_steps);
-
- dof_handler.distribute_dofs (fe);
+
+ dof_handler.distribute_dofs (fe);
DoFRenumbering::component_wise (dof_handler);
-
+
std::vector<unsigned int> dofs_per_component (dim+2);
- DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
+ DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
const unsigned int n_u = dofs_per_component[0],
n_p = dofs_per_component[dim],
n_s = dofs_per_component[dim+1];
<< " (" << n_u << '+' << n_p << '+'<< n_s <<')'
<< std::endl
<< std::endl;
-
+
const unsigned int
n_couplings = dof_handler.max_couplings_between_dofs();
-
+
sparsity_pattern.reinit (3,3);
sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings);
sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings);
sparsity_pattern.block(0,2).reinit (n_u, n_s, n_couplings);
sparsity_pattern.block(1,2).reinit (n_p, n_s, n_couplings);
sparsity_pattern.block(2,2).reinit (n_s, n_s, n_couplings);
-
+
sparsity_pattern.collect_sizes();
-
+
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress();
-
+
system_matrix.reinit (sparsity_pattern);
-
+
solution.reinit (3);
solution.block(0).reinit (n_u);
solution.block(1).reinit (n_p);
solution.block(2).reinit (n_s);
solution.collect_sizes ();
-
+
old_solution.reinit (3);
old_solution.block(0).reinit (n_u);
old_solution.block(1).reinit (n_p);
old_solution.block(2).reinit (n_s);
old_solution.collect_sizes ();
-
+
system_rhs.reinit (3);
system_rhs.block(0).reinit (n_u);
system_rhs.block(1).reinit (n_p);
// crack permeability function is as simple
// as just changing the namespace name.
template <int dim>
-void TwoPhaseFlowProblem<dim>::assemble_system ()
-{
+void TwoPhaseFlowProblem<dim>::assemble_system ()
+{
system_matrix=0;
system_rhs=0;
- QGauss<dim> quadrature_formula(degree+2);
+ QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
- FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
+ FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values | update_normal_vectors |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
+
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
+
const PressureRightHandSide<dim> pressure_right_hand_side;
const PressureBoundaryValues<dim> pressure_boundary_values;
- const RandomMedium::KInverse<dim> k_inverse;
-
+ const RandomMedium::KInverse<dim> k_inverse;
+
std::vector<double> pressure_rhs_values (n_q_points);
std::vector<double> boundary_values (n_face_q_points);
std::vector<Tensor<2,dim> > k_inverse_values (n_q_points);
-
+
std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
std::vector<std::vector<Tensor<1,dim> > > old_solution_grads(n_q_points,
std::vector<Tensor<1,dim> > (dim+2));
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
- {
+ {
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
// self-explanatory given the explicit
// form of the bilinear form stated in
// the introduction:
- for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double old_s = old_solution_values[q](dim+1);
const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
const double div_phi_i_u = fe_values[velocities].divergence (i, q);
const double phi_i_p = fe_values[pressure].value (i, q);
- const double phi_i_s = fe_values[saturation].value (i, q);
- const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient(i, q);
-
+ const double phi_i_s = fe_values[saturation].value (i, q);
+
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q);
const double div_phi_j_u = fe_values[velocities].divergence (j, q);
const double phi_j_p = fe_values[pressure].value (j, q);
const double phi_j_s = fe_values[saturation].value (j, q);
-
+
local_matrix(i,j) += (phi_i_u * k_inverse_values[q] *
mobility_inverse(old_s,viscosity) * phi_j_u
- div_phi_i_u * phi_j_p
- phi_i_p * div_phi_j_u
+ phi_i_s * phi_j_s)
- * fe_values.JxW(q);
+ * fe_values.JxW(q);
}
local_rhs(i) += (-phi_i_p * pressure_rhs_values[q])*
fe_values.JxW(q);
}
-
+
// Next, we also have to deal with the
// pressure boundary values. This,
if (cell->at_boundary(face_no))
{
fe_face_values.reinit (cell, face_no);
-
+
pressure_boundary_values
.value_list (fe_face_values.get_quadrature_points(),
boundary_values);
- for (unsigned int q=0; q<n_face_q_points; ++q)
+ for (unsigned int q=0; q<n_face_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim>
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
local_matrix(i,j));
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
// been computed. We therefore have this
// separate function to this end.
template <int dim>
-void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
-{
- QGauss<dim> quadrature_formula(degree+2);
- QGauss<dim-1> face_quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
+void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
+{
+ QGauss<dim> quadrature_formula(degree+2);
+ QGauss<dim-1> face_quadrature_formula(degree+2);
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
- FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
+ FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values | update_normal_vectors |
update_quadrature_points | update_JxW_values);
- FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula,
+ FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula,
update_values);
-
+
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
-
+
Vector<double> local_rhs (dofs_per_cell);
std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
SaturationBoundaryValues<dim> saturation_boundary_values;
const FEValuesExtractors::Scalar saturation (dim+1);
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
// \mathbf{v}^{n+1},\nabla sigma)$,
// where $\sigma$ is the saturation
// component of the test function:
- for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double old_s = old_solution_values[q](dim+1);
const double phi_i_s = fe_values[saturation].value (i, q);
const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient (i, q);
-
+
local_rhs(i) += (time_step *
f_saturation(old_s,viscosity) *
present_u *
neighbor_face = cell->neighbor_of_neighbor(face_no);
fe_face_values_neighbor.reinit (neighbor, neighbor_face);
-
+
fe_face_values_neighbor
.get_function_values (old_solution,
old_solution_values_face_neighbor);
-
+
for (unsigned int q=0; q<n_face_q_points; ++q)
neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
}
-
+
for (unsigned int q=0; q<n_face_q_points; ++q)
{
fe_face_values.normal_vector(q);
const bool is_outflow_q_point = (normal_flux >= 0);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_rhs(i) -= time_step *
normal_flux *
fe_face_values.JxW(q);
}
}
-
+
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
- system_rhs(local_dof_indices[i]) += local_rhs(i);
+ system_rhs(local_dof_indices[i]) += local_rhs(i);
}
-}
+}
// step-20. After that, we have to deal with
// the saturation equation (see below):
template <int dim>
-void TwoPhaseFlowProblem<dim>::solve ()
+void TwoPhaseFlowProblem<dim>::solve ()
{
const InverseMatrix<SparseMatrix<double> >
m_inverse (system_matrix.block(0,0));
Vector<double> tmp (solution.block(0).size());
Vector<double> schur_rhs (solution.block(1).size());
Vector<double> tmp2 (solution.block(2).size());
-
+
// First the pressure, using the pressure
// Schur complement of the first two
system_matrix.block(1,0).vmult (schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
-
+
SchurComplement
schur_complement (system_matrix, m_inverse);
-
+
ApproximateSchurComplement
approximate_schur_complement (system_matrix);
-
+
InverseMatrix<ApproximateSchurComplement>
preconditioner (approximate_schur_complement);
-
+
SolverControl solver_control (solution.block(1).size(),
1e-12*schur_rhs.l2_norm());
SolverCG<> cg (solver_control);
cg.solve (schur_complement, solution.block(1), schur_rhs,
preconditioner);
-
+
std::cout << " "
<< solver_control.last_step()
<< " CG Schur complement iterations for pressure."
system_matrix.block(0,1).vmult (tmp, solution.block(1));
tmp *= -1;
tmp += system_rhs.block(0);
-
+
m_inverse.vmult (solution.block(0), tmp);
}
// reasonable range:
assemble_rhs_S ();
{
-
+
SolverControl solver_control (system_matrix.block(2,2).m(),
1e-8*system_rhs.block(2).l2_norm());
SolverCG<> cg (solver_control);
cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2),
PreconditionIdentity());
-
+
project_back_saturation ();
-
+
std::cout << " "
<< solver_control.last_step()
<< " CG iterations for saturation."
- << std::endl;
- }
+ << std::endl;
+ }
+
-
- old_solution = solution;
+ old_solution = solution;
}
-
+
// @sect4{TwoPhaseFlowProblem::output_results}
{
if (timestep_number % 5 != 0)
return;
-
+
std::vector<std::string> solution_names;
switch (dim)
{
solution_names.push_back ("p");
solution_names.push_back ("S");
break;
-
+
case 3:
solution_names.push_back ("u");
solution_names.push_back ("v");
solution_names.push_back ("p");
solution_names.push_back ("S");
break;
-
+
default:
Assert (false, ExcNotImplemented());
}
-
+
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names);
data_out.build_patches (degree+1);
-
+
std::ostringstream filename;
filename << "solution-" << timestep_number << ".vtk";
double
TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
{
- QGauss<dim> quadrature_formula(degree+2);
+ QGauss<dim> quadrature_formula(degree+2);
const unsigned int n_q_points
= quadrature_formula.size();
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values);
std::vector<Vector<double> > solution_values(n_q_points,
Vector<double>(dim+2));
double max_velocity = 0;
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
{
Tensor<1,dim> velocity;
for (unsigned int i=0; i<dim; ++i)
- velocity[i] = solution_values[q](i);
-
+ velocity[i] = solution_values[q](i);
+
max_velocity = std::max (max_velocity,
velocity.norm());
}
// end time of a time step only at the end of
// the time step.
template <int dim>
-void TwoPhaseFlowProblem<dim>::run ()
+void TwoPhaseFlowProblem<dim>::run ()
{
make_grid_and_dofs();
-
+
{
ConstraintMatrix constraints;
constraints.close();
-
+
VectorTools::project (dof_handler,
constraints,
QGauss<dim>(degree+2),
InitialValues<dim>(),
old_solution);
}
-
+
timestep_number = 1;
double time = 0;
-
+
do
- {
+ {
std::cout << "Timestep " << timestep_number
- << std::endl;
+ << std::endl;
assemble_system ();
solve ();
-
+
output_results ();
time += time_step;
while (time <= 250);
}
-
+
// @sect3{The <code>main</code> function}
// That's it. In the main function, we pass
// elements, i.e. $RT_0\times DQ_0 \times
// DQ_0$. The rest is as in all the other
// programs.
-int main ()
+int main ()
{
try
{
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
-
+
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"