From: bangerth Date: Wed, 18 Dec 2013 20:18:21 +0000 (+0000) Subject: Finish documenting step-26. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d17a76088f1800458a0864ccd37d74686e08f7db;p=dealii-svn.git Finish documenting step-26. git-svn-id: https://svn.dealii.org/trunk@32053 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-26/doc/builds-on b/deal.II/examples/step-26/doc/builds-on index 48a0f73876..17402734c7 100644 --- a/deal.II/examples/step-26/doc/builds-on +++ b/deal.II/examples/step-26/doc/builds-on @@ -1 +1 @@ -step-4 +step-6 diff --git a/deal.II/examples/step-26/doc/kind b/deal.II/examples/step-26/doc/kind index c1d9154931..86a44aa1ef 100644 --- a/deal.II/examples/step-26/doc/kind +++ b/deal.II/examples/step-26/doc/kind @@ -1 +1 @@ -techniques +time dependent diff --git a/deal.II/examples/step-26/doc/results.dox b/deal.II/examples/step-26/doc/results.dox index 6a401f850d..798cc8f298 100644 --- a/deal.II/examples/step-26/doc/results.dox +++ b/deal.II/examples/step-26/doc/results.dox @@ -38,7 +38,7 @@ Number of active cells: 1803 Number of degrees of freedom: 2109 @endcode -Maybe of more interest is a visualization of the solution and the mesh on wich +Maybe of more interest is a visualization of the solution and the mesh on which it was computed: HeatEquation implementation} // - // The next step then is the implementation of the main class. + // It is time now for the implementation of the main class. Let's + // start with the constructor which selects a linear element, a time + // step constant at 1/500 (remember that one period of the source + // on the right hand side was set to 0.2 above, so we resolve each + // period with 100 time steps) and chooses the Crank Nicolson method + // by setting $\theta=1/2$. template HeatEquation::HeatEquation () : @@ -202,11 +207,30 @@ namespace Step26 dof_handler(triangulation), time_step(1. / 500), theta(0.5) - { - } + {} + // @sect4{HeatEquation::setup_system} + // + // The next function is the one that sets up the DoFHandler object, + // computes the constraints, and sets the linear algebra objects + // to their correct sizes. We also compute the mass and Laplace + // matrix here by simply calling two functions in the library. + // + // Note that we compute these matrices taking into account already the + // constraints due to hanging nodes. These are all homogenous, i.e., + // they only consist of constraints of the form $U_i = \alpha_{ij} U_j + // + \alpha_{ik} U_k$ (whereas inhomogenous constraints would also + // have a term not proportional to $U$, i.e., $U_i = \alpha_{ij} U_j + // + \alpha_{ik} U_k + c_i$). For this kind of constraint, we can + // eliminate hanging nodes independently in the matrix and the + // right hand side vectors, but this is not the case for inhomogenous + // constraints for which we can eliminate constrained degrees of freedom + // only by looking at both the system matrix and corresponding right + // right hand side at the same time. This may become a problem when + // dealing with non-zero Dirichlet boundary conditions, though we + // do not do this here in the current program. template void HeatEquation::setup_system() { @@ -254,7 +278,10 @@ namespace Step26 } - + // @sect4{HeatEquation::solve_time_step} + // + // The next function is the one that solves the actual linear system + // for a single time step. There is nothing surprising here: template void HeatEquation::solve_time_step() { @@ -264,7 +291,8 @@ namespace Step26 PreconditionSSOR<> preconditioner; preconditioner.initialize(system_matrix, 1.0); - cg.solve(system_matrix, solution, system_rhs, preconditioner); + cg.solve(system_matrix, solution, system_rhs, + preconditioner); constraints.distribute(solution); @@ -274,6 +302,9 @@ namespace Step26 + // @sect4{HeatEquation::output_results} + // + // Neither is there anything new in generating graphical output: template void HeatEquation::output_results() const { @@ -292,40 +323,29 @@ namespace Step26 } - // @sect4{BoussinesqFlowProblem::refine_mesh} + // @sect4{HeatEquation::refine_mesh} // - // This function takes care of the adaptive mesh refinement. The three tasks + // This function is the interesting part of the program. It 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 + // estimator on the solution. 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 + // with the largest estimated error that together make up 60 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. + // a combined 40 per cent of the error. Note that for problems such as the + // current one where the areas where something is going on are shifting + // around, we want to aggressively coarsen so that we can move cells + // around to where it is necessary. // - // 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. + // As already discussed in the introduction, too small a mesh leads to + // too small a time step, whereas too large a mesh leads to too little + // resolution. Consequently, after the first two steps, we have two + // loops that limit refinement and coarsening to an allowable range of + // cells: template void HeatEquation::refine_mesh (const unsigned int min_grid_level, const unsigned int max_grid_level) @@ -341,6 +361,7 @@ namespace Step26 GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, estimated_error_per_cell, 0.6, 0.4); + if (triangulation.n_levels() > max_grid_level) for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(max_grid_level); @@ -357,53 +378,48 @@ namespace Step26 // SolutionTransfer class and we have to prepare the solution vectors that // should be transferred to the new grid (we will lose the old grid once // we have done the refinement so the transfer has to happen concurrently - // with refinement). What we definitely need are the current and the old - // temperature (BDF-2 time stepping requires two old solutions). Since the - // SolutionTransfer objects only support to transfer one object per dof - // handler, we need to collect the two temperature solutions in one data - // structure. Moreover, we choose to transfer the Stokes solution, too, - // since we need the velocity at two previous time steps, of which only - // one is calculated on the fly. + // with refinement). At the point where we call this function, we will + // have just computed the solution, so we no longer need the old_solution + // variable (it will be overwritten by the solution just after the mesh + // may have been refined, i.e., at the end of the time step; see below). + // In other words, we only need the one solution vector, and we copy it + // to a temporary object where it is safe from being reset when we further + // down below call setup_system(). // - // Consequently, we initialize two SolutionTransfer objects for the Stokes - // and temperature DoFHandler objects, by attaching them to the old dof - // handlers. With this at place, we can prepare the triangulation and the - // data vectors for refinement (in this order). - std::vector > x_solution (2); - x_solution[0] = solution; - x_solution[1] = old_solution; - + // Consequently, we initialize a SolutionTransfer object by attaching + // it to the old DoF handler. We then prepare the triangulation and the + // data vector for refinement (in this order). SolutionTransfer solution_trans(dof_handler); + Vector previous_solution; + previous_solution = solution; triangulation.prepare_coarsening_and_refinement(); - solution_trans.prepare_for_coarsening_and_refinement(x_solution); + solution_trans.prepare_for_coarsening_and_refinement(previous_solution); // Now everything is ready, so do the refinement and recreate the dof // structure on the new grid, and initialize the matrix structures and the - // new vectors in the setup_dofs function. Next, we actually - // perform the interpolation of the solutions between the grids. We create - // another copy of temporary vectors for temperature (now corresponding to - // the new grid), and let the interpolate function do the job. Then, the - // resulting array of vectors is written into the respective vector member - // variables. For the Stokes vector, everything is just the same – - // except that we do not need another temporary vector since we just - // interpolate a single vector. In the end, we have to tell the program - // that the matrices and preconditioners need to be regenerated, since the - // mesh has changed. + // new vectors in the setup_system function. Next, we actually + // perform the interpolation of the solution from old to new grid. triangulation.execute_coarsening_and_refinement (); setup_system (); - std::vector > tmp (2); - tmp[0].reinit (solution); - tmp[1].reinit (solution); - solution_trans.interpolate(x_solution, tmp); - - solution = tmp[0]; - old_solution = tmp[1]; + solution_trans.interpolate(previous_solution, solution); } + // @sect4{HeatEquation::run} + // + // This is the main driver of the program, where we loop over all + // time steps. At the top of the function, we set the number of + // initial global mesh refinements and the number of initial cycles of + // adaptive mesh refinement by repeating the first time step a few + // times. Then we create a mesh, initialize the various objects we will + // work with, set a label for where we should start when re-running + // the first time step, and interpolate the initial solution onto + // out mesh (we choose the zero function here, which of course we could + // do in a simpler way by just setting the solution vector to zero). We + // also output the initial time step once. template void HeatEquation::run() { @@ -422,6 +438,10 @@ namespace Step26 start_time_iteration: + tmp.reinit (solution.size()); + forcing_terms.reinit (solution.size()); + + VectorTools::interpolate(dof_handler, ZeroFunction(), old_solution); @@ -432,6 +452,12 @@ start_time_iteration: output_results(); + // Then we start the main loop until the computed time exceeds our + // end time of 0.5. The first task is to build the right hand + // side of the linear system we need to solve in each time step. + // Recall that it contains the term $MU^{n-1}-(1-\theta)k_n AU^{n-1}$. + // We put these terms into the variable system_rhs, with the + // help of a temporary vector: while (time <= 0.5) { time += time_step; @@ -440,14 +466,18 @@ start_time_iteration: std::cout << "Time step " << timestep_number << " at t=" << time << std::endl; - tmp.reinit (solution.size()); - forcing_terms.reinit (solution.size()); - mass_matrix.vmult(system_rhs, old_solution); laplace_matrix.vmult(tmp, old_solution); system_rhs.add(-(1 - theta) * time_step, tmp); + // The second piece is to compute the contributions of the source + // terms. This corresponds to the term $k_n + // \left[ (1-\theta)F^{n-1} + \theta F^n \right]$. The following + // code calls VectorTools::create_right_hand_side to compute the + // vectors $F$, where we set the time of the right hand side + // (source) function before we evaluate it. The result of this + // all ends up in the forcing_terms variable: RightHandSide rhs_function; rhs_function.set_time(time); VectorTools::create_right_hand_side(dof_handler, @@ -465,8 +495,25 @@ start_time_iteration: forcing_terms.add(time_step * (1 - theta), tmp); + // Next, we add the forcing terms to the ones that + // come from the time stepping, and also build the matrix + // $M+k_n\theta A$ that we have to invert in each time step. + // The final piece of these operations is to eliminate + // hanging node constrained degrees of freedom from the + // linear system: system_rhs += forcing_terms; + system_matrix.copy_from(mass_matrix); + system_matrix.add(theta * time_step, laplace_matrix); + + constraints.condense (system_matrix, system_rhs); + + // There is one more operation we need to do before we + // can solve it: boundary values. To this end, we create + // a boundary value object, set the proper time to the one + // of the current time step, and evaluate it as we have + // done many times before. The result is used to also + // set the correct boundary values in the linear system: { BoundaryValues boundary_values_function; boundary_values_function.set_time(time); @@ -477,20 +524,27 @@ start_time_iteration: boundary_values_function, boundary_values); - system_matrix.copy_from(mass_matrix); - system_matrix.add(theta * time_step, laplace_matrix); MatrixTools::apply_boundary_values(boundary_values, system_matrix, solution, system_rhs); } - constraints.condense (system_matrix, system_rhs); - + // With this out of the way, all we have to do is solve the + // system, generate graphical data, and... solve_time_step(); output_results(); + // ...take care of mesh refinement. Here, what we want to do is + // (i) refine the requested number of times at the very beginning + // of the solution procedure, after which we jump to the top to + // restart the time iteration, (ii) refine every fifth time + // step after that. + // + // The time loop and, indeed, the main part of the program ends + // with starting into the next time step by setting old_solution + // to the solution we have just computed. if ((timestep_number == 1) && (pre_refinement_step < n_adaptive_pre_refinement_steps)) { @@ -498,19 +552,32 @@ start_time_iteration: initial_global_refinement + n_adaptive_pre_refinement_steps); ++pre_refinement_step; + tmp.reinit (solution.size()); + forcing_terms.reinit (solution.size()); + std::cout << std::endl; goto start_time_iteration; } else if ((timestep_number > 0) && (timestep_number % 5 == 0)) - refine_mesh (initial_global_refinement, - initial_global_refinement + n_adaptive_pre_refinement_steps); + { + refine_mesh (initial_global_refinement, + initial_global_refinement + n_adaptive_pre_refinement_steps); + tmp.reinit (solution.size()); + forcing_terms.reinit (solution.size()); + } old_solution = solution; } } } + +// @sect3{The main function} +// +// Having made it this far, there is, again, nothing +// much to discuss for the main function of this +// program: it looks like all such functions since step-6. int main() { try