From: heister Date: Fri, 8 Feb 2013 16:56:00 +0000 (+0000) Subject: add benchmark test X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bca92a3879dd9df3876050f965447acf8bee27cc;p=dealii-svn.git add benchmark test git-svn-id: https://svn.dealii.org/trunk@28279 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/benchmarks/test_assembly/Makefile b/tests/benchmarks/test_assembly/Makefile new file mode 100644 index 0000000000..f35c808260 --- /dev/null +++ b/tests/benchmarks/test_assembly/Makefile @@ -0,0 +1,144 @@ +# $Id: Makefile 25724 2012-07-24 23:35:36Z bangerth $ + + +# 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 = step-22 + +# 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 = off + + +# 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 = ../deal.II + + +# 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 *vtk *ucd *.d2 + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tells 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. deal.II has two +# libraries: one for the debug mode version of the +# application and one for optimized mode. +libs.g := $(lib-deal2.g) +libs.o := $(lib-deal2.o) + + +# We now use the variable defined above to 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 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. +all: $(target)$(EXEEXT) +$(target)$(EXEEXT) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@ $^ $(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)$(EXEEXT) + @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========= $( $@" + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +./%.$(OBJEXT) : + @echo "==============optimized===== $( $@" + @$(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: all 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 creates 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/deal.II/*/*.h) + @echo ============================ Remaking $@ + @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \ + > $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then true ; else rm $@ ; false ; fi + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/tests/benchmarks/test_assembly/step-22.cc b/tests/benchmarks/test_assembly/step-22.cc new file mode 100644 index 0000000000..7ae7707e5a --- /dev/null +++ b/tests/benchmarks/test_assembly/step-22.cc @@ -0,0 +1,1360 @@ +/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */ + +/* $Id: step-22.cc 25840 2012-08-09 20:22:00Z bangerth $ */ +/* */ +/* Copyright (C) 2008, 2009, 2010, 2011, 2012 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. */ + + + // @sect3{Include files} + + // As usual, we start by including + // some well-known files: +#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 + + // Then we need to include the header file + // for the sparse direct solver UMFPACK: +#include +#include + + // This includes the library for the + // incomplete LU factorization that will + // be used as a preconditioner in 3D: +#include + + // This is C++: +#include +#include + + // As in all programs, the namespace dealii + // is included: +namespace Step22 +{ + using namespace dealii; + + // @sect3{Defining the inner preconditioner type} + + // As explained in the introduction, we are + // going to use different preconditioners for + // two and three space dimensions, + // respectively. We distinguish between + // them by the use of the spatial dimension + // as a template parameter. See step-4 for + // details on templates. We are not going to + // create any preconditioner object here, all + // we do is to create class that holds a + // local typedef determining the + // preconditioner class so we can write our + // program in a dimension-independent way. + template + struct InnerPreconditioner; + + // In 2D, we are going to use a sparse direct + // solver as preconditioner: + template <> + struct InnerPreconditioner<2> + { + typedef SparseILU type; +// typedef SparseDirectUMFPACK type; + }; + + // And the ILU preconditioning in 3D, called + // by SparseILU: + template <> + struct InnerPreconditioner<3> + { + typedef SparseILU type; + }; + + + // @sect3{The StokesProblem class template} + + // This is an adaptation of step-20, so the + // main class and the data types are the + // same as used there. In this example we + // also use adaptive grid refinement, which + // is handled in analogy to + // step-6. According to the discussion in + // the introduction, we are also going to + // use the ConstraintMatrix for + // implementing Dirichlet boundary + // conditions. Hence, we change the name + // hanging_node_constraints + // into constraints. + template + class StokesProblem + { + public: + StokesProblem (const unsigned int degree); + void run (); + + private: + void setup_dofs (); + void assemble_system (); + void solve (); + void output_results (const unsigned int refinement_cycle) const; + void refine_mesh (); + + const unsigned int degree; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + + // This one is new: We shall use a + // so-called shared pointer structure to + // access the preconditioner. Shared + // pointers are essentially just a + // convenient form of pointers. Several + // shared pointers can point to the same + // object (just like regular pointers), + // but when the last shared pointer + // object to point to a preconditioner + // object is deleted (for example if a + // shared pointer object goes out of + // scope, if the class of which it is a + // member is destroyed, or if the pointer + // is assigned a different preconditioner + // object) then the preconditioner object + // pointed to is also destroyed. This + // ensures that we don't have to manually + // track in how many places a + // preconditioner object is still + // referenced, it can never create a + // memory leak, and can never produce a + // dangling pointer to an already + // destroyed object: + std_cxx1x::shared_ptr::type> A_preconditioner; + + TimerOutput timer; + }; + + // @sect3{Boundary values and right hand side} + + // As in step-20 and most other + // example programs, the next task is + // to define the data for the PDE: + // For the Stokes problem, we are + // going to use natural boundary + // values on parts of the boundary + // (i.e. homogenous Neumann-type) for + // which we won't have to do anything + // special (the homogeneity implies + // that the corresponding terms in + // the weak form are simply zero), + // and boundary conditions on the + // velocity (Dirichlet-type) on the + // rest of the boundary, as described + // in the introduction. + // + // In order to enforce the Dirichlet + // boundary values on the velocity, + // we will use the + // VectorTools::interpolate_boundary_values + // function as usual which requires + // us to write a function object with + // as many components as the finite + // element has. In other words, we + // have to define the function on the + // $(u,p)$-space, but we are going to + // filter out the pressure component + // when interpolating the boundary + // values. + + // The following function object is a + // representation of the boundary + // values described in the + // introduction: + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + }; + + + template + double + BoundaryValues::value (const Point &p, + const unsigned int component) const + { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + + if (component == 0) + return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0)); + return 0; + } + + + template + void + BoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + + + + // We implement similar functions for + // the right hand side which for the + // current example is simply zero: + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + + }; + + + template + double + RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + template + void + RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } + + + // @sect3{Linear solvers and preconditioners} + + // The linear solvers and preconditioners are + // discussed extensively in the + // introduction. Here, we create the + // respective objects that will be used. + + // @sect4{The InverseMatrix class template} + + // The InverseMatrix + // class represents the data + // structure for an inverse + // matrix. It is derived from the one + // in step-20. The only difference is + // that we now do include a + // preconditioner to the matrix since + // we will apply this class to + // different kinds of matrices that + // will require different + // preconditioners (in step-20 we did + // not use a preconditioner in this + // class at all). The types of matrix + // and preconditioner are passed to + // this class via template + // parameters, and matrix and + // preconditioner objects of these + // types will then be passed to the + // constructor when an + // InverseMatrix object + // is created. The member function + // vmult is, as in + // step-20, a multiplication with a + // vector, obtained by solving a + // linear system: + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + const SmartPointer preconditioner; + }; + + + template + InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (&preconditioner) + {} + + + // This is the implementation of the + // vmult function. + + // In this class we use a rather large + // tolerance for the solver control. The + // reason for this is that the function is + // used very frequently, and hence, any + // additional effort to make the residual + // in the CG solve smaller makes the + // solution more expensive. Note that we do + // not only use this class as a + // preconditioner for the Schur complement, + // but also when forming the inverse of the + // Laplace matrix – which is hence + // directly responsible for the accuracy of + // the solution itself, so we can't choose + // a too large tolerance, either. + template + void InverseMatrix::vmult (Vector &dst, + const Vector &src) const + { + SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); + SolverCG<> cg (solver_control); + + dst = 0; + + cg.solve (*matrix, dst, src, *preconditioner); + } + + + // @sect4{The SchurComplement class template} + + // This class implements the Schur complement + // discussed in the introduction. It is in + // analogy to step-20. Though, we now call + // it with a template parameter + // Preconditioner in order to + // access that when specifying the respective + // type of the inverse matrix class. As a + // consequence of the definition above, the + // declaration InverseMatrix now + // contains the second template parameter + // for a preconditioner class as above, which + // affects the SmartPointer + // object m_inverse as well. + template + class SchurComplement : public Subscriptor + { + public: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix, Preconditioner> &A_inverse); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer, Preconditioner> > A_inverse; + + mutable Vector tmp1, tmp2; + }; + + + + template + SchurComplement:: + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix,Preconditioner> &A_inverse) + : + system_matrix (&system_matrix), + A_inverse (&A_inverse), + tmp1 (system_matrix.block(0,0).m()), + tmp2 (system_matrix.block(0,0).m()) + {} + + + template + void SchurComplement::vmult (Vector &dst, + const Vector &src) const + { + system_matrix->block(0,1).vmult (tmp1, src); + A_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); + } + + + // @sect3{StokesProblem class implementation} + + // @sect4{StokesProblem::StokesProblem} + + // The constructor of this class + // looks very similar to the one of + // step-20. The constructor + // initializes the variables for the + // polynomial degree, triangulation, + // finite element system and the dof + // handler. The underlying polynomial + // functions are of order + // degree+1 for the + // vector-valued velocity components + // and of order degree + // for the pressure. This gives the + // LBB-stable element pair + // $Q_{degree+1}^d\times Q_{degree}$, + // often referred to as the + // Taylor-Hood element. + // + // Note that we initialize the triangulation + // with a MeshSmoothing argument, which + // ensures that the refinement of cells is + // done in a way that the approximation of + // the PDE solution remains well-behaved + // (problems arise if grids are too + // unstructered), see the documentation of + // Triangulation::MeshSmoothing + // for details. + template + StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + triangulation (Triangulation::maximum_smoothing), + fe (FE_Q(degree+1), dim, + FE_Q(degree), 1), + dof_handler (triangulation), + timer (std::cout, TimerOutput::summary, TimerOutput::cpu_times) + {} + + + // @sect4{StokesProblem::setup_dofs} + + // Given a mesh, this function + // associates the degrees of freedom + // with it and creates the + // corresponding matrices and + // vectors. At the beginning it also + // releases the pointer to the + // preconditioner object (if the + // shared pointer pointed at anything + // at all at this point) since it + // will definitely not be needed any + // more after this point and will + // have to be re-computed after + // assembling the matrix, and unties + // the sparse matrix from its + // sparsity pattern object. + // + // We then proceed with distributing + // degrees of freedom and renumbering + // them: In order to make the ILU + // preconditioner (in 3D) work + // efficiently, it is important to + // enumerate the degrees of freedom + // in such a way that it reduces the + // bandwidth of the matrix, or maybe + // more importantly: in such a way + // that the ILU is as close as + // possible to a real LU + // decomposition. On the other hand, + // we need to preserve the block + // structure of velocity and pressure + // already seen in in step-20 and + // step-21. This is done in two + // steps: First, all dofs are + // renumbered to improve the ILU and + // then we renumber once again by + // components. Since + // DoFRenumbering::component_wise + // does not touch the renumbering + // within the individual blocks, the + // basic renumbering from the first + // step remains. As for how the + // renumber degrees of freedom to + // improve the ILU: deal.II has a + // number of algorithms that attempt + // to find orderings to improve ILUs, + // or reduce the bandwidth of + // matrices, or optimize some other + // aspect. The DoFRenumbering + // namespace shows a comparison of + // the results we obtain with several + // of these algorithms based on the + // testcase discussed here in this + // tutorial program. Here, we will + // use the traditional Cuthill-McKee + // algorithm already used in some of + // the previous tutorial programs. + // In the + // section on improved ILU + // we're going to discuss this issue + // in more detail. + + // There is one more change compared + // to previous tutorial programs: + // There is no reason in sorting the + // dim velocity + // components individually. In fact, + // rather than first enumerating all + // $x$-velocities, then all + // $y$-velocities, etc, we would like + // to keep all velocities at the same + // location together and only + // separate between velocities (all + // components) and pressures. By + // default, this is not what the + // DoFRenumbering::component_wise + // function does: it treats each + // vector component separately; what + // we have to do is group several + // components into "blocks" and pass + // this block structure to that + // function. Consequently, we + // allocate a vector + // block_component with + // as many elements as there are + // components and describe all + // velocity components to correspond + // to block 0, while the pressure + // component will form block 1: + template + void StokesProblem::setup_dofs () + { + A_preconditioner.reset (); + system_matrix.clear (); + + timer.enter_section("distribute dofs"); + dof_handler.distribute_dofs (fe); + timer.exit_section("distribute dofs"); + + timer.enter_section("renumbering"); + DoFRenumbering::Cuthill_McKee (dof_handler); + + std::vector block_component (dim+1,0); + block_component[dim] = 1; + DoFRenumbering::component_wise (dof_handler, block_component); + timer.exit_section("renumbering"); + + timer.enter_section("setup constraints"); + { + constraints.clear (); + std::vector component_mask (dim+1, true); + component_mask[dim] = false; + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 1, + BoundaryValues(), + constraints, + component_mask); + } + + constraints.close (); + timer.exit_section("setup constraints"); + + // In analogy to step-20, we count the dofs + // in the individual components. We could + // do this in the same way as there, but we + // want to operate on the block structure + // we used already for the renumbering: The + // function + // DoFTools::count_dofs_per_block + // does the same as + // DoFTools::count_dofs_per_component, + // but now grouped as velocity and pressure + // block via block_component. + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); + const unsigned int n_u = dofs_per_block[0], + n_p = dofs_per_block[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 << ')' + << std::endl; + + // The next task is to allocate a + // sparsity pattern for the system matrix + // we will create. We could do this in + // the same way as in step-20, + // i.e. directly build an object of type + // SparsityPattern through + // DoFTools::make_sparsity_pattern. However, + // there is a major reason not to do so: + // In 3D, the function + // DoFTools::max_couplings_between_dofs + // yields a conservative but rather large + // number for the coupling between the + // individual dofs, so that the memory + // initially provided for the creation of + // the sparsity pattern of the matrix is + // far too much -- so much actually that + // the initial sparsity pattern won't + // even fit into the physical memory of + // most systems already for + // moderately-sized 3D problems, see also + // the discussion in step-18. Instead, + // we first build a temporary object that + // uses a different data structure that + // doesn't require allocating more memory + // than necessary but isn't suitable for + // use as a basis of SparseMatrix or + // BlockSparseMatrix objects; in a second + // step we then copy this object into an + // object of BlockSparsityPattern. This + // is entirely analgous to what we + // already did in step-11 and step-18. + // + // There is one snag again here, though: + // it turns out that using the + // CompressedSparsityPattern (or the + // block version + // BlockCompressedSparsityPattern we + // would use here) has a bottleneck that + // makes the algorithm to build the + // sparsity pattern be quadratic in the + // number of degrees of freedom. This + // doesn't become noticeable until we get + // well into the range of several 100,000 + // degrees of freedom, but eventually + // dominates the setup of the linear + // system when we get to more than a + // million degrees of freedom. This is + // due to the data structures used in the + // CompressedSparsityPattern class, + // nothing that can easily be + // changed. Fortunately, there is an easy + // solution: the + // CompressedSimpleSparsityPattern class + // (and its block variant + // BlockCompressedSimpleSparsityPattern) + // has exactly the same interface, uses a + // different %internal data structure and + // is linear in the number of degrees of + // freedom and therefore much more + // efficient for large problems. As + // another alternative, we could also + // have chosen the class + // BlockCompressedSetSparsityPattern that + // uses yet another strategy for %internal + // memory management. Though, that class + // turns out to be more memory-demanding + // than + // BlockCompressedSimpleSparsityPattern + // for this example. + // + // Consequently, this is the class that + // we will use for our intermediate + // sparsity representation. All this is + // done inside a new scope, which means + // that the memory of csp + // will be released once the information + // has been copied to + // sparsity_pattern. + { + timer.enter_section("make csp"); + BlockCompressedSimpleSparsityPattern csp (2,2); + + csp.block(0,0).reinit (n_u, n_u); + csp.block(1,0).reinit (n_p, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(1,1).reinit (n_p, n_p); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + timer.exit_section("make csp"); + timer.enter_section("copy sp"); + sparsity_pattern.copy_from (csp); + timer.exit_section("copy sp"); + } + + // Finally, the system matrix, + // solution and right hand side are + // created from the block + // structure as in step-20: + timer.enter_section("create matrix and vectors"); + system_matrix.reinit (sparsity_pattern); + + solution.reinit (2); + solution.block(0).reinit (n_u); + solution.block(1).reinit (n_p); + solution.collect_sizes (); + + system_rhs.reinit (2); + system_rhs.block(0).reinit (n_u); + system_rhs.block(1).reinit (n_p); + system_rhs.collect_sizes (); + timer.exit_section("create matrix and vectors"); + } + + + // @sect4{StokesProblem::assemble_system} + + // The assembly process follows the + // discussion in step-20 and in the + // introduction. We use the well-known + // abbreviations for the data structures + // that hold the local matrix, right + // hand side, and global + // numbering of the degrees of freedom + // for the present cell. + template + void StokesProblem::assemble_system () + { + system_matrix=0; + system_rhs=0; + + QGauss quadrature_formula(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + // Next, we need two objects that work as + // extractors for the FEValues + // object. Their use is explained in detail + // in the report on @ref vector_valued : + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + // As an extension over step-20 and + // step-21, we include a few + // optimizations that make assembly + // much faster for this particular + // problem. The improvements are + // based on the observation that we + // do a few calculations too many + // times when we do as in step-20: + // The symmetric gradient actually + // has dofs_per_cell + // different values per quadrature + // point, but we extract it + // dofs_per_cell*dofs_per_cell + // times from the FEValues object - + // for both the loop over + // i and the inner + // loop over j. In 3d, + // that means evaluating it + // $89^2=7921$ instead of $89$ + // times, a not insignificant + // difference. + // + // So what we're + // going to do here is to avoid + // such repeated calculations by + // getting a vector of rank-2 + // tensors (and similarly for + // the divergence and the basis + // function value on pressure) + // at the quadrature point prior + // to starting the loop over the + // dofs on the cell. First, we + // create the respective objects + // that will hold these + // values. Then, we start the + // loop over all cells and the loop + // over the quadrature points, + // where we first extract these + // values. There is one more + // optimization we implement here: + // the local matrix (as well as + // the global one) is going to + // be symmetric, since all + // the operations involved are + // symmetric with respect to $i$ + // and $j$. This is implemented by + // simply running the inner loop + // not to dofs_per_cell, + // but only up to i, + // the index of the outer loop. + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + 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; + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; q phi_p[i] * + // phi_p[j] , yielding a + // pressure mass matrix in the + // $(1,1)$ block of the matrix as + // discussed in the + // introduction. That this term only + // ends up in the $(1,1)$ block stems + // from the fact that both of the + // factors in phi_p[i] * + // phi_p[j] are only non-zero + // when all the other terms vanish + // (and the other way around). + // + // Note also that operator* is + // overloaded for symmetric + // tensors, yielding the scalar + // product between the two + // tensors in the first line of + // the local matrix + // contribution. + + // Before we can write the local data + // into the global matrix (and + // simultaneously use the + // ConstraintMatrix object to apply + // Dirichlet boundary conditions and + // eliminate hanging node + // constraints, as we discussed in + // the introduction), we have to be + // careful about one thing, + // though. We have only build up half + // of the local matrix because of + // symmetry, but we're going to save + // the full system matrix in order to + // use the standard functions for + // solution. This is done by flipping + // the indices in case we are + // pointing into the empty part of + // the local matrix. + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, local_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + + // Before we're going to solve this + // linear system, we generate a + // preconditioner for the + // velocity-velocity matrix, i.e., + // block(0,0) in the + // system matrix. As mentioned + // above, this depends on the + // spatial dimension. Since the two + // classes described by the + // InnerPreconditioner::type + // typedef have the same interface, + // we do not have to do anything + // different whether we want to use + // a sparse direct solver or an + // ILU: + std::cout << " Computing preconditioner..." << std::endl << std::flush; + + A_preconditioner + = std_cxx1x::shared_ptr::type>(new typename InnerPreconditioner::type()); + A_preconditioner->initialize (system_matrix.block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + + } + + + + // @sect4{StokesProblem::solve} + + // After the discussion in the introduction + // and the definition of the respective + // classes above, the implementation of the + // solve function is rather + // straigt-forward and done in a similar way + // as in step-20. To start with, we need an + // object of the InverseMatrix + // class that represents the inverse of the + // matrix A. As described in the + // introduction, the inverse is generated + // with the help of an inner preconditioner + // of type + // InnerPreconditioner::type. + template + void StokesProblem::solve () + { + const InverseMatrix, + typename InnerPreconditioner::type> + A_inverse (system_matrix.block(0,0), *A_preconditioner); + Vector tmp (solution.block(0).size()); + + // This is as in step-20. We generate the + // right hand side $B A^{-1} F - G$ for the + // Schur complement and an object that + // represents the respective linear + // operation $B A^{-1} B^T$, now with a + // template parameter indicating the + // preconditioner - in accordance with the + // definition of the class. + { + Vector schur_rhs (solution.block(1).size()); + A_inverse.vmult (tmp, system_rhs.block(0)); + system_matrix.block(1,0).vmult (schur_rhs, tmp); + schur_rhs -= system_rhs.block(1); + + SchurComplement::type> + schur_complement (system_matrix, A_inverse); + + // The usual control structures for + // the solver call are created... + SolverControl solver_control (solution.block(1).size(), + 1e-6*schur_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + // Now to the preconditioner to the + // Schur complement. As explained in + // the introduction, the + // preconditioning is done by a mass + // matrix in the pressure variable. It + // is stored in the $(1,1)$ block of + // the system matrix (that is not used + // anywhere else but in + // preconditioning). + // + // Actually, the solver needs to have + // the preconditioner in the form + // $P^{-1}$, so we need to create an + // inverse operation. Once again, we + // use an object of the class + // InverseMatrix, which + // implements the vmult + // operation that is needed by the + // solver. In this case, we have to + // invert the pressure mass matrix. As + // it already turned out in earlier + // tutorial programs, the inversion of + // a mass matrix is a rather cheap and + // straight-forward operation (compared + // to, e.g., a Laplace matrix). The CG + // method with ILU preconditioning + // converges in 5-10 steps, + // independently on the mesh size. + // This is precisely what we do here: + // We choose another ILU preconditioner + // and take it along to the + // InverseMatrix object via the + // corresponding template parameter. A + // CG solver is then called within the + // vmult operation of the inverse + // matrix. + // + // An alternative that is cheaper to + // build, but needs more iterations + // afterwards, would be to choose a + // SSOR preconditioner with factor + // 1.2. It needs about twice the number + // of iterations, but the costs for its + // generation are almost neglible. + SparseILU preconditioner; + preconditioner.initialize (system_matrix.block(1,1), + SparseILU::AdditionalData()); + + InverseMatrix,SparseILU > + m_inverse (system_matrix.block(1,1), preconditioner); + + // With the Schur complement and an + // efficient preconditioner at hand, we + // can solve the respective equation + // for the pressure (i.e. block 0 in + // the solution vector) in the usual + // way: + cg.solve (schur_complement, solution.block(1), schur_rhs, + m_inverse); + + // After this first solution step, the + // hanging node constraints have to be + // distributed to the solution in order + // to achieve a consistent pressure + // field. + constraints.distribute (solution); + + std::cout << " " + << solver_control.last_step() + << " outer CG Schur complement iterations for pressure" + << std::endl; + } + + // As in step-20, we finally need to + // solve for the velocity equation where + // we plug in the solution to the + // pressure equation. This involves only + // objects we already know - so we simply + // multiply $p$ by $B^T$, subtract the + // right hand side and multiply by the + // inverse of $A$. At the end, we need to + // distribute the constraints from + // hanging nodes in order to obtain a + // constistent flow field: + { + system_matrix.block(0,1).vmult (tmp, solution.block(1)); + tmp *= -1; + tmp += system_rhs.block(0); + + A_inverse.vmult (solution.block(0), tmp); + + constraints.distribute (solution); + } + } + + + // @sect4{StokesProblem::output_results} + + // The next function generates graphical + // output. In this example, we are going to + // use the VTK file format. We attach + // names to the individual variables in the + // problem: velocity to the + // dim components of velocity + // and pressure to the + // pressure. + // + // Not all visualization programs have the + // ability to group individual vector + // components into a vector to provide + // vector plots; in particular, this holds + // for some VTK-based visualization + // programs. In this case, the logical + // grouping of components into vectors + // should already be described in the file + // containing the data. In other words, + // what we need to do is provide our output + // writers with a way to know which of the + // components of the finite element + // logically form a vector (with $d$ + // components in $d$ space dimensions) + // rather than letting them assume that we + // simply have a bunch of scalar fields. + // This is achieved using the members of + // the + // DataComponentInterpretation + // namespace: as with the filename, we + // create a vector in which the first + // dim components refer to the + // velocities and are given the tag + // DataComponentInterpretation::component_is_part_of_vector; + // we finally push one tag + // DataComponentInterpretation::component_is_scalar + // to describe the grouping of the pressure + // variable. + + // The rest of the function is then + // the same as in step-20. + template + void + StokesProblem::output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (); + + std::ostringstream filename; + filename << "solution-" + << Utilities::int_to_string (refinement_cycle, 2) + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); + } + + + // @sect4{StokesProblem::refine_mesh} + + // This is the last interesting function of + // the StokesProblem class. + // As indicated by its name, it takes the + // solution to the problem and refines the + // mesh where this is needed. The procedure + // is the same as in the respective step in + // step-6, with the exception that we base + // the refinement only on the change in + // pressure, i.e., we call the Kelly error + // estimator with a mask + // object. Additionally, we do not coarsen + // the grid again: + template + void + StokesProblem::refine_mesh () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + srand (0); + for (unsigned int i=0;i component_mask (dim+1, false); + component_mask[dim] = true; + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::type(), + solution, + estimated_error_per_cell, + component_mask); +*/ + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.0); + triangulation.execute_coarsening_and_refinement (); + } + + + // @sect4{StokesProblem::run} + + // The last step in the Stokes class is, as + // usual, the function that generates the + // initial grid and calls the other + // functions in the respective order. + // + // We start off with a rectangle of size $4 + // \times 1$ (in 2d) or $4 \times 1 \times + // 1$ (in 3d), placed in $R^2/R^3$ as + // $(-2,2)\times(-1,0)$ or + // $(-2,2)\times(0,1)\times(-1,0)$, + // respectively. It is natural to start + // with equal mesh size in each direction, + // so we subdivide the initial rectangle + // four times in the first coordinate + // direction. To limit the scope of the + // variables involved in the creation of + // the mesh to the range where we actually + // need them, we put the entire block + // between a pair of braces: + template + void StokesProblem::run () + { + { + std::vector subdivisions (dim, 1); + subdivisions[0] = 4; + + const Point bottom_left = (dim == 2 ? + Point(-2,-1) : + Point(-2,0,-1)); + const Point top_right = (dim == 2 ? + Point(2,0) : + Point(2,1,0)); + + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + bottom_left, + top_right); + } + + // A boundary indicator of 1 is set to all + // boundaries that are subject to Dirichlet + // boundary conditions, i.e. to faces that + // are located at 0 in the last coordinate + // direction. See the example description + // above for details. + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->center()[dim-1] == 0) + cell->face(f)->set_all_boundary_indicators(1); + + + // We then apply an initial refinement + // before solving for the first time. In + // 3D, there are going to be more degrees + // of freedom, so we refine less there: + triangulation.refine_global (8-dim); + + // As first seen in step-6, we cycle over + // the different refinement levels and + // refine (except for the first cycle), + // setup the degrees of freedom and + // matrices, assemble, solve and create + // output: + for (unsigned int refinement_cycle = 0; refinement_cycle<1; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; + + //if (refinement_cycle > 0) + { + timer.enter_section("refine"); + refine_mesh (); + refine_mesh (); + refine_mesh (); + //triangulation.refine_global(1); + timer.exit_section("refine"); + } + + //timer.enter_section("setup"); + setup_dofs (); + //timer.exit_section("setup"); + + std::cout << " Assembling..." << std::endl << std::flush; + timer.enter_section("assembly"); + assemble_system (); + timer.exit_section("assembly"); + + std::cout << " Solving..." << std::flush; + //timer.enter_section("solver"); + //solve (); + //timer.exit_section("solver"); + + //timer.enter_section("results"); + //output_results (refinement_cycle); + //timer.exit_section("results"); + + std::cout << std::endl; + } + + } +} + + + // @sect3{The main function} + + // The main function is the same as in + // step-20. We pass the element degree as a + // parameter and choose the space dimension + // at the well-known template slot. +int main () +{ + try + { + using namespace dealii; + using namespace Step22; + + deallog.depth_console (0); + + StokesProblem<2> flow_problem(1); + 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; +} diff --git a/tests/benchmarks/testlist.sh b/tests/benchmarks/testlist.sh index 0d5e065195..92675beef1 100644 --- a/tests/benchmarks/testlist.sh +++ b/tests/benchmarks/testlist.sh @@ -1,2 +1,2 @@ #!/bin/bash -export TESTS="step-22 tablehandler" +export TESTS="step-22 tablehandler test_assembly"