]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the beginnings of a Boussinesq flow problem
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 May 2007 02:27:16 +0000 (02:27 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 May 2007 02:27:16 +0000 (02:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@14711 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-22/Makefile [new file with mode: 0644]
deal.II/examples/step-22/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-22/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-22/step-22.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-22/Makefile b/deal.II/examples/step-22/Makefile
new file mode 100644 (file)
index 0000000..4aee8cd
--- /dev/null
@@ -0,0 +1,156 @@
+# $Id: Makefile,v 1.4 2006/02/10 17:53:05 wolf Exp $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-deal2-3d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-deal2-3d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/*/include/*/*.h)
+       @echo ============================ Remaking $@
+       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+               > $@ \
+         || (rm -f $@ ; false)
+       @if test -s $@ ; then : else rm $@ ; fi
+
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-22/doc/intro.dox b/deal.II/examples/step-22/doc/intro.dox
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/deal.II/examples/step-22/doc/results.dox b/deal.II/examples/step-22/doc/results.dox
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc
new file mode 100644 (file)
index 0000000..44b8374
--- /dev/null
@@ -0,0 +1,1170 @@
+/* $Id: step-22.cc 14685 2007-05-22 06:06:41Z bangerth $ */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2007 */
+
+/*    $Id: step-21.cc 14685 2007-05-22 06:06:41Z bangerth $       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 2007 by the deal.II authors */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the  text  and     */
+/*    further information on this license.                        */
+
+
+
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <fstream>
+#include <sstream>
+
+#include <base/tensor_function.h>
+
+using namespace dealii;
+
+
+                                 
+template <int dim>
+class BoussinesqFlowProblem 
+{
+  public:
+    BoussinesqFlowProblem (const unsigned int degree);
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void assemble_rhs_S ();
+    double get_maximal_velocity () const;
+    void solve ();
+    void project_back_saturation ();
+    void output_results () const;
+    
+    const unsigned int   degree;
+    
+    Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+
+    const unsigned int n_refinement_steps;
+    
+    double time_step;
+    unsigned int timestep_number;
+    double viscosity;    
+    BlockVector<double> solution;
+    BlockVector<double> old_solution;
+    BlockVector<double> system_rhs;
+};
+
+
+
+template <int 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 
+{
+  return 0;
+}
+
+
+
+template <int dim>
+class Buoyancy : public Function<dim> 
+{
+  public:
+    Buoyancy () : Function<dim>(dim) {}
+    
+    virtual void vector_value (const Point<dim>   &p,
+                              Vector<double>     &values) const;
+};
+
+
+
+template <int dim>
+void
+Buoyancy<dim>::vector_value (const Point<dim>  &p,
+                            Vector<double>    &values) const 
+{
+  values = 0;
+  values(dim-1) = p[dim-1];
+}
+
+
+template <int 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 
+{
+  return 0;
+}
+
+
+
+template <int 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 
+{
+  if (p[0] == 0)
+    return 1;
+  else
+    return 0;
+}
+
+
+
+
+template <int 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, 
+                               Vector<double>   &value) const;
+
+};
+
+
+template <int dim>
+double
+InitialValues<dim>::value (const Point<dim>  &p,
+                           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 
+{
+  ZeroFunction<dim>(dim+2).vector_value (p, values);
+}
+
+
+
+
+
+
+
+namespace SingleCurvingCrack
+{
+  template <int dim>
+  class KInverse : public TensorFunction<2,dim>
+  {
+    public:
+      virtual void value_list (const std::vector<Point<dim> > &points,
+                               std::vector<Tensor<2,dim> >    &values) const;
+  };
+
+
+  template <int dim>
+  void
+  KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+                             std::vector<Tensor<2,dim> >    &values) const
+  {
+    Assert (points.size() == values.size(),
+            ExcDimensionMismatch (points.size(), values.size()));
+
+    for (unsigned int p=0; p<points.size(); ++p)
+      {
+        values[p].clear ();
+
+        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;
+      }
+  }
+}
+
+
+
+namespace RandomMedium
+{
+  template <int dim>
+  class KInverse : public TensorFunction<2,dim>
+  {
+    public:
+      virtual void value_list (const std::vector<Point<dim> > &points,
+                               std::vector<Tensor<2,dim> >    &values) const;
+
+    private:
+      static std::vector<Point<dim> > centers;
+
+      static std::vector<Point<dim> > get_centers ();
+  };
+
+
+
+  template <int dim>
+  std::vector<Point<dim> >
+  KInverse<dim>::centers = KInverse<dim>::get_centers();
+
+
+  template <int dim>
+  std::vector<Point<dim> >
+  KInverse<dim>::get_centers ()
+  {
+    const unsigned int N = (dim == 2 ?
+                            40 :
+                            (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)
+        centers_list[i][d] = static_cast<double>(rand())/RAND_MAX;
+
+    return centers_list;
+  }
+
+
+
+  template <int dim>
+  void
+  KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+                             std::vector<Tensor<2,dim> >    &values) const
+  {
+    Assert (points.size() == values.size(),
+            ExcDimensionMismatch (points.size(), values.size()));
+
+    for (unsigned int p=0; p<points.size(); ++p)
+      {
+        values[p].clear ();
+
+        double permeability = 0;
+        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 mobility_inverse (const double S,
+                         const double viscosity)
+{
+  return 1.0 /(1.0/viscosity * S * S + (1-S) * (1-S));
+}
+
+double f_saturation (const double S,
+                     const double viscosity)
+{   
+  return S*S /( S * S +viscosity * (1-S) * (1-S));
+}
+
+
+
+
+
+
+template <int dim>
+Tensor<1,dim>
+extract_u (const FEValuesBase<dim> &fe_values,
+           const unsigned int i,
+           const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+
+  for (unsigned int d=0; d<dim; ++d)
+    tmp[d] = fe_values.shape_value_component (i,q,d);
+
+  return tmp;
+}
+
+
+
+template <int dim>
+double
+extract_div_u (const FEValuesBase<dim> &fe_values,
+               const unsigned int i,
+               const unsigned int q)
+{
+  double divergence = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    divergence += fe_values.shape_grad_component (i,q,d)[d];
+
+  return divergence;
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+extract_grad_s_u (const FEValuesBase<dim> &fe_values,
+                 const unsigned int i,
+                 const unsigned int q)
+{
+  Tensor<2,dim> tmp;
+
+  for (unsigned int d=0; d<dim; ++d)
+    for (unsigned int e=0; e<dim; ++e)
+      tmp[d][e] += (fe_values.shape_grad_component (i,q,d)[e] +
+                   fe_values.shape_grad_component (i,q,e)[d]) / 2;
+
+  return tmp;
+}
+
+
+  
+template <int dim>
+double extract_p (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  return fe_values.shape_value_component (i,q,dim);
+}
+
+
+
+template <int dim>
+double extract_s (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  return fe_values.shape_value_component (i,q,dim+1);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+extract_grad_s (const FEValuesBase<dim> &fe_values,
+                const unsigned int i,
+                const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+  for (unsigned int d=0; d<dim; ++d)
+    tmp[d] = fe_values.shape_grad_component (i,q,dim+1)[d];
+
+  return tmp;
+}
+
+
+
+
+template <class Matrix>
+class InverseMatrix : public Subscriptor
+{
+  public:
+    InverseMatrix (const Matrix &m);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const Matrix> matrix;
+
+    mutable GrowingVectorMemory<> vector_memory;    
+};
+
+
+template <class Matrix>
+InverseMatrix<Matrix>::InverseMatrix (const Matrix &m)
+                :
+                matrix (&m)
+{}
+
+
+                                 
+template <class Matrix>
+void InverseMatrix<Matrix>::vmult (Vector<double>       &dst,
+                                   const Vector<double> &src) const
+{
+  SolverControl solver_control (src.size(), 1e-8*src.l2_norm());
+  SolverCG<> cg (solver_control, vector_memory);
+
+  dst = 0;
+  
+  cg.solve (*matrix, dst, src, PreconditionIdentity());        
+}
+
+
+
+class SchurComplement : public Subscriptor
+{
+  public:
+    SchurComplement (const BlockSparseMatrix<double> &A,
+                     const InverseMatrix<SparseMatrix<double> > &Minv);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double> > > m_inverse;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+
+SchurComplement::
+SchurComplement (const BlockSparseMatrix<double> &A,
+                 const InverseMatrix<SparseMatrix<double> > &Minv)
+                :
+                system_matrix (&A),
+                m_inverse (&Minv),
+                tmp1 (A.block(0,0).m()),
+                tmp2 (A.block(0,0).m())
+{}
+
+
+void SchurComplement::vmult (Vector<double>       &dst,
+                             const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  m_inverse->vmult (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+
+class ApproximateSchurComplement : public Subscriptor
+{
+  public:
+    ApproximateSchurComplement (const BlockSparseMatrix<double> &A);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+ApproximateSchurComplement::
+ApproximateSchurComplement (const BlockSparseMatrix<double> &A)
+                :
+                system_matrix (&A),
+                tmp1 (A.block(0,0).m()),
+                tmp2 (A.block(0,0).m())
+{}
+
+
+void ApproximateSchurComplement::vmult (Vector<double>       &dst,
+                                        const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+
+
+
+
+
+template <int dim>
+BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
+                :
+                degree (degree),
+                fe (FE_Q<dim>(degree+1), dim,
+                    FE_DGQ<dim>(degree), 1,
+                    FE_DGQ<dim>(degree), 1),
+                dof_handler (triangulation),
+                n_refinement_steps (5),
+                time_step (0),
+                viscosity (0.2)
+{}
+
+
+
+
+template <int dim>
+void BoussinesqFlowProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, 0, 1);  
+  triangulation.refine_global (n_refinement_steps);
+  
+  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);  
+  const unsigned int n_u = dofs_per_component[0] * dim,
+                     n_p = dofs_per_component[dim],
+                     n_s = dofs_per_component[dim+1];
+
+  std::cout << "Number of active cells: "
+            << triangulation.n_active_cells()
+            << std::endl
+            << "Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << " (" << 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(2,0).reinit (n_s, n_u, n_couplings);
+  sparsity_pattern.block(0,1).reinit (n_u, n_p, n_couplings);
+  sparsity_pattern.block(1,1).reinit (n_p, n_p, n_couplings);
+  sparsity_pattern.block(2,1).reinit (n_s, n_p, 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);
+  system_rhs.block(2).reinit (n_s);
+  system_rhs.collect_sizes ();
+}
+
+
+template <int dim>
+double
+scalar_product (const Tensor<2,dim> &a,
+               const Tensor<2,dim> &b)
+{
+  double tmp;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      tmp += a[i][j] * b[i][j];
+  return tmp;
+}
+
+
+
+template <int dim>
+void BoussinesqFlowProblem<dim>::assemble_system () 
+{  
+  system_matrix=0;
+  system_rhs=0;
+
+  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_q_points  | update_JxW_values);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                    update_values    | update_normal_vectors |
+                                    update_q_points  | update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  
+  const unsigned int   n_q_points      = quadrature_formula.n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+
+  FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  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 Buoyancy<dim>               buoyancy;
+  const PressureBoundaryValues<dim> pressure_boundary_values;
+  const RandomMedium::KInverse<dim> k_inverse;   
+  
+  std::vector<double>               pressure_rhs_values (n_q_points);
+  std::vector<Vector<double> >      buoyancy_values (n_q_points,
+                                                    Vector<double>(dim));
+  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));
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    { 
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+
+      fe_values.get_function_values (old_solution, old_solution_values);
+
+      pressure_right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                           pressure_rhs_values);
+      buoyancy.vector_value_list (fe_values.get_quadrature_points(),
+                                 buoyancy_values);
+      k_inverse.value_list (fe_values.get_quadrature_points(),
+                            k_inverse_values);
+
+      for (unsigned int q=0; q<n_q_points; ++q)            
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+
+            const Tensor<1,dim> phi_i_u      = extract_u (fe_values, i, q);
+            const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q);
+            const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
+            const double        phi_i_p      = extract_p (fe_values, i, q);
+            const double        phi_i_s      = extract_s (fe_values, i, q); 
+            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+            
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              {
+                const Tensor<2,dim> phi_j_grads_u     = extract_grad_s_u (fe_values, j, q);
+                const double        div_phi_j_u = extract_div_u (fe_values, j, q);
+                const double        phi_j_p     = extract_p (fe_values, j, q);
+                const double        phi_j_s     = extract_s (fe_values, j, q);
+                
+                local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_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);     
+              }
+
+           Assert (dim == 2, ExcInternalError());
+            local_rhs(i) += (-phi_i_p * pressure_rhs_values[q]
+                            +phi_i_u[0] * buoyancy_values[q](0)
+                            +phi_i_u[1] * buoyancy_values[q](1))*
+                            fe_values.JxW(q);
+          }
+      
+
+      for (unsigned int face_no=0;
+           face_no<GeometryInfo<dim>::faces_per_cell;
+           ++face_no)
+        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 i=0; i<dofs_per_cell; ++i)
+                {
+                  const Tensor<1,dim>
+                    phi_i_u = extract_u (fe_face_values, i, q);
+
+                  local_rhs(i) += -(phi_i_u *
+                                    fe_face_values.normal_vector(q) *
+                                    boundary_values[q] *
+                                    fe_face_values.JxW(q));
+                }
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+         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);
+    }
+  
+  {
+    std::vector<bool> component_mask (dim+2, true);
+    component_mask[dim] = component_mask[dim+1] = false;
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             ZeroFunction<dim>(dim+2),
+                                             boundary_values,
+                                             component_mask);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       system_matrix,
+                                       solution,
+                                       system_rhs);  
+  }
+}
+
+
+
+
+
+
+template <int dim>
+void BoussinesqFlowProblem<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_q_points  | update_JxW_values);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                    update_values    | update_normal_vectors |
+                                    update_q_points  | update_JxW_values);
+  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.n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+  
+  Vector<double>       local_rhs (dofs_per_cell);
+
+  std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > old_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > old_solution_values_face_neighbor(n_face_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > present_solution_values(n_q_points, Vector<double>(dim+2));
+  std::vector<Vector<double> > present_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+
+  std::vector<double> neighbor_saturation (n_face_q_points);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  SaturationBoundaryValues<dim> saturation_boundary_values;
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      local_rhs = 0;
+      fe_values.reinit (cell);
+
+      fe_values.get_function_values (old_solution, old_solution_values);
+      fe_values.get_function_values (solution, present_solution_values);
+
+      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);
+            Tensor<1,dim> present_u;
+            for (unsigned int d=0; d<dim; ++d)
+              present_u[d] = present_solution_values[q](d);
+
+            const double        phi_i_s      = extract_s(fe_values, i, q);
+            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+                     
+            local_rhs(i) += (time_step *
+                             f_saturation(old_s,viscosity) *
+                             present_u *
+                             grad_phi_i_s
+                             +
+                             old_s * phi_i_s)
+                            *
+                            fe_values.JxW(q);
+          }
+
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
+           ++face_no)
+        {
+          fe_face_values.reinit (cell, face_no);
+
+          fe_face_values.get_function_values (old_solution, old_solution_values_face);
+          fe_face_values.get_function_values (solution, present_solution_values_face);
+
+          if (cell->at_boundary(face_no))
+            saturation_boundary_values
+              .value_list (fe_face_values.get_quadrature_points(),
+                           neighbor_saturation);
+          else
+            {
+              const typename DoFHandler<dim>::active_cell_iterator
+                neighbor = cell->neighbor(face_no);
+              const unsigned int
+                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)
+            {
+              Tensor<1,dim> present_u_face;
+              for (unsigned int d=0; d<dim; ++d)
+                present_u_face[d] = present_solution_values_face[q](d);
+
+              const double normal_flux = present_u_face *
+                                         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 *
+                                f_saturation((is_outflow_q_point == true
+                                              ?
+                                              old_solution_values_face[q](dim+1)
+                                              :
+                                              neighbor_saturation[q]),
+                                             viscosity) *
+                                extract_s(fe_face_values,i,q) *
+                                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);       
+    }
+} 
+
+
+
+
+template <int dim>
+void BoussinesqFlowProblem<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());
+  
+
+  {
+    m_inverse.vmult (tmp, system_rhs.block(0));
+    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 (system_matrix.block(0,0).m(),
+                                  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."
+              << std::endl;
+  }
+
+  {
+    system_matrix.block(0,1).vmult (tmp, solution.block(1));
+    tmp *= -1;
+    tmp += system_rhs.block(0);
+    
+    m_inverse.vmult (solution.block(0), tmp);
+  }
+
+  time_step = std::pow(0.5, double(n_refinement_steps)) /
+              std::max(get_maximal_velocity(),1.);
+
+  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;             
+  } 
+
+   
+  old_solution = solution; 
+}
+                                 
+
+
+template <int dim>
+void BoussinesqFlowProblem<dim>::output_results ()  const
+{
+  if (timestep_number % 5 != 0)
+    return;
+  
+  std::vector<std::string> solution_names;
+  switch (dim)
+    {
+      case 2:
+            solution_names.push_back ("u");
+            solution_names.push_back ("v");
+            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 ("w");
+            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";
+
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+}
+
+
+
+
+template <int dim>
+void
+BoussinesqFlowProblem<dim>::project_back_saturation ()
+{
+  for (unsigned int i=0; i<solution.block(2).size(); ++i)
+    if (solution.block(2)(i) < 0)
+      solution.block(2)(i) = 0;
+    else
+      if (solution.block(2)(i) > 1)
+        solution.block(2)(i) = 1;
+}
+
+
+
+template <int dim>
+double
+BoussinesqFlowProblem<dim>::get_maximal_velocity () const
+{
+  QGauss<dim>   quadrature_formula(degree+2); 
+  const unsigned int   n_q_points
+    = quadrature_formula.n_quadrature_points;
+
+  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();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      fe_values.get_function_values (solution, solution_values);
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          Tensor<1,dim> velocity;
+          for (unsigned int i=0; i<dim; ++i)
+            velocity[i] = solution_values[q](i);          
+          
+          max_velocity = std::max (max_velocity,
+                                   velocity.norm());
+        }
+    }
+
+  return max_velocity;
+}
+
+
+
+template <int dim>
+void BoussinesqFlowProblem<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; 
+
+      assemble_system ();
+
+      solve ();
+      
+      output_results ();
+
+      time += time_step;
+      ++timestep_number;
+      std::cout << "   Now at t=" << time
+                << ", dt=" << time_step << '.'
+                << std::endl
+                << std::endl;
+    }
+  while (time <= 250);
+}
+
+    
+
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      BoussinesqFlowProblem<2> flow_problem(0);
+      flow_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.