From 7c2579d8a812b1970705ef2be7c64372864dea12 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 27 Oct 2008 14:57:31 +0000 Subject: [PATCH] Further cleanups and documentation work. git-svn-id: https://svn.dealii.org/trunk@17350 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 807 ++++++++++++++-------------- deal.II/examples/step-32/step-32.cc | 63 +-- 2 files changed, 444 insertions(+), 426 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index b0a50ccfe0..35a0cf3b19 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1187,9 +1187,9 @@ void BoussinesqFlowProblem::setup_dofs () // degrees of freedom, so we can let the // DoFTools::make_sparsity_pattern function // omit these entries by setting the last - // boolean flag to false. Once the - // sparsity pattern is ready, we can use it - // to initialize the Trilinos + // boolean flag to false. Once + // the sparsity pattern is ready, we can + // use it to initialize the Trilinos // matrices. Note that the Trilinos // matrices store the sparsity pattern // internally, so there is no need to keep @@ -1310,25 +1310,25 @@ void BoussinesqFlowProblem::setup_dofs () // we create data structures for the cell // matrix and the relation between local and // global DoFs. The vectors - // phi_grad_u and phi_p are - // going to hold the values of the basis - // functions in order to faster build up the - // local matrices, as was already done in - // step-22. Before we start the loop over all - // active cells, we have to specify which - // components are pressure and which are - // velocity. + // phi_grad_u and + // phi_p are going to hold the + // values of the basis functions in order to + // faster build up the local matrices, as was + // already done in step-22. Before we start + // the loop over all active cells, we have to + // specify which components are pressure and + // which are velocity. template void BoussinesqFlowProblem::assemble_stokes_preconditioner () { stokes_preconditioner_matrix = 0; - QGauss quadrature_formula(stokes_degree+2); - FEValues stokes_fe_values (stokes_fe, quadrature_formula, - update_JxW_values | - update_values | - update_gradients); + const QGauss quadrature_formula(stokes_degree+2); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_JxW_values | + update_values | + update_gradients); const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -1397,19 +1397,22 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () - // @sect4{BoussinesqFlowProblem::assemble_stokes_preconditioner} + // @sect4{BoussinesqFlowProblem::build_stokes_preconditioner} // // This function generates the inner - // preconditioners that are going to be - // used for the Schur complement block - // preconditioner. Since the - // preconditioners need only to be - // regenerated when the matrices change, - // this function does not have to do - // anything in case the matrices have not - // changed (i.e., the flag - // rebuild_stokes_preconditioner - // has the value false). + // preconditioners that are going to be used + // for the Schur complement block + // preconditioner. Since the preconditioners + // need only to be regenerated when the + // matrices change, this function does not + // have to do anything in case the matrices + // have not changed (i.e., the flag + // rebuild_stokes_preconditioner + // has the value + // false). Otherwise its first + // task is to call + // assemble_stokes_preconditioner + // to generate the preconditioner matrices. // // Next, we set up the preconditioner for // the velocity-velocity matrix @@ -1430,8 +1433,8 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () // vector component. We do this using the // function // DoFTools::extract_constant_modes, a - // function that generates a bunch of - // dim vectors, where each one + // function that generates a set of + // dim vectors, where each one // has ones in the respective component // of the vector problem and zeros // elsewhere. Hence, these are the @@ -1461,46 +1464,52 @@ BoussinesqFlowProblem::build_stokes_preconditioner () amg_data.constant_modes = constant_modes; // Next, we set some more options of the - // AMG preconditioner. In particular, + // AMG preconditioner. In particular, we // need to tell the AMG setup that we use // quadratic basis functions for the // velocity matrix (this implies more - // nonzero elements in the matrix, so - // that a more rubust algorithm needs to - // be chosen internally). Moreover, we - // want to be able to control how the - // coarsening structure is build up. The - // way AMG does this is to look which - // matrix entries are of similar size - // than the diagonal entry in order to - // algebraically build a coarse-grid - // structure. By setting the parameter - // aggregation_threshold to + // nonzero elements in the matrix, so that + // a more rubust algorithm needs to be + // chosen internally). Moreover, we want to + // be able to control how the coarsening + // structure is build up. The way AMG does + // this is to look which matrix entries are + // of similar size as the diagonal entry in + // order to algebraically build a + // coarse-grid structure. By setting the + // parameter + // aggregation_threshold to // 0.05, we specify that all entries that // are more than five precent of size of - // some diagonal pivots in that row - // should form one coarse grid - // point. This parameter is rather - // ad-hoc, and some fine-tuning of it can - // influence the performance of the - // preconditioner. As a rule of thumb, - // larger values of - // aggregation_threshold will + // some diagonal pivots in that row should + // form one coarse grid point. This + // parameter is rather ad-hoc, and some + // fine-tuning of it can influence the + // performance of the preconditioner. As a + // rule of thumb, larger values of + // aggregation_threshold will // decrease the number of iterations, but - // increase the costs per iteration. + // increase the costs per iteration. A look + // at the Trilinos documentation will + // provide more information on these + // parameters. With this data set, we then + // initialize the preconditioner with the + // matrix we want it to apply to. // - // Eventually, we initialize the - // preconditioner for the inversion of - // the pressure mass matrix. This matrix - // is symmetric and well-behaved, so we - // can chose a simple preconditioner. We - // stick with an incomple Cholesky (IC) + // Finally, we also initialize the + // preconditioner for the inversion of the + // pressure mass matrix. This matrix is + // symmetric and well-behaved, so we can + // chose a simple preconditioner. We stick + // with an incomple Cholesky (IC) // factorization preconditioner, which is - // designed for symmetric matrices. We - // wrap the preconditioners into a - // boost::shared_ptr pointer, which makes - // it easier to recreate the - // preconditioner. + // designed for symmetric matrices. We wrap + // the preconditioners into a + // boost::shared_ptr pointer, + // which makes it easier to recreate the + // preconditioner next time around since we + // do not have to care about destroying the + // previously used object. amg_data.elliptic = true; amg_data.higher_order_elements = true; amg_data.aggregation_threshold = 5e-2; @@ -1528,53 +1537,50 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // Stokes system matrix and right hand // side, and the second is to create matrix // and right hand sides for the temperature - // dofs, which depends on the result for - // the velocity. + // dofs, which depends on the result of the + // linear system for the velocity. // - // This function does the first of these - // two tasks. There are two different - // situations for calling this - // function. The first one is when we reset - // the mesh, and both the matrix and the - // right hand side have to be - // generated. The second situation only - // sets up the right hand side. The reason - // for having two different accesses is - // that the matrix of the Stokes system - // does not change in time unless the mesh - // is changed, so we can save a - // considerable amount of work by doing the - // full assembly only when it is needed. + // This function is called at the beginning + // of each time step. In the first time step + // or if the mesh has changed, indicated by + // the rebuild_stokes_matrix, we + // need to assemble the Stokes matrix; on the + // other hand, if the mesh hasn't changed and + // the matrix is already available, this is + // not necessary and all we need to do is + // assemble the right hand side vector which + // changes in each time step. // // Regarding the technical details of - // implementation, not much has changed - // from step-22. We reset matrix and - // vector, create a quadrature formula on - // the cells and one on cell faces (for - // implementing Neumann boundary - // conditions). Then, we create a - // respective FEValues object for both the - // cell and the face integration. For the - // the update flags of the first, we - // perform the calculations of basis - // function derivatives only in case of a - // full assembly, since they are not needed - // otherwise, which makes the call of the - // FEValues::reinit function further down - // in the program more efficient. + // implementation, not much has changed from + // step-22. We reset matrix and vector, + // create a quadrature formula on the cells, + // and then create the respective FEValues + // object. For the update flags, we require + // basis function derivatives only in case of + // a full assembly, since they are not needed + // for the right hand side; as always, + // choosing the minimal set of flags + // depending on what is currently needed + // makes the call to FEValues::reinit further + // down in the program more efficient. // // There is one thing that needs to be - // commented – since we have a - // individual finite element and DoFHandler - // for the temperature, we need to generate - // a second FEValues object for the proper - // evaluation of the temperature - // solution. This isn't too complicated to - // realize here: just use the temperature - // structures and set an update flag for - // the basis function values which we need - // for evaluation of the temperature - // solution. + // commented – since we have a separate + // finite element and DoFHandler for the + // temperature, we need to generate a second + // FEValues object for the proper evaluation + // of the temperature solution. This isn't + // too complicated to realize here: just use + // the temperature structures and set an + // update flag for the basis function values + // which we need for evaluation of the + // temperature solution. The only important + // part to remember here is that the same + // quadrature formula is used for both + // FEValues objects to ensure that we get + // matching information when we loop over the + // quadrature points of the two objects. // // The declarations proceed with some // shortcuts for array sizes, the creation @@ -1592,20 +1598,19 @@ void BoussinesqFlowProblem::assemble_stokes_system () stokes_rhs=0; - QGauss quadrature_formula (stokes_degree+2); - - FEValues stokes_fe_values (stokes_fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - (rebuild_stokes_matrix == true - ? - update_gradients - : - UpdateFlags(0))); - - FEValues temperature_fe_values (temperature_fe, quadrature_formula, - update_values); + const QGauss quadrature_formula (stokes_degree+2); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + (rebuild_stokes_matrix == true + ? + update_gradients + : + UpdateFlags(0))); + + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values); const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -1615,21 +1620,24 @@ void BoussinesqFlowProblem::assemble_stokes_system () std::vector local_dof_indices (dofs_per_cell); - // The vector - // old_solution_values - // evaluates the temperature solution at - // the old time level at the quadrature - // points, which is needed for building - // the source term in the right hand side - // of the momentum equation. + // Next we need a vector that will contain + // the values of the temperature solution + // at the previous time level at the + // quadrature points to assemble the source + // term in the right hand side of the + // momentum equation. Let's call this vector + // old_solution_values. // // The set of vectors we create next hold // the evaluations of the basis functions - // that will be used for creating the - // matrices. This gives faster access to - // that data, which increases the - // performance of the assembly. See - // step-22 for details. + // as well as their gradients and + // symmetrized gradients that will be used + // for creating the matrices. Putting these + // into their own arrays rather than asking + // the FEValues object for this information + // each time it is needed is an + // optimization to accelerate the assembly + // process, see step-22 for details. // // The last two declarations are used to // extract the individual blocks @@ -1655,13 +1663,13 @@ void BoussinesqFlowProblem::assemble_stokes_system () // and the temperature system use the // same grid, but that's the only way to // keep degrees of freedom in sync. The - // first commands within the loop are + // first statements within the loop are // again all very familiar, doing the // update of the finite element data as // specified by the update flags, zeroing // out the local arrays and getting the // values of the old solution at the - // quadrature point. Then we are ready to + // quadrature points. Then we are ready to // loop over the quadrature points on the // cell. typename DoFHandler::active_cell_iterator @@ -1685,23 +1693,29 @@ void BoussinesqFlowProblem::assemble_stokes_system () { const double old_temperature = old_temperature_values[q]; - // Extract the basis relevant terms in - // the inner products once in advance as - // shown in step-22 in order to - // accelerate assembly. + // Next we extract the values and + // gradients of basis functions + // relevant to the terms in the + // inner products. As shown in + // step-22 this helps accelerate + // assembly. // - // Once this is done, we start the loop - // over the rows and columns of the local - // matrix and feed the matrix with the - // relevant products. The right hand side - // is filled with the forcing term driven - // by temperature in direction of gravity - // (which is vertical in our example). - // Note that the right hand side term is - // always generated, whereas the matrix - // contributions are only updated when it - // is requested by the - // rebuild_matrices flag. + // Once this is done, we start the + // loop over the rows and columns + // of the local matrix and feed the + // matrix with the relevant + // products. The right hand side is + // filled with the forcing term + // driven by temperature in + // direction of gravity (which is + // vertical in our example). Note + // that the right hand side term is + // always generated, whereas the + // matrix contributions are only + // updated when it is requested by + // the + // rebuild_matrices + // flag. for (unsigned int k=0; k::assemble_stokes_system () // The last step in the loop over all // cells is to enter the local // contributions into the global matrix - // and vector structures to the positions - // specified in + // and vector structures to the + // positions specified in // local_dof_indices. - // Again, we only add the matrix data - // when it is requested. Again, we let - // the ConstraintMatrix class do the - // insertion of the cell matrix elements - // to the global matrix, which already - // condenses the hanging node - // constraints. + // Again, we let the ConstraintMatrix + // class do the insertion of the cell + // matrix elements to the global + // matrix, which already condenses the + // hanging node constraints. cell->get_dof_indices (local_dof_indices); if (rebuild_stokes_matrix == true) @@ -1786,7 +1798,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () // the stiffness (diffusion) matrix. We // will then sum up the matrix plus the // stiffness matrix times the time step - // size. + // size once we know the actual time step. // // So the details for this first step are // very simple. In case we need to @@ -1806,7 +1818,7 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () temperature_mass_matrix = 0; temperature_stiffness_matrix = 0; - QGauss quadrature_formula(temperature_degree+2); + QGauss quadrature_formula (temperature_degree+2); FEValues temperature_fe_values (temperature_fe, quadrature_formula, update_values | update_gradients | update_JxW_values); @@ -1819,25 +1831,23 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () std::vector local_dof_indices (dofs_per_cell); - std::vector gamma_values (n_q_points); - std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); - // Now, let's start the loop over all - // cells in the triangulation. We need to - // zero out the local matrices, update - // the finite element evaluations, and - // then loop over the rows and columns of - // the matrices on each quadrature point, - // where we then create the mass matrix - // and the stiffness matrix (Laplace - // terms times the diffusion - // EquationData::kappa. Finally, - // we let the hanging node constraints - // insert these values into the global - // matrix, and directly condense the - // constraints into the matrix. + // Now, let's start the loop over all cells + // in the triangulation. We need to zero + // out the local matrices, update the + // finite element evaluations, and then + // loop over the rows and columns of the + // matrices on each quadrature point, where + // we then create the mass matrix and the + // stiffness matrix (Laplace terms times + // the diffusion + // EquationData::kappa. Finally, + // we let the constraints object insert + // these values into the global matrix, and + // directly condense the constraints into + // the matrix. typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); @@ -1904,15 +1914,15 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () // system, which means that we have to // evaluate second derivatives, specified // by the update flag - // update_hessians. The - // temperature equation is coupled to the + // update_hessians. + // + // The temperature equation is coupled to the // Stokes system by means of the fluid - // velocity, and these two parts of the - // solution are associated with different - // dof handlers. So we need to create a - // second FEValues object for the - // evaluation of the velocity at the - // quadrature points. + // velocity. These two parts of the solution + // are associated with different DoFHandlers, + // so we again need to create a second + // FEValues object for the evaluation of the + // velocity at the quadrature points. template void BoussinesqFlowProblem::assemble_temperature_system () { @@ -1933,13 +1943,15 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_rhs = 0; - QGauss quadrature_formula(temperature_degree+2); - FEValues temperature_fe_values (temperature_fe, quadrature_formula, - update_values | update_gradients | - update_hessians | - update_quadrature_points | update_JxW_values); - FEValues stokes_fe_values (stokes_fe, quadrature_formula, - update_values); + const QGauss quadrature_formula(temperature_degree+2); + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values | + update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_values); const unsigned int dofs_per_cell = temperature_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -1949,21 +1961,19 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector local_dof_indices (dofs_per_cell); - // Here comes the declaration of vectors - // to hold the old and present solution - // values and gradients for both the cell - // as well as faces to the cell, that - // will be generated from the global - // solution vectors. Next comes the - // declaration of an object to hold the - // temperature right hande side values, - // and we again use shortcuts for the - // temperature basis - // functions. Eventually, we need to find - // the maximum of velocity, temperature - // and the diameter of the computational - // domain which will be used for the - // definition of the stabilization + // Next comes the declaration of vectors to + // hold the old and present solution values + // and gradients at quadrature points of + // the current cell. We also declarate an + // object to hold the temperature right + // hande side values + // (gamma_values), and we + // again use shortcuts for the temperature + // basis functions. Eventually, we need to + // find the maximum of velocity, + // temperature and the diameter of the + // computational domain which will be used + // for the definition of the stabilization // parameter. std::vector > present_stokes_values (n_q_points, Vector(dim+1)); @@ -1985,16 +1995,21 @@ void BoussinesqFlowProblem::assemble_temperature_system () global_T_range = get_extrapolated_temperature_range(); const double global_Omega_diameter = GridTools::diameter (triangulation); - // Now, let's start the loop over all - // cells in the triangulation. First set - // the local rhs to zero, and then get - // the values of the old solution - // functions (and the current velocity) - // at the quadrature points, since they - // are going to be needed for the - // definition of the stabilization - // parameters and as coefficients in the - // equation, respectively. + // Now, let's start the loop over all cells + // in the triangulation. Again, we need two + // cell iterators that walk in parallel + // through the cells of the two involved + // DoFHandler objects for the Stokes and + // temperature part. Within the loop, we + // first set the local rhs to zero, and + // then get the values and derivatives of + // the old solution functions (and the + // current velocity) at the quadrature + // points, since they are going to be + // needed for the definition of the + // stabilization parameters and as + // coefficients in the equation, + // respectively. typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); @@ -2035,7 +2050,7 @@ void BoussinesqFlowProblem::assemble_temperature_system () // discussion in the introduction // using the dedicated // function. With that at hand, we - // can define get into the loop + // can get into the loop // over quadrature points and local // rhs vector components. The terms // here are quite lenghty, but @@ -2146,19 +2161,19 @@ void BoussinesqFlowProblem::assemble_temperature_system () // @sect4{BoussinesqFlowProblem::solve} // // This function solves the linear - // equation systems. According to + // systems of equations. Following to // the introduction, we start with // the Stokes system, where we need // to generate our block Schur // preconditioner. Since all the // relevant actions are implemented // in the class - // BlockSchurPreconditioner, + // BlockSchurPreconditioner, // all we have to do is to // initialize the class // appropriately. What we need to // pass down is an - // InverseMatrix object + // InverseMatrix object // for the pressure mass matrix, // which we set up using the // respective class together with @@ -2167,10 +2182,10 @@ void BoussinesqFlowProblem::assemble_temperature_system () // preconditioner for the // velocity-velocity matrix. Note // that both - // Mp_preconditioner and - // Amg_preconditioner are + // Mp_preconditioner and + // Amg_preconditioner are // only pointers, so we use - // * to pass down the + // * to pass down the // actual preconditioner objects. // // Once the preconditioner is @@ -2182,12 +2197,12 @@ void BoussinesqFlowProblem::assemble_temperature_system () // the solver. GMRES needs to // internally store temporary // vectors for each iteration (see - // even the discussion in the + // the discussion in the // results section of step-22) // – the more vectors it can // use, the better it will - // generally perform. To let memory - // demands not increase to much, we + // generally perform. To keep memory + // demands in check, we // set the number of vectors to // 100. This means that up to 100 // solver iterations, every @@ -2197,16 +2212,15 @@ void BoussinesqFlowProblem::assemble_temperature_system () // specified tolerance, it will // work on a reduced set of vectors // by restarting at every 100 - // iterations. Then, we solve the - // system and distribute the - // constraints in the Stokes - // system, i.e. hanging nodes and - // no-flux boundary condition, in - // order to have the appropriate - // solution values even at - // constrained dofs. Finally, we - // write the number of iterations - // to the screen. + // iterations. + // + // With this all set up, we solve the system + // and distribute the constraints in the + // Stokes system, i.e. hanging nodes and + // no-flux boundary condition, in order to + // have the appropriate solution values even + // at constrained dofs. Finally, we write the + // number of iterations to the screen. template void BoussinesqFlowProblem::solve () { @@ -2237,39 +2251,47 @@ void BoussinesqFlowProblem::solve () << std::endl; } - // Once we know the Stokes - // solution, we can determine the - // new time step from the maximal - // velocity. We have to do this to + // Once we know the Stokes solution, we can + // determine the new time step from the + // maximal velocity. We have to do this to // satisfy the CFL condition since - // convection terms are treated - // explicitly in the temperature - // equation, as discussed in the - // introduction. Next we set up the - // temperature system and the right - // hand side using the function - // assemble_temperature_system(). Knowing - // the matrix and right hand side - // of the temperature equation, we - // set up a preconditioner and a - // solver. The temperature matrix - // is a mass matrix plus a Laplace - // matrix times a small number, the - // time step. Hence, the mass - // matrix dominates and we get a - // reasonable good preconditioner - // by simple means, namely SSOR. We - // set the relaxation parameter to - // 1.2. As a solver, we choose the - // conjugate gradient method CG. As - // before, we tell the solver to - // use Trilinos vectors via the - // template argument - // TrilinosWrappers::Vector - // at construction. Finally, we - // solve, distribute the hanging - // node constraints and write out - // the number of iterations. + // convection terms are treated explicitly + // in the temperature equation, as + // discussed in the introduction. The exact + // form of the formula used here for the + // time step is discussed in the results + // section of this program. + // + // Next we set up the temperature system + // and the right hand side using the + // function + // assemble_temperature_system(). Knowing + // the matrix and right hand side of the + // temperature equation, we set up a + // preconditioner and a solver. The + // temperature matrix is a mass matrix + // (with eigenvalues around one) plus a + // Laplace matrix (with eigenvalues between + // zero and $ch^{-2}$) times a small number + // proportional to the time step + // $k_n$. Hence, the resulting symmetric + // and positive definite matrix has + // eigenvalues in the range + // $[1,1+k_nh^{-2}]$ (up to + // constants). This matrix is only + // moderately ill conditioned even for + // small mesh sizes and we get a reasonable + // good preconditioner by simple means, for + // example SSOR. We set the relaxation + // parameter to 1.2. As a solver, we choose + // the conjugate gradient method CG. As + // before, we tell the solver to use + // Trilinos vectors via the template + // argument + // TrilinosWrappers::Vector at + // construction. Finally, we solve, + // distribute the hanging node constraints + // and write out the number of iterations. old_time_step = time_step; time_step = 1./(1.6*dim*std::sqrt(1.*dim)) / temperature_degree * @@ -2299,11 +2321,14 @@ void BoussinesqFlowProblem::solve () << " CG iterations for temperature." << std::endl; - // In the end of this function, we - // step through the vector and read - // out the maximum and minimum - // temperature value, which we also - // want to output. + // At the end of this function, we step + // through the vector and read out the + // maximum and minimum temperature value, + // which we also want to output. This + // will come in handy when determining + // the correct constant in the choice of + // time step as discuss in the results + // section of this program. double min_temperature = temperature_solution(0), max_temperature = temperature_solution(0); for (unsigned int i=0; i::solve () // @sect4{BoussinesqFlowProblem::output_results} // - // This function writes the - // solution to a vtk output file - // for visualization, which is done - // every tenth time step. This is - // usually a quite simple task, - // since the deal.II library - // provides functions that do - // almost all the job for us. In - // this case, the situation is a - // bit more complicated, since we - // want to visualize both the - // Stokes solution and the - // temperature as one data set, but - // we have done all the - // calculations based on two - // different. The way we're going - // to achieve this recombination is - // to create a joint DoFHandler - // that collects both components, - // the Stokes solution and the - // temperature solution. This can - // be nicely done by combining the - // finite elements from the two - // systems to form one FESystem, - // and let this collective system - // define a new DoFHandler - // object. To be sure that - // everything was done correctly, - // we perform a sanity check that - // ensures that we got all the dofs - // from both Stokes and temperature - // even in the combined system. + // This function writes the solution to a VTK + // output file for visualization, which is + // done every tenth time step. This is + // usually quite a simple task, since the + // deal.II library provides functions that do + // almost all the job for us. In this case, + // the situation is a bit more complicated, + // since we want to visualize both the Stokes + // solution and the temperature as one data + // set, but we have done all the calculations + // based on two different DoFHandler objects, + // a situation the DataOut class usually used + // for output is not prepared to deal + // with. The way we're going to achieve this + // recombination is to create a joint + // DoFHandler that collects both components, + // the Stokes solution and the temperature + // solution. This can be nicely done by + // combining the finite elements from the two + // systems to form one FESystem, and let this + // collective system define a new DoFHandler + // object. To be sure that everything was + // done correctly, we perform a sanity check + // that ensures that we got all the dofs from + // both Stokes and temperature even in the + // combined system. // - // Next, we create a vector that - // collects the actual solution - // values (up to now, we've just - // provided the tools for it - // without reading any data. Since - // this vector is only going to be - // used for output, we create it as - // a deal.II vector that nicely - // cooperate with the data output - // classes. Remember that we used - // Trilinos vectors for assembly - // and solving. + // Next, we create a vector that will collect + // the actual solution values. Since this + // vector is only going to be used for + // output, we create it as a deal.II vector + // that nicely cooperate with the data output + // classes. Remember that we used Trilinos + // vectors for assembly and solving. template void BoussinesqFlowProblem::output_results () const { @@ -2386,41 +2400,33 @@ void BoussinesqFlowProblem::output_results () const Vector joint_solution (joint_dof_handler.n_dofs()); // Unfortunately, there is no - // straight-forward relation that - // tells us how to sort Stokes and - // temperature vector into the - // joint vector. The way we can get - // around this trouble is to rely - // on the information collected in - // the FESystem. For each dof in a - // cell, the joint finite element - // knows to which equation - // component (velocity component, - // pressure, or temperature) it - // belongs – that's the - // information we need! So we step - // through all cells (as a - // complication, we need to create - // iterations for the cells in the - // Stokes system and the - // temperature system, too, even - // though they are the same in all - // the three cases), and for each - // joint cell dof, we read out that - // component using the function - // joint_fe.system_to_base_index(i).second. We - // also need to keep track whether - // we're on a Stokes dof or a - // temperature dof, which is - // contained in - // joint_fe.system_to_base_index(i).first.first. Eventually, - // the dof_indices data structures - // on either of the three systems - // tell us how the relation between - // global vector and local dofs - // looks like on the present cell, - // which concludes this tedious - // work. + // straight-forward relation that tells us + // how to sort Stokes and temperature + // vector into the joint vector. The way we + // can get around this trouble is to rely + // on the information collected in the + // FESystem. For each dof in a cell, the + // joint finite element knows to which + // equation component (velocity component, + // pressure, or temperature) it belongs + // – that's the information we need! + // So we step through all cells (with + // iterators into all three DoFHandlers + // moving in synch), and for each joint + // cell dof, we read out that component + // using the + // FiniteElement::system_to_base_index + // function (see there for a description of + // what the various parts of its return + // value contain). We also need to keep + // track whether we're on a Stokes dof or a + // temperature dof, which is contained in + // joint_fe.system_to_base_index(i).first.first. Eventually, + // the dof_indices data structures on + // either of the three systems tell us how + // the relation between global vector and + // local dofs looks like on the present + // cell, which concludes this tedious work. { std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); std::vector local_stokes_dof_indices (stokes_fe.dofs_per_cell); @@ -2468,7 +2474,7 @@ void BoussinesqFlowProblem::output_results () const // the individual components), and // attach the joint dof handler to // a DataOut object. The first - // dim components are the + // dim components are the // vector velocity, and then we // have pressure and // temperature. This information is @@ -2517,42 +2523,57 @@ void BoussinesqFlowProblem::output_results () const // @sect4{BoussinesqFlowProblem::refine_mesh} // - // This function takes care of the - // adaptive mesh refinement. The - // three tasks this function - // performs is to first find out - // which cells to refine/coarsen, - // then to actually do the - // refinement and eventually - // transfer the solution vectors - // between the two different - // grids. The first task is simply - // achieved by using the - // well-established Kelly error - // estimator on the temperature (it - // is the temperature we're mainly - // interested in for this program, - // and we need to be accurate in - // regions of high temperature - // gradients, also to not have too - // much numerical diffusion). The - // second task is to actually do - // the remeshing. That involves - // only basic functions as well, - // such as the - // refine_and_coarsen_fixed_fraction - // that refines the 80 precent of - // the cells which have the largest - // estimated error and coarsens the - // 10 precent with the smallest - // error. For reasons of limited - // computer ressources, we have to - // set a limit on the maximum - // refinement level. We do this - // after the refinement indicator - // has been applied to the cells, - // and simply unselect cells with - // too high grid level. + // This function takes care of the adaptive + // mesh refinement. The three tasks this + // function performs is to first find out + // which cells to refine/coarsen, then to + // actually do the refinement and eventually + // transfer the solution vectors between the + // two different grids. The first task is + // simply achieved by using the + // well-established Kelly error estimator on + // the temperature (it is the temperature + // we're mainly interested in for this + // program, and we need to be accurate in + // regions of high temperature gradients, + // also to not have too much numerical + // diffusion). The second task is to actually + // do the remeshing. That involves only basic + // functions as well, such as the + // refine_and_coarsen_fixed_fraction + // that refines those cells with the largest + // estimated error that together make up 80 + // per cent of the error, and coarsens those + // cells with the smallest error that make up + // for a combined 10 per cent of the + // error. + // + // If implemented like this, we would get a + // program that will not make much progress: + // Remember that we expect temperature fields + // that are nearly discontinuous (the + // diffusivity $\kappa$ is very small after + // all) and consequently we can expect that a + // freely adapted mesh will refine further + // and further into the areas of large + // gradients. This decrease in mesh size will + // then be accompanied by a decrease in time + // step, requiring an exceedingly large + // number of time steps to solve to a given + // final time. It will also lead to meshes + // that are much better at resolving + // discontinuities after several mesh + // refinement cycles than in the beginning. + // + // In particular to prevent the decrease in + // time step size and the correspondingly + // large number of time steps, we limit the + // maximal refinement depth of the mesh. To + // this end, after the refinement indicator + // has been applied to the cells, we simply + // loop over all cells on the finest level + // and unselect them from refinement if they + // would result in too high a mesh level. template void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) { @@ -2614,8 +2635,10 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) TrilinosWrappers::BlockVector x_stokes(2); x_stokes = stokes_solution; - SolutionTransfer temperature_trans(temperature_dof_handler); - SolutionTransfer stokes_trans(stokes_dof_handler); + SolutionTransfer + temperature_trans(temperature_dof_handler); + SolutionTransfer + stokes_trans(stokes_dof_handler); triangulation.prepare_coarsening_and_refinement(); temperature_trans.prepare_for_coarsening_and_refinement(x_temperature); @@ -2626,7 +2649,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // dof structure on the new grid, // and initialize the matrix // structures and the new vectors - // in the setup_dofs + // in the setup_dofs // function. Next, we actually // perform the interpolation of the // solutions between the grids. We @@ -2677,7 +2700,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // refinement and additional // adative refinement steps, and // then create a cube in - // dim dimensions and set + // dim dimensions and set // up the dofs for the first // time. Since we want to start the // time stepping already with an @@ -2691,7 +2714,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // Before we start, we project the // initial values to the grid and // obtain the first data for the - // old_temperature_solution + // old_temperature_solution // vector. Then, we initialize time // step number and time step and // start the time loop. diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 3d60697308..73415567ba 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -723,11 +723,11 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () { stokes_preconditioner_matrix = 0; - QGauss quadrature_formula(stokes_degree+2); - FEValues stokes_fe_values (stokes_fe, quadrature_formula, - update_JxW_values | - update_values | - update_gradients); + const QGauss quadrature_formula (stokes_degree+2); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_JxW_values | + update_values | + update_gradients); const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -826,34 +826,29 @@ void BoussinesqFlowProblem::assemble_stokes_system () stokes_rhs=0; - QGauss quadrature_formula(stokes_degree+2); - QGauss face_quadrature_formula(stokes_degree+2); - - FEValues stokes_fe_values (stokes_fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - (rebuild_stokes_matrix == true - ? - update_gradients - : - UpdateFlags(0))); - - FEValues temperature_fe_values (temperature_fe, quadrature_formula, - update_values); + const QGauss quadrature_formula(stokes_degree+2); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + (rebuild_stokes_matrix == true + ? + update_gradients + : + UpdateFlags(0))); + + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values); const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_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); - std::vector boundary_values (n_face_q_points); - std::vector old_temperature_values(n_q_points); std::vector > phi_u (dofs_per_cell); @@ -950,10 +945,10 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () temperature_mass_matrix = 0; temperature_stiffness_matrix = 0; - QGauss quadrature_formula(temperature_degree+2); - FEValues temperature_fe_values (temperature_fe, quadrature_formula, - update_values | update_gradients | - update_JxW_values); + const QGauss quadrature_formula(temperature_degree+2); + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values | update_gradients | + update_JxW_values); const unsigned int dofs_per_cell = temperature_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -963,8 +958,6 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () std::vector local_dof_indices (dofs_per_cell); - std::vector gamma_values (n_q_points); - std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); @@ -1041,11 +1034,13 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_rhs = 0; - QGauss quadrature_formula(temperature_degree+2); - FEValues temperature_fe_values (temperature_fe, quadrature_formula, - update_values | update_gradients | - update_hessians | - update_quadrature_points | update_JxW_values); + const QGauss quadrature_formula(temperature_degree+2); + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values | + update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values); FEValues stokes_fe_values (stokes_fe, quadrature_formula, update_values); -- 2.39.5