From: bangerth Date: Sun, 27 May 2007 02:27:16 +0000 (+0000) Subject: Add the beginnings of a Boussinesq flow problem X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=35c54304fec68baa12770b68615e14fdec2afc2b;p=dealii-svn.git Add the beginnings of a Boussinesq flow problem git-svn-id: https://svn.dealii.org/trunk@14711 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-22/Makefile b/deal.II/examples/step-22/Makefile new file mode 100644 index 0000000000..4aee8cd9cc --- /dev/null +++ b/deal.II/examples/step-22/Makefile @@ -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========= $( $@ \ + || (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 index 0000000000..e69de29bb2 diff --git a/deal.II/examples/step-22/doc/results.dox b/deal.II/examples/step-22/doc/results.dox new file mode 100644 index 0000000000..e69de29bb2 diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc new file mode 100644 index 0000000000..44b8374e20 --- /dev/null +++ b/deal.II/examples/step-22/step-22.cc @@ -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 +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +using namespace dealii; + + + +template +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 triangulation; + FESystem fe; + DoFHandler dof_handler; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + const unsigned int n_refinement_steps; + + double time_step; + unsigned int timestep_number; + double viscosity; + + BlockVector solution; + BlockVector old_solution; + BlockVector system_rhs; +}; + + + +template +class PressureRightHandSide : public Function +{ + public: + PressureRightHandSide () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +double +PressureRightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + + + +template +class Buoyancy : public Function +{ + public: + Buoyancy () : Function(dim) {} + + virtual void vector_value (const Point &p, + Vector &values) const; +}; + + + +template +void +Buoyancy::vector_value (const Point &p, + Vector &values) const +{ + values = 0; + values(dim-1) = p[dim-1]; +} + + +template +class PressureBoundaryValues : public Function +{ + public: + PressureBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double +PressureBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return 0; +} + + + +template +class SaturationBoundaryValues : public Function +{ + public: + SaturationBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +double +SaturationBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + if (p[0] == 0) + return 1; + else + return 0; +} + + + + +template +class InitialValues : public Function +{ + public: + InitialValues () : Function(dim+2) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + +}; + + +template +double +InitialValues::value (const Point &p, + const unsigned int component) const +{ + return ZeroFunction(dim+2).value (p, component); +} + + +template +void +InitialValues::vector_value (const Point &p, + Vector &values) const +{ + ZeroFunction(dim+2).vector_value (p, values); +} + + + + + + + +namespace SingleCurvingCrack +{ + template + class KInverse : public TensorFunction<2,dim> + { + public: + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + }; + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p + class KInverse : public TensorFunction<2,dim> + { + public: + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + + private: + static std::vector > centers; + + static std::vector > get_centers (); + }; + + + + template + std::vector > + KInverse::centers = KInverse::get_centers(); + + + template + std::vector > + KInverse::get_centers () + { + const unsigned int N = (dim == 2 ? + 40 : + (dim == 3 ? + 100 : + throw ExcNotImplemented())); + + std::vector > centers_list (N); + for (unsigned int i=0; i(rand())/RAND_MAX; + + return centers_list; + } + + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p +Tensor<1,dim> +extract_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + + for (unsigned int d=0; d +double +extract_div_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + double divergence = 0; + for (unsigned int d=0; d +Tensor<2,dim> +extract_grad_s_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<2,dim> tmp; + + for (unsigned int d=0; d +double extract_p (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + return fe_values.shape_value_component (i,q,dim); +} + + + +template +double extract_s (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + return fe_values.shape_value_component (i,q,dim+1); +} + + + +template +Tensor<1,dim> +extract_grad_s (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + for (unsigned int d=0; d +class InverseMatrix : public Subscriptor +{ + public: + InverseMatrix (const Matrix &m); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + + mutable GrowingVectorMemory<> vector_memory; +}; + + +template +InverseMatrix::InverseMatrix (const Matrix &m) + : + matrix (&m) +{} + + + +template +void InverseMatrix::vmult (Vector &dst, + const Vector &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 &A, + const InverseMatrix > &Minv); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer > > m_inverse; + + mutable Vector tmp1, tmp2; +}; + + + +SchurComplement:: +SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix > &Minv) + : + system_matrix (&A), + m_inverse (&Minv), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void SchurComplement::vmult (Vector &dst, + const Vector &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 &A); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + + mutable Vector tmp1, tmp2; +}; + + +ApproximateSchurComplement:: +ApproximateSchurComplement (const BlockSparseMatrix &A) + : + system_matrix (&A), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void ApproximateSchurComplement::vmult (Vector &dst, + const Vector &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 +BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) + : + degree (degree), + fe (FE_Q(degree+1), dim, + FE_DGQ(degree), 1, + FE_DGQ(degree), 1), + dof_handler (triangulation), + n_refinement_steps (5), + time_step (0), + viscosity (0.2) +{} + + + + +template +void BoussinesqFlowProblem::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 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 +double +scalar_product (const Tensor<2,dim> &a, + const Tensor<2,dim> &b) +{ + double tmp; + for (unsigned int i=0; i +void BoussinesqFlowProblem::assemble_system () +{ + system_matrix=0; + system_rhs=0; + + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const PressureRightHandSide pressure_right_hand_side; + const Buoyancy buoyancy; + const PressureBoundaryValues pressure_boundary_values; + const RandomMedium::KInverse k_inverse; + + std::vector pressure_rhs_values (n_q_points); + std::vector > buoyancy_values (n_q_points, + Vector(dim)); + std::vector boundary_values (n_face_q_points); + std::vector > k_inverse_values (n_q_points); + + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + std::vector > > old_solution_grads(n_q_points, + std::vector > (dim+2)); + + typename DoFHandler::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 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 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::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 + 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 component_mask (dim+2, true); + component_mask[dim] = component_mask[dim+1] = false; + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim+2), + boundary_values, + component_mask); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + } +} + + + + + + +template +void BoussinesqFlowProblem::assemble_rhs_S () +{ + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues fe_face_values (fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_q_points | update_JxW_values); + FEFaceValues 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 local_rhs (dofs_per_cell); + + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + std::vector > old_solution_values_face(n_face_q_points, Vector(dim+2)); + std::vector > old_solution_values_face_neighbor(n_face_q_points, Vector(dim+2)); + std::vector > present_solution_values(n_q_points, Vector(dim+2)); + std::vector > present_solution_values_face(n_face_q_points, Vector(dim+2)); + + std::vector neighbor_saturation (n_face_q_points); + std::vector local_dof_indices (dofs_per_cell); + + SaturationBoundaryValues saturation_boundary_values; + + typename DoFHandler::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 present_u; + for (unsigned int d=0; d 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::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::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 present_u_face; + for (unsigned int d=0; d= 0); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i +void BoussinesqFlowProblem::solve () +{ + const InverseMatrix > + m_inverse (system_matrix.block(0,0)); + Vector tmp (solution.block(0).size()); + Vector schur_rhs (solution.block(1).size()); + Vector 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 + 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 +void BoussinesqFlowProblem::output_results () const +{ + if (timestep_number % 5 != 0) + return; + + std::vector 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 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 +void +BoussinesqFlowProblem::project_back_saturation () +{ + for (unsigned int i=0; i 1) + solution.block(2)(i) = 1; +} + + + +template +double +BoussinesqFlowProblem::get_maximal_velocity () const +{ + QGauss quadrature_formula(degree+2); + const unsigned int n_q_points + = quadrature_formula.n_quadrature_points; + + FEValues fe_values (fe, quadrature_formula, + update_values); + std::vector > solution_values(n_q_points, + Vector(dim+2)); + double max_velocity = 0; + + typename DoFHandler::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 velocity; + for (unsigned int i=0; i +void BoussinesqFlowProblem::run () +{ + make_grid_and_dofs(); + + { + ConstraintMatrix constraints; + constraints.close(); + + VectorTools::project (dof_handler, + constraints, + QGauss(degree+2), + InitialValues(), + 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; +}