From c4d08de39e8feb56fbd83318199ebeaa27e07674 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 30 Jun 2023 11:07:17 -0600 Subject: [PATCH] Create the framework for step-86. --- examples/step-86/CMakeLists.txt | 39 ++ examples/step-86/doc/builds-on | 1 + examples/step-86/doc/intro.dox | 2 + examples/step-86/doc/kind | 1 + examples/step-86/doc/results.dox | 1 + examples/step-86/doc/tooltip | 1 + examples/step-86/step-86.cc | 697 +++++++++++++++++++++++++++++++ 7 files changed, 742 insertions(+) create mode 100644 examples/step-86/CMakeLists.txt create mode 100644 examples/step-86/doc/builds-on create mode 100644 examples/step-86/doc/intro.dox create mode 100644 examples/step-86/doc/kind create mode 100644 examples/step-86/doc/results.dox create mode 100644 examples/step-86/doc/tooltip create mode 100644 examples/step-86/step-86.cc diff --git a/examples/step-86/CMakeLists.txt b/examples/step-86/CMakeLists.txt new file mode 100644 index 0000000000..209c9a0d5e --- /dev/null +++ b/examples/step-86/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-86 tutorial program: +## + +# Set the name of the project and target: +set(TARGET "step-86") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# file(GLOB_RECURSE TARGET_SRC "source/*.cc") +# file(GLOB_RECURSE TARGET_INC "include/*.h") +# set(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +set(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +cmake_minimum_required(VERSION 3.13.4) + +find_package(deal.II 9.6.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +endif() + +deal_ii_initialize_cached_variables() +project(${TARGET}) +deal_ii_invoke_autopilot() diff --git a/examples/step-86/doc/builds-on b/examples/step-86/doc/builds-on new file mode 100644 index 0000000000..1aabbdfdd6 --- /dev/null +++ b/examples/step-86/doc/builds-on @@ -0,0 +1 @@ +step-26 diff --git a/examples/step-86/doc/intro.dox b/examples/step-86/doc/intro.dox new file mode 100644 index 0000000000..c59cfdf609 --- /dev/null +++ b/examples/step-86/doc/intro.dox @@ -0,0 +1,2 @@ + +

Introduction

diff --git a/examples/step-86/doc/kind b/examples/step-86/doc/kind new file mode 100644 index 0000000000..86a44aa1ef --- /dev/null +++ b/examples/step-86/doc/kind @@ -0,0 +1 @@ +time dependent diff --git a/examples/step-86/doc/results.dox b/examples/step-86/doc/results.dox new file mode 100644 index 0000000000..f4c6feefb5 --- /dev/null +++ b/examples/step-86/doc/results.dox @@ -0,0 +1 @@ +

Results

diff --git a/examples/step-86/doc/tooltip b/examples/step-86/doc/tooltip new file mode 100644 index 0000000000..c5653ac287 --- /dev/null +++ b/examples/step-86/doc/tooltip @@ -0,0 +1 @@ +Using better time steppers for the heat equation. diff --git a/examples/step-86/step-86.cc b/examples/step-86/step-86.cc new file mode 100644 index 0000000000..058889812a --- /dev/null +++ b/examples/step-86/step-86.cc @@ -0,0 +1,697 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2013 - 2021 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Author: Wolfgang Bangerth, Texas A&M University, 2023 + */ + + +// The program starts with the usual include files, all of which you should +// have seen before by now: +#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 the usual placing of all content of this program into a namespace and +// the importation of the deal.II namespace into the one we will work in: +namespace Step86 +{ + using namespace dealii; + + + // @sect3{The HeatEquation class} + // + // The next piece is the declaration of the main class of this program. It + // follows the well trodden path of previous examples. If you have looked at + // step-6, for example, the only thing worth noting here is that we need to + // build two matrices (the mass and Laplace matrix) and keep the current and + // previous time step's solution. We then also need to store the current + // time, the size of the time step, and the number of the current time + // step. The last of the member variables denotes the theta parameter + // discussed in the introduction that allows us to treat the explicit and + // implicit Euler methods as well as the Crank-Nicolson method and other + // generalizations all in one program. + // + // As far as member functions are concerned, the only possible surprise is + // that the refine_mesh function takes arguments for the + // minimal and maximal mesh refinement level. The purpose of this is + // discussed in the introduction. + template + class HeatEquation + { + public: + HeatEquation(); + void run(); + + private: + void setup_system(); + void solve_time_step(); + void output_results() const; + void refine_mesh(const unsigned int min_grid_level, + const unsigned int max_grid_level); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix mass_matrix; + SparseMatrix laplace_matrix; + SparseMatrix system_matrix; + + Vector solution; + Vector old_solution; + Vector system_rhs; + + double time; + double time_step; + unsigned int timestep_number; + + const double theta; + }; + + + + // @sect3{Equation data} + + // In the following classes and functions, we implement the various pieces + // of data that define this problem (right hand side and boundary values) + // that are used in this program and for which we need function objects. The + // right hand side is chosen as discussed at the end of the + // introduction. For boundary values, we choose zero values, but this is + // easily changed below. + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function() + , period(0.2) + {} + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + const double period; + }; + + + + template + double RightHandSide::value(const Point & p, + const unsigned int component) const + { + (void)component; + AssertIndexRange(component, 1); + Assert(dim == 2, ExcNotImplemented()); + + const double time = this->get_time(); + const double point_within_period = + (time / period - std::floor(time / period)); + + if ((point_within_period >= 0.0) && (point_within_period <= 0.2)) + { + if ((p[0] > 0.5) && (p[1] > -0.5)) + return 1; + else + return 0; + } + else if ((point_within_period >= 0.5) && (point_within_period <= 0.7)) + { + if ((p[0] > -0.5) && (p[1] > 0.5)) + return 1; + else + return 0; + } + else + return 0; + } + + + + template + class BoundaryValues : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + + + template + double BoundaryValues::value(const Point & /*p*/, + const unsigned int component) const + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + return 0; + } + + + + // @sect3{The HeatEquation implementation} + // + // 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() + : fe(1) + , 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 do not take the hanging node constraints into account when + // assembling the matrices (both functions have an AffineConstraints argument + // that defaults to an empty object). This is because we are going to + // condense the constraints in run() after combining the matrices for the + // current time-step. + template + void HeatEquation::setup_system() + { + dof_handler.distribute_dofs(fe); + + std::cout << std::endl + << "===========================================" << std::endl + << "Number of active cells: " << triangulation.n_active_cells() + << std::endl + << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << std::endl; + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /*keep_constrained_dofs = */ true); + sparsity_pattern.copy_from(dsp); + + mass_matrix.reinit(sparsity_pattern); + laplace_matrix.reinit(sparsity_pattern); + system_matrix.reinit(sparsity_pattern); + + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree + 1), + mass_matrix); + MatrixCreator::create_laplace_matrix(dof_handler, + QGauss(fe.degree + 1), + laplace_matrix); + + solution.reinit(dof_handler.n_dofs()); + old_solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + + // @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() + { + SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); + SolverCG> cg(solver_control); + + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.0); + + cg.solve(system_matrix, solution, system_rhs, preconditioner); + + constraints.distribute(solution); + + std::cout << " " << solver_control.last_step() << " CG iterations." + << std::endl; + } + + + + // @sect4{HeatEquation::output_results} + // + // Neither is there anything new in generating graphical output other than the + // fact that we tell the DataOut object what the current time and time step + // number is, so that this can be written into the output file: + template + void HeatEquation::output_results() const + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "U"); + + data_out.build_patches(); + + data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number)); + + const std::string filename = + "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk"; + std::ofstream output(filename); + data_out.write_vtk(output); + } + + + // @sect4{HeatEquation::refine_mesh} + // + // 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 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 60 per cent of the + // error, and coarsens those cells with the smallest error that make up for + // 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. + // + // 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) + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, + estimated_error_per_cell, + 0.6, + 0.4); + + if (triangulation.n_levels() > max_grid_level) + for (const auto &cell : + triangulation.active_cell_iterators_on_level(max_grid_level)) + cell->clear_refine_flag(); + for (const auto &cell : + triangulation.active_cell_iterators_on_level(min_grid_level)) + cell->clear_coarsen_flag(); + // These two loops above are slightly different but this is easily + // explained. In the first loop, instead of calling + // triangulation.end() we may as well have called + // triangulation.end_active(max_grid_level). The two + // calls should yield the same iterator since iterators are sorted + // by level and there should not be any cells on levels higher than + // on level max_grid_level. In fact, this very piece + // of code makes sure that this is the case. + + // As part of mesh refinement we need to transfer the solution vectors + // from the old mesh to the new one. To this end we use the + // 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). 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 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(previous_solution); + + // Now everything is ready, so do the refinement and recreate the DoF + // structure on the new grid, and finally initialize the matrix structures + // and the new vectors in the setup_system function. Next, we + // actually perform the interpolation of the solution from old to new + // grid. The final step is to apply the hanging node constraints to the + // solution vector, i.e., to make sure that the values of degrees of + // freedom located on hanging nodes are so that the solution is + // continuous. This is necessary since SolutionTransfer only operates on + // cells locally, without regard to the neighborhood. + triangulation.execute_coarsening_and_refinement(); + setup_system(); + + solution_trans.interpolate(previous_solution, solution); + constraints.distribute(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. + // + // @note If you're an experienced programmer, you may be surprised + // that we use a goto statement in this piece of code! + // goto statements are not particularly well liked any + // more since Edsgar Dijkstra, one of the greats of computer science, + // wrote a letter in 1968 called "Go To Statement considered harmful" + // (see here). + // The author of this code subscribes to this notion whole-heartedly: + // goto is hard to understand. In fact, deal.II contains + // virtually no occurrences: excluding code that was essentially + // transcribed from books and not counting duplicated code pieces, + // there are 3 locations in about 600,000 lines of code at the time + // this note is written; we also use it in 4 tutorial programs, in + // exactly the same context as here. Instead of trying to justify + // the occurrence here, let's first look at the code and we'll come + // back to the issue at the end of function. + template + void HeatEquation::run() + { + const unsigned int initial_global_refinement = 2; + const unsigned int n_adaptive_pre_refinement_steps = 4; + + GridGenerator::hyper_L(triangulation); + triangulation.refine_global(initial_global_refinement); + + setup_system(); + + unsigned int pre_refinement_step = 0; + + Vector tmp; + Vector forcing_terms; + + start_time_iteration: + + time = 0.0; + timestep_number = 0; + + tmp.reinit(solution.size()); + forcing_terms.reinit(solution.size()); + + + VectorTools::interpolate(dof_handler, + Functions::ZeroFunction(), + old_solution); + solution = old_solution; + + 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; + ++timestep_number; + + std::cout << "Time step " << timestep_number << " at t=" << time + << std::endl; + + 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, + QGauss(fe.degree + 1), + rhs_function, + tmp); + forcing_terms = tmp; + forcing_terms *= time_step * theta; + + rhs_function.set_time(time - time_step); + VectorTools::create_right_hand_side(dof_handler, + QGauss(fe.degree + 1), + rhs_function, + tmp); + + 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); + + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + boundary_values_function, + boundary_values); + + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + 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)) + { + refine_mesh(initial_global_refinement, + 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); + tmp.reinit(solution.size()); + forcing_terms.reinit(solution.size()); + } + + old_solution = solution; + } + } +} // namespace Step86 +// Now that you have seen what the function does, let us come back to the issue +// of the goto. In essence, what the code does is +// something like this: +// @code +// void run () +// { +// initialize; +// start_time_iteration: +// for (timestep=1...) +// { +// solve timestep; +// if (timestep==1 && not happy with the result) +// { +// adjust some data structures; +// goto start_time_iteration; // simply try again +// } +// postprocess; +// } +// } +// @endcode +// Here, the condition "happy with the result" is whether we'd like to keep +// the current mesh or would rather refine the mesh and start over on the +// new mesh. We could of course replace the use of the goto +// by the following: +// @code +// void run () +// { +// initialize; +// while (true) +// { +// solve timestep; +// if (not happy with the result) +// adjust some data structures; +// else +// break; +// } +// postprocess; +// +// for (timestep=2...) +// { +// solve timestep; +// postprocess; +// } +// } +// @endcode +// This has the advantage of getting rid of the goto +// but the disadvantage of having to duplicate the code that implements +// the "solve timestep" and "postprocess" operations in two different +// places. This could be countered by putting these parts of the code +// (sizable chunks in the actual implementation above) into their +// own functions, but a while(true) loop with a +// break statement is not really all that much easier +// to read or understand than a goto. +// +// In the end, one might simply agree that in general +// goto statements are a bad idea but be pragmatic and +// state that there may be occasions where they can help avoid code +// duplication and awkward control flow. This may be one of these +// places, and it matches the position Steve McConnell takes in his +// excellent book "Code Complete" @cite CodeComplete about good +// programming practices (see the mention of this book in the +// introduction of step-1) that spends a surprising ten pages on the +// question of goto in general. + + +// @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 + { + using namespace Step86; + + HeatEquation<2> heat_equation_solver; + heat_equation_solver.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; +} -- 2.39.5