From 921458fb7587f545013045021d6282665f26b806 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 26 Jul 2024 11:33:09 -0400 Subject: [PATCH] Remove step-52 --- doc/doxygen/tutorial/tutorial.h.in | 12 - doc/news/8.1.0-vs-8.2.0.h | 2 +- doc/news/9.1.1-vs-9.2.0.h | 2 +- examples/step-26/doc/intro.dox | 2 +- examples/step-26/doc/results.dox | 2 +- examples/step-52/CMakeLists.txt | 55 -- examples/step-52/doc/builds-on | 1 - examples/step-52/doc/intro.dox | 260 ------ examples/step-52/doc/kind | 1 - examples/step-52/doc/results.dox | 39 - examples/step-52/doc/tooltip | 1 - examples/step-52/step-52.cc | 836 ------------------ tests/examples/step-52.diff | 74 -- .../examples/step-52.with_umfpack=true.output | 16 - 14 files changed, 4 insertions(+), 1299 deletions(-) delete mode 100644 examples/step-52/CMakeLists.txt delete mode 100644 examples/step-52/doc/builds-on delete mode 100644 examples/step-52/doc/intro.dox delete mode 100644 examples/step-52/doc/kind delete mode 100644 examples/step-52/doc/results.dox delete mode 100644 examples/step-52/doc/tooltip delete mode 100644 examples/step-52/step-52.cc delete mode 100644 tests/examples/step-52.diff delete mode 100644 tests/examples/step-52.with_umfpack=true.output diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 9b8f63bfc7..a8a0051454 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -477,12 +477,6 @@ * * * - * step-52 - * Solving the time dependent neutron diffusion equation using - * Runge-Kutta methods. Explicit and implicit time stepping. - * - * - * * step-53 * Describing the geometry of complex domains and curved boundaries. * @@ -1506,12 +1500,6 @@ * * * - * Time dependent neutron diffusion equation - * - * step-52 - * - * - * * Time dependent fluid structure interaction problems * * step-70 diff --git a/doc/news/8.1.0-vs-8.2.0.h b/doc/news/8.1.0-vs-8.2.0.h index 708f8e840d..5a6dddc369 100644 --- a/doc/news/8.1.0-vs-8.2.0.h +++ b/doc/news/8.1.0-vs-8.2.0.h @@ -163,7 +163,7 @@ inconvenience this causes. (Luca Heltai, 2014/09/13) -
  • New: The new tutorial program step-52 explains how to use the +
  • New: The new tutorial program step 52 explains how to use the new time stepping methods.
    (Bruno Turcksin, Damien Lebrun-Grandie, 2014/09/12) diff --git a/doc/news/9.1.1-vs-9.2.0.h b/doc/news/9.1.1-vs-9.2.0.h index cbcc535151..aa1c1e9651 100644 --- a/doc/news/9.1.1-vs-9.2.0.h +++ b/doc/news/9.1.1-vs-9.2.0.h @@ -989,7 +989,7 @@ inconvenience this causes.
  • - Fixed: step-26 and step-52 now also output the simulation time + Fixed: step-26 and step 52 now also output the simulation time in the VTK/VTU output files.
    (Wolfgang Bangerth, 2020/03/27) diff --git a/examples/step-26/doc/intro.dox b/examples/step-26/doc/intro.dox index 44115e6998..ae3564e246 100644 --- a/examples/step-26/doc/intro.dox +++ b/examples/step-26/doc/intro.dox @@ -52,7 +52,7 @@ the explicit Euler ($\theta=0$), implicit Euler ($\theta=1$) and Crank-Nicolson ($\theta=\frac 12$) time discretizations. Since the latter has the highest convergence order, we will choose $\theta=\frac 12$ in the program below, but make it so that playing with this parameter remains simple. (If you -are interested in playing with higher order methods, take a look at step-52.) +are interested in playing with higher order methods, take a look at step-86.) Given this time discretization, space discretization happens as it always does, by multiplying with test functions, integrating by parts, and then diff --git a/examples/step-26/doc/results.dox b/examples/step-26/doc/results.dox index 36cf16c28c..d8b4f2172b 100644 --- a/examples/step-26/doc/results.dox +++ b/examples/step-26/doc/results.dox @@ -90,7 +90,7 @@ more widely used (and even less accurate) implicit (=backward) Euler method, but many other, more accurate methods such as BDF or Runge-Kutta methods are available and should be used as they do not represent much additional effort. It is not difficult to implement this for the current -program, if one wanted; a more systematic treatment is also given in step-52. +program, if one wanted; a more systematic treatment is also given in step-86. As a general rule, however, one should not be implementing time stepping methods by hand, as we do here, for problems that do not require diff --git a/examples/step-52/CMakeLists.txt b/examples/step-52/CMakeLists.txt deleted file mode 100644 index 4d39b4173d..0000000000 --- a/examples/step-52/CMakeLists.txt +++ /dev/null @@ -1,55 +0,0 @@ -## -# CMake script for the step-52 tutorial program: -## - -# Set the name of the project and target: -set(TARGET "step-52") - -# 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 - ) - -# Define the output that should be cleaned: -set(CLEAN_UP_FILES *.vtu) - -# Usually, you will not need to modify anything beyond this point... - -cmake_minimum_required(VERSION 3.13.4) - -find_package(deal.II 9.7.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() - -# -# Are all dependencies fulfilled? -# -if(NOT DEAL_II_WITH_UMFPACK) - message(FATAL_ERROR " -Error! This tutorial requires a deal.II library that was configured with the following options: - DEAL_II_WITH_UMFPACK = ON -However, the deal.II library found at ${DEAL_II_PATH} was configured with these options: - DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK} -This conflicts with the requirements." - ) -endif() - -deal_ii_initialize_cached_variables() -project(${TARGET}) -deal_ii_invoke_autopilot() diff --git a/examples/step-52/doc/builds-on b/examples/step-52/doc/builds-on deleted file mode 100644 index 1aabbdfdd6..0000000000 --- a/examples/step-52/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-26 diff --git a/examples/step-52/doc/intro.dox b/examples/step-52/doc/intro.dox deleted file mode 100644 index dc7b1d4587..0000000000 --- a/examples/step-52/doc/intro.dox +++ /dev/null @@ -1,260 +0,0 @@ -
    - -This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie. - -@note In order to run this program, deal.II must be configured to use -the UMFPACK sparse direct solver. Refer to the ReadMe for instructions how to do this. - - -

    Introduction

    - -This program shows how to use Runge-Kutta methods to solve a time-dependent -problem. It solves a small variation of the heat equation discussed first in -step-26 but, since the purpose of this program is only to demonstrate using -more advanced ways to interface with deal.II's time stepping algorithms, only -solves a simple problem on a uniformly refined mesh. - -@note At the end of the day, time stepping is a problem that is only - solved efficiently if you use adaptive time step selection and error - control. This is implemented in many external libraries, but not in - deal.II itself. You may want to consider looking at step-86 for a - worked out example. In the meantime, this program is deprecated and - will be removed after the deal.II 9.6 release. - - - -

    Problem statement

    - -In this example, we solve the one-group time-dependent diffusion -approximation of the neutron transport equation (see step-28 for the -time-independent multigroup diffusion). This is a model for how neutrons move -around highly scattering media, and consequently it is a variant of the -time-dependent diffusion equation -- which is just a different name for the -heat equation discussed in step-26, plus some extra terms. -We assume that the medium is not -fissible and therefore, the neutron flux satisfies the following equation: -@f{eqnarray*}{ -\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla \cdot D(x) \nabla \phi(x,t) -- \Sigma_a(x) \phi(x,t) + S(x,t) -@f} -augmented by appropriate boundary conditions. Here, $v$ is the velocity of -neutrons (for simplicity we assume it is equal to 1 which can be achieved by -simply scaling the time variable), $D$ is the diffusion coefficient, -$\Sigma_a$ is the absorption cross section, and $S$ is a source. Because we are -only interested in the time dependence, we assume that $D$ and $\Sigma_a$ are -constant. - -Since this program only intends to demonstrate how to use advanced time -stepping algorithms, we will only look for the solutions of relatively simple -problems. Specifically, we are looking for a solution on a square domain -$[0,b]\times[0,b]$ of the form -@f{eqnarray*}{ -\phi(x,t) = A\sin(\omega t)(bx-x^2). -@f} -By using quadratic finite elements, we can represent this function exactly at -any particular time, and all the error will be due to the time -discretization. We do this because it is then easy to observe the order of -convergence of the various time stepping schemes we will consider, without -having to separate spatial and temporal errors. - -We impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and -$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. We choose the -source term so that the corresponding solution is -in fact of the form stated above: -@f{eqnarray*}{ -S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t) -\left(\Sigma_a (bx-x^2)+2D\right) \right). -@f} -Because the solution is a sine in time, we know that the exact solution -satisfies $\phi\left(x,\frac{\pi}{\omega}\right) = 0$. -Therefore, the error at time $t=\frac{\pi}{\omega}$ is simply the norm of the numerical -solution, i.e., $\|e(\cdot,t=\frac{\pi}{\omega})\|_{L_2} = \|\phi_h(\cdot,t=\frac{\pi}{\omega})\|_{L_2}$, -and is particularly easily evaluated. In the code, we evaluate the $l_2$ norm -of the vector of nodal values of $\phi_h$ instead of the $L_2$ norm of the -associated spatial function, since the former is simpler to compute; however, -on uniform meshes, the two are just related by a constant and we can -consequently observe the temporal convergence order with either. - - -

    Runge-Kutta methods

    - -The Runge-Kutta methods implemented in deal.II assume that the equation to be -solved can be written as: -@f{eqnarray*}{ -\frac{dy}{dt} = g(t,y). -@f} -On the other hand, when using finite elements, discretized time derivatives always result in the -presence of a @ref GlossMassMatrix "mass matrix" on the left hand side. This can easily be seen by -considering that if the solution vector $y(t)$ in the equation above is in fact the vector -of nodal coefficients $U(t)$ for a variable of the form -@f{eqnarray*}{ - u_h(x,t) = \sum_j U_j(t) \varphi_j(x) -@f} -with spatial shape functions $\varphi_j(x)$, then multiplying an equation of -the form -@f{eqnarray*}{ - \frac{\partial u(x,t)}{\partial t} = q(t,u(x,t)) -@f} -by test functions, integrating over $\Omega$, substituting $u\rightarrow u_h$ -and restricting the test functions to the $\varphi_i(x)$ from above, then this -spatially discretized equation has the form -@f{eqnarray*}{ -M\frac{dU}{dt} = f(t,U), -@f} -where $M$ is the mass matrix and $f(t,U)$ is the spatially discretized version -of $q(t,u(x,t))$ (where $q$ is typically the place where spatial -derivatives appear, but this is not of much concern for the moment given that -we only consider time derivatives). In other words, this form fits the general -scheme above if we write -@f{eqnarray*}{ -\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y). -@f} - -Runke-Kutta methods are time stepping schemes that approximate $y(t_n)\approx -y_{n}$ through a particular one-step approach. They are typically written in the form -@f{eqnarray*}{ -y_{n+1} = y_n + \sum_{i=1}^s b_i k_i -@f} -where for the form of the right hand side above -@f{eqnarray*}{ -k_i = \Delta t \, M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right). -@f} -Here $a_{ij}$, $b_i$, and $c_i$ are known coefficients that identify which -particular Runge-Kutta scheme you want to use, and $\Delta t=t_{n+1}-t_n$ is the time step -used. Different time stepping methods of the Runge-Kutta class differ in the -number of stages $s$ and the values they use for the coefficients $a_{ij}$, -$b_i$, and $c_i$ but are otherwise easy to implement since one can look up -tabulated values for these coefficients. (These tables are often called -Butcher tableaus.) - -At the time of the writing of this tutorial, the methods implemented in -deal.II can be divided in three categories: -
      -
    1. Explicit Runge-Kutta; in order for a method to be explicit, it is -necessary that in the formula above defining $k_i$, $k_i$ does not appear -on the right hand side. In other words, these methods have to satisfy -$a_{ii}=0, i=1,\ldots,s$. -
    2. Embedded (or adaptive) Runge-Kutta; we will discuss their properties below. -
    3. Implicit Runge-Kutta; this class of methods require the solution of a -possibly nonlinear system the stages $k_i$ above, i.e., they have -$a_{ii}\neq 0$ for at least one of the stages $i=1,\ldots,s$. -
    -Many well known time stepping schemes that one does not typically associate -with the names Runge or Kutta can in fact be written in a way so that they, -too, can be expressed in these categories. They oftentimes represent the -lowest-order members of these families; one example is the simple explicit -Euler method. - - -

    Explicit Runge-Kutta methods

    - -These methods, only require a function to evaluate $M^{-1}f(t,y)$ but not -(as implicit methods) to solve an equation that involves -$f(t,y)$ for $y$. As all explicit time stepping methods, they become unstable -when the time step chosen is too large. - -Well known methods in this class include forward Euler, third order -Runge-Kutta, and fourth order Runge-Kutta (often abbreviated as RK4). - - -

    Embedded Runge-Kutta methods

    - -These methods use both a lower and a higher order method to -estimate the error and decide if the time step needs to be shortened or can be -increased. The term "embedded" refers to the fact that the lower-order method -does not require additional evaluates of the function $M^{-1}f(\cdot,\cdot)$ -but reuses data that has to be computed for the high order method anyway. It -is, in other words, essentially free, and we get the error estimate as a side -product of using the higher order method. - -This class of methods include Heun-Euler, Bogacki-Shampine, Dormand-Prince (ode45 in -Matlab and often abbreviated as RK45 to indicate that the lower and higher order methods -used here are 4th and 5th order Runge-Kutta methods, respectively), Fehlberg, -and Cash-Karp. - -At the time of the writing, only embedded explicit methods have been implemented. - - -

    Implicit Runge-Kutta methods

    - -Implicit methods require the solution of (possibly nonlinear) systems of the -form $\alpha y = f(t,y)$ -for $y$ in each (sub-)timestep. Internally, this is -done using a Newton-type method and, consequently, they require that the user -provide functions that can evaluate $M^{-1}f(t,y)$ and -$\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently -$\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$. - -The particular form of this operator results from the fact that each Newton -step requires the solution of an equation of the form -@f{align*}{ - \left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y - = -M h(t,y) -@f} -for some (given) $h(t,y)$. Implicit methods are -always stable, regardless of the time step size, but too large time steps of -course affect the accuracy of the solution, even if the numerical -solution remains stable and bounded. - -Methods in this class include backward Euler, implicit midpoint, -Crank-Nicolson, and the two stage SDIRK method (short for "singly diagonally -implicit Runge-Kutta", a term coined to indicate that the diagonal elements -$a_{ii}$ defining the time stepping method are all equal; this property -allows for the Newton matrix $I-\tau M^{-1}\frac{\partial f}{\partial y}$ to -be re-used between stages because $\tau$ is the same every time). - - -

    Spatially discrete formulation

    - -By expanding the solution of our model problem -as always using shape functions $\psi_j$ and writing -@f{eqnarray*}{ -\phi_h(x,t) = \sum_j U_j(t) \psi_j(x), -@f} -we immediately get the spatially discretized version of the diffusion equation as -@f{eqnarray*}{ - M \frac{dU(t)}{dt} - = -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t) -@f} -where -@f{eqnarray*}{ - M_{ij} &=& (\psi_i,\psi_j), \\ - {\cal D}_{ij} &=& (D\nabla\psi_i,\nabla\psi_j)_\Omega, \\ - {\cal A}_{ij} &=& (\Sigma_a\psi_i,\psi_j)_\Omega, \\ - {\cal S}_{i}(t) &=& (\psi_i,S(x,t))_\Omega. -@f} -See also step-24 and step-26 to understand how we arrive here. -Boundary terms are not necessary due to the chosen boundary conditions for -the current problem. To use the Runge-Kutta methods, we recast this -as follows: -@f{eqnarray*}{ -f(y) = -{\cal D}y - {\cal A}y + {\cal S}. -@f} -In the code, we will need to be able to evaluate this function $f(U)$ along -with its derivative, -@f{eqnarray*}{ -\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}. -@f} - - -

    Notes on the testcase

    - -To simplify the problem, the domain is two dimensional and the mesh is -uniformly refined (there is no need to adapt the mesh since we use quadratic -finite elements and the exact solution is quadratic). Going from a two -dimensional domain to a three dimensional domain is not very -challenging. However if you intend to solve more complex problems where the -mesh must be adapted (as is done, for example, in step-26), then it is -important to remember the following issues: - -
      -
    1. You will need to project the solution to the new mesh when the mesh is changed. Of course, - the mesh - used should be the same from the beginning to the end of each time step, - a question that arises because Runge-Kutta methods use multiple - evaluations of the equations within each time step. -
    2. You will need to update the mass matrix and its inverse every time the - mesh is changed. -
    -The techniques for these steps are readily available by looking at step-26. diff --git a/examples/step-52/doc/kind b/examples/step-52/doc/kind deleted file mode 100644 index 86a44aa1ef..0000000000 --- a/examples/step-52/doc/kind +++ /dev/null @@ -1 +0,0 @@ -time dependent diff --git a/examples/step-52/doc/results.dox b/examples/step-52/doc/results.dox deleted file mode 100644 index 03ea68c15e..0000000000 --- a/examples/step-52/doc/results.dox +++ /dev/null @@ -1,39 +0,0 @@ -

    Results

    - -The point of this program is less to show particular results, but instead to -show how it is done. This we have already demonstrated simply by discussing -the code above. Consequently, the output the program yields is relatively -sparse and consists only of the console output and the solutions given in VTU -format for visualization. - -The console output contains both errors and, for some of the methods, the -number of steps they performed: -@code -Explicit methods: - Forward Euler: error=1.00883 - Third order Runge-Kutta: error=0.000227982 - Fourth order Runge-Kutta: error=1.90541e-06 - -Implicit methods: - Backward Euler: error=1.03428 - Implicit Midpoint: error=0.00862702 - Crank-Nicolson: error=0.00862675 - SDIRK: error=0.0042349 - -Embedded explicit methods: - Heun-Euler: error=0.0073012 - steps performed=284 - Bogacki-Shampine: error=0.000408407 - steps performed=181 - Dopri: error=0.000836695 - steps performed=120 - Fehlberg: error=0.00248922 - steps performed=106 - Cash-Karp: error=0.0787735 - steps performed=106 -@endcode - -As expected the higher order methods give (much) more accurate solutions. We -also see that the (rather inaccurate) Heun-Euler method increased the number of -time steps in order to satisfy the tolerance. On the other hand, the other -embedded methods used a lot less time steps than what was prescribed. diff --git a/examples/step-52/doc/tooltip b/examples/step-52/doc/tooltip deleted file mode 100644 index 7754b44e6b..0000000000 --- a/examples/step-52/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -Time-dependent diffusion equation. Time stepping methods. diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc deleted file mode 100644 index 6c610891e7..0000000000 --- a/examples/step-52/step-52.cc +++ /dev/null @@ -1,836 +0,0 @@ -/* ------------------------------------------------------------------------ - * - * SPDX-License-Identifier: LGPL-2.1-or-later - * Copyright (C) 2014 - 2024 by the deal.II authors - * - * This file is part of the deal.II library. - * - * Part of the source code is dual licensed under Apache-2.0 WITH - * LLVM-exception OR LGPL-2.1-or-later. Detailed license information - * governing the source code and code contributions can be found in - * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. - * - * ------------------------------------------------------------------------ - * - * Authors: Damien Lebrun-Grandie, Bruno Turcksin, 2014 - */ - -// @sect3{Include files} - -// The first task as usual is to include the functionality of these well-known -// deal.II library files and some C++ header files. -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include - -#include -#include - -#include -#include - -#include -#include -#include -#include - -// This is the only include file that is new: It includes all the Runge-Kutta -// methods. -#include - - -// The next step is like in all previous tutorial programs: We put everything -// into a namespace of its own and then import the deal.II classes and functions -// into it. -namespace Step52 -{ - using namespace dealii; - - // @sect3{The Diffusion class} - - // The next piece is the declaration of the main class. Most of the - // functions in this class are not new and have been explained in previous - // tutorials. The only interesting functions are - // evaluate_diffusion() and - // id_minus_tau_J_inverse(). evaluate_diffusion() - // evaluates the diffusion equation, $M^{-1}(f(t,y))$, at a given time and a - // given $y$. id_minus_tau_J_inverse() evaluates $\left(I-\tau - // M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or equivalently - // $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at a given - // time, for a given $\tau$ and $y$. This function is needed when an - // implicit method is used. - class Diffusion - { - public: - Diffusion(); - - void run(); - - private: - void setup_system(); - - void assemble_system(); - - double get_source(const double time, const Point<2> &point) const; - - Vector evaluate_diffusion(const double time, - const Vector &y) const; - - Vector id_minus_tau_J_inverse(const double time, - const double tau, - const Vector &y); - - void output_results(const double time, - const unsigned int time_step, - TimeStepping::runge_kutta_method method) const; - - // The next three functions are the drivers for the explicit methods, the - // implicit methods, and the embedded explicit methods respectively. The - // driver function for embedded explicit methods returns the number of - // steps executed given that it only takes the number of time steps passed - // as an argument as a hint, but internally computed the optimal time step - // itself. - void explicit_method(const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - void implicit_method(const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - unsigned int - embedded_explicit_method(const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - - const unsigned int fe_degree; - - const double diffusion_coefficient; - const double absorption_cross_section; - - Triangulation<2> triangulation; - - const FE_Q<2> fe; - - DoFHandler<2> dof_handler; - - AffineConstraints constraint_matrix; - - SparsityPattern sparsity_pattern; - - SparseMatrix system_matrix; - SparseMatrix mass_matrix; - SparseMatrix mass_minus_tau_Jacobian; - - SparseDirectUMFPACK inverse_mass_matrix; - - Vector solution; - }; - - - - // We choose quadratic finite elements and we initialize the parameters. - Diffusion::Diffusion() - : fe_degree(2) - , diffusion_coefficient(1. / 30.) - , absorption_cross_section(1.) - , fe(fe_degree) - , dof_handler(triangulation) - {} - - - - // @sect4{Diffusion::setup_system} - // Now, we create the constraint matrix and the sparsity pattern. Then, we - // initialize the matrices and the solution vector. - void Diffusion::setup_system() - { - dof_handler.distribute_dofs(fe); - - VectorTools::interpolate_boundary_values(dof_handler, - 1, - Functions::ZeroFunction<2>(), - constraint_matrix); - constraint_matrix.close(); - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp, constraint_matrix); - sparsity_pattern.copy_from(dsp); - - system_matrix.reinit(sparsity_pattern); - mass_matrix.reinit(sparsity_pattern); - mass_minus_tau_Jacobian.reinit(sparsity_pattern); - solution.reinit(dof_handler.n_dofs()); - } - - - - // @sect4{Diffusion::assemble_system} - // In this function, we compute $-\int D \nabla b_i \cdot \nabla b_j - // d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}$ and the mass - // matrix $\int b_i b_j d\boldsymbol{r}$. The @ref GlossMassMatrix "mass matrix" is then - // inverted using a direct solver; the inverse_mass_matrix - // variable will then store the inverse of the mass matrix so that - // $M^{-1}$ can be applied to a vector using the vmult() - // function of that object. (Internally, UMFPACK does not really store - // the inverse of the matrix, but its LU factors; applying the inverse - // matrix is then equivalent to doing one forward and one backward solves - // with these two factors, which has the same complexity as applying an - // explicit inverse of the matrix). - void Diffusion::assemble_system() - { - system_matrix = 0.; - mass_matrix = 0.; - - const QGauss<2> quadrature_formula(fe_degree + 1); - - FEValues<2> fe_values(fe, - quadrature_formula, - update_values | update_gradients | update_JxW_values); - - - const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - FullMatrix cell_mass_matrix(dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices(dofs_per_cell); - - for (const auto &cell : dof_handler.active_cell_iterators()) - { - cell_matrix = 0.; - cell_mass_matrix = 0.; - - fe_values.reinit(cell); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - { - cell_matrix(i, j) += - ((-diffusion_coefficient * // (-D - fe_values.shape_grad(i, q_point) * // * grad phi_i - fe_values.shape_grad(j, q_point) // * grad phi_j - - absorption_cross_section * // -Sigma - fe_values.shape_value(i, q_point) * // * phi_i - fe_values.shape_value(j, q_point)) // * phi_j) - * fe_values.JxW(q_point)); // * dx - cell_mass_matrix(i, j) += fe_values.shape_value(i, q_point) * - fe_values.shape_value(j, q_point) * - fe_values.JxW(q_point); - } - - cell->get_dof_indices(local_dof_indices); - - constraint_matrix.distribute_local_to_global(cell_matrix, - local_dof_indices, - system_matrix); - constraint_matrix.distribute_local_to_global(cell_mass_matrix, - local_dof_indices, - mass_matrix); - } - - inverse_mass_matrix.initialize(mass_matrix); - } - - - - // @sect4{Diffusion::get_source} - // - // In this function, the source term of the equation for a given time and a - // given point is computed. - double Diffusion::get_source(const double time, const Point<2> &point) const - { - const double intensity = 10.; - const double frequency = numbers::PI / 10.; - const double b = 5.; - const double x = point[0]; - - return intensity * - (frequency * std::cos(frequency * time) * (b * x - x * x) + - std::sin(frequency * time) * - (absorption_cross_section * (b * x - x * x) + - 2. * diffusion_coefficient)); - } - - - - // @sect4{Diffusion::evaluate_diffusion} - // - // Next, we evaluate the weak form of the diffusion equation at a given time - // $t$ and for a given vector $y$. In other words, as outlined in the - // introduction, we evaluate $M^{-1}(-{\cal D}y - {\cal A}y + {\cal - // S})$. For this, we have to apply the matrix $-{\cal D} - {\cal A}$ - // (previously computed and stored in the variable - // system_matrix) to $y$ and then add the source term which we - // integrate as we usually do. (Integrating up the solution could be done - // using VectorTools::create_right_hand_side() if you wanted to save a few - // lines of code, or wanted to take advantage of doing the integration in - // parallel.) The result is then multiplied by $M^{-1}$. - Vector Diffusion::evaluate_diffusion(const double time, - const Vector &y) const - { - Vector tmp(dof_handler.n_dofs()); - tmp = 0.; - system_matrix.vmult(tmp, y); - - const QGauss<2> quadrature_formula(fe_degree + 1); - - FEValues<2> fe_values(fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values); - - - const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); - const unsigned int n_q_points = quadrature_formula.size(); - - Vector cell_source(dofs_per_cell); - - std::vector local_dof_indices(dofs_per_cell); - - for (const auto &cell : dof_handler.active_cell_iterators()) - { - cell_source = 0.; - - fe_values.reinit(cell); - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - const double source = - get_source(time, fe_values.quadrature_point(q_point)); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - cell_source(i) += fe_values.shape_value(i, q_point) * // phi_i(x) - source * // * S(x) - fe_values.JxW(q_point); // * dx - } - - cell->get_dof_indices(local_dof_indices); - - constraint_matrix.distribute_local_to_global(cell_source, - local_dof_indices, - tmp); - } - - Vector value(dof_handler.n_dofs()); - inverse_mass_matrix.vmult(value, tmp); - - return value; - } - - - // @sect4{Diffusion::id_minus_tau_J_inverse} - // - // We compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$. This - // is done in several steps: - // - compute $M-\tau \frac{\partial f}{\partial y}$ - // - invert the matrix to get $\left(M-\tau \frac{\partial f} - // {\partial y}\right)^{-1}$ - // - compute $tmp=My$ - // - compute $z=\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = - // \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$ - // - return z. - Vector Diffusion::id_minus_tau_J_inverse(const double /*time*/, - const double tau, - const Vector &y) - { - SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian; - - mass_minus_tau_Jacobian.copy_from(mass_matrix); - mass_minus_tau_Jacobian.add(-tau, system_matrix); - - inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian); - - Vector tmp(dof_handler.n_dofs()); - mass_matrix.vmult(tmp, y); - - Vector result(y); - inverse_mass_minus_tau_Jacobian.vmult(result, tmp); - - return result; - } - - - - // @sect4{Diffusion::output_results} - // - // The following function then outputs the solution in vtu files indexed by - // the number of the time step and the name of the time stepping method. Of - // course, the (exact) result should really be the same for all time - // stepping method, but the output here at least allows us to compare them. - void Diffusion::output_results(const double time, - const unsigned int time_step, - TimeStepping::runge_kutta_method method) const - { - std::string method_name; - - switch (method) - { - case TimeStepping::FORWARD_EULER: - { - method_name = "forward_euler"; - break; - } - case TimeStepping::RK_THIRD_ORDER: - { - method_name = "rk3"; - break; - } - case TimeStepping::RK_CLASSIC_FOURTH_ORDER: - { - method_name = "rk4"; - break; - } - case TimeStepping::BACKWARD_EULER: - { - method_name = "backward_euler"; - break; - } - case TimeStepping::IMPLICIT_MIDPOINT: - { - method_name = "implicit_midpoint"; - break; - } - case TimeStepping::SDIRK_TWO_STAGES: - { - method_name = "sdirk"; - break; - } - case TimeStepping::HEUN_EULER: - { - method_name = "heun_euler"; - break; - } - case TimeStepping::BOGACKI_SHAMPINE: - { - method_name = "bogacki_shampine"; - break; - } - case TimeStepping::DOPRI: - { - method_name = "dopri"; - break; - } - case TimeStepping::FEHLBERG: - { - method_name = "fehlberg"; - break; - } - case TimeStepping::CASH_KARP: - { - method_name = "cash_karp"; - break; - } - default: - { - break; - } - } - - DataOut<2> data_out; - - data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "solution"); - - data_out.build_patches(); - - data_out.set_flags(DataOutBase::VtkFlags(time, time_step)); - - const std::string filename = "solution_" + method_name + "-" + - Utilities::int_to_string(time_step, 3) + - ".vtu"; - std::ofstream output(filename); - data_out.write_vtu(output); - - static std::vector> times_and_names; - - static std::string method_name_prev = ""; - static std::string pvd_filename; - if (method_name_prev != method_name) - { - times_and_names.clear(); - method_name_prev = method_name; - pvd_filename = "solution_" + method_name + ".pvd"; - } - times_and_names.emplace_back(time, filename); - std::ofstream pvd_output(pvd_filename); - DataOutBase::write_pvd_record(pvd_output, times_and_names); - } - - - // @sect4{Diffusion::explicit_method} - // - // This function is the driver for all the explicit methods. At the - // top it initializes the time stepping and the solution (by setting - // it to zero and then ensuring that boundary value and hanging node - // constraints are respected; of course, with the mesh we use here, - // hanging node constraints are not in fact an issue). It then calls - // evolve_one_time_step which performs one time step. - // Time is stored and incremented through a DiscreteTime object. - // - // For explicit methods, evolve_one_time_step needs to - // evaluate $M^{-1}(f(t,y))$, i.e, it needs - // evaluate_diffusion. Because - // evaluate_diffusion is a member function, it needs to - // be bound to this. After each evolution step, we - // again apply the correct boundary values and hanging node - // constraints. - // - // Finally, the solution is output - // every 10 time steps. - void Diffusion::explicit_method(const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) - { - const double time_step = - (final_time - initial_time) / static_cast(n_time_steps); - - solution = 0.; - constraint_matrix.distribute(solution); - - TimeStepping::ExplicitRungeKutta> explicit_runge_kutta( - method); - output_results(initial_time, 0, method); - DiscreteTime time(initial_time, final_time, time_step); - while (time.is_at_end() == false) - { - explicit_runge_kutta.evolve_one_time_step( - [this](const double time, const Vector &y) { - return this->evaluate_diffusion(time, y); - }, - time.get_current_time(), - time.get_next_step_size(), - solution); - time.advance_time(); - - constraint_matrix.distribute(solution); - - if (time.get_step_number() % 10 == 0) - output_results(time.get_current_time(), - time.get_step_number(), - method); - } - } - - - - // @sect4{Diffusion::implicit_method} - // This function is equivalent to explicit_method but for - // implicit methods. When using implicit methods, we need to evaluate - // $M^{-1}(f(t,y))$ and $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial - // y}\right)^{-1}$ for which we use the two member functions previously - // introduced. - void Diffusion::implicit_method(const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) - { - const double time_step = - (final_time - initial_time) / static_cast(n_time_steps); - - solution = 0.; - constraint_matrix.distribute(solution); - - TimeStepping::ImplicitRungeKutta> implicit_runge_kutta( - method); - output_results(initial_time, 0, method); - DiscreteTime time(initial_time, final_time, time_step); - while (time.is_at_end() == false) - { - implicit_runge_kutta.evolve_one_time_step( - [this](const double time, const Vector &y) { - return this->evaluate_diffusion(time, y); - }, - [this](const double time, const double tau, const Vector &y) { - return this->id_minus_tau_J_inverse(time, tau, y); - }, - time.get_current_time(), - time.get_next_step_size(), - solution); - time.advance_time(); - - constraint_matrix.distribute(solution); - - if (time.get_step_number() % 10 == 0) - output_results(time.get_current_time(), - time.get_step_number(), - method); - } - } - - - - // @sect4{Diffusion::embedded_explicit_method} - // This function is the driver for the embedded explicit methods. It requires - // more parameters: - // - coarsen_param: factor multiplying the current time step when the error - // is below the threshold. - // - refine_param: factor multiplying the current time step when the error - // is above the threshold. - // - min_delta: smallest time step acceptable. - // - max_delta: largest time step acceptable. - // - refine_tol: threshold above which the time step is refined. - // - coarsen_tol: threshold below which the time step is coarsen. - // - // Embedded methods use a guessed time step. If the error using this time step - // is too large, the time step will be reduced. If the error is below the - // threshold, a larger time step will be tried for the next time step. - // delta_t_guess is the guessed time step produced by the - // embedded method. In summary, time step size is potentially modified in - // three ways: - // - Reducing or increasing time step size within - // TimeStepping::EmbeddedExplicitRungeKutta::evolve_one_time_step(). - // - Using the calculated delta_t_guess. - // - Automatically adjusting the step size of the last time step to ensure - // simulation ends precisely at final_time. This adjustment - // is handled inside the DiscreteTime instance. - unsigned int Diffusion::embedded_explicit_method( - const TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) - { - const double time_step = - (final_time - initial_time) / static_cast(n_time_steps); - const double coarsen_param = 1.2; - const double refine_param = 0.8; - const double min_delta = 1e-8; - const double max_delta = 10 * time_step; - const double refine_tol = 1e-1; - const double coarsen_tol = 1e-5; - - solution = 0.; - constraint_matrix.distribute(solution); - - TimeStepping::EmbeddedExplicitRungeKutta> - embedded_explicit_runge_kutta(method, - coarsen_param, - refine_param, - min_delta, - max_delta, - refine_tol, - coarsen_tol); - output_results(initial_time, 0, method); - DiscreteTime time(initial_time, final_time, time_step); - while (time.is_at_end() == false) - { - const double new_time = - embedded_explicit_runge_kutta.evolve_one_time_step( - [this](const double time, const Vector &y) { - return this->evaluate_diffusion(time, y); - }, - time.get_current_time(), - time.get_next_step_size(), - solution); - time.set_next_step_size(new_time - time.get_current_time()); - time.advance_time(); - - constraint_matrix.distribute(solution); - - if (time.get_step_number() % 10 == 0) - output_results(time.get_current_time(), - time.get_step_number(), - method); - - time.set_desired_next_step_size( - embedded_explicit_runge_kutta.get_status().delta_t_guess); - } - - return time.get_step_number(); - } - - - - // @sect4{Diffusion::run} - // - // The following is the main function of the program. At the top, we create - // the grid (a $[0,5]\times [0,5]$ square) and refine it four times to get a - // mesh that has 16 by 16 cells, for a total of 256. We then set the boundary - // indicator to 1 for those parts of the boundary where $x=0$ and $x=5$. - void Diffusion::run() - { - GridGenerator::hyper_cube(triangulation, 0., 5.); - triangulation.refine_global(4); - - for (const auto &cell : triangulation.active_cell_iterators()) - for (const auto &face : cell->face_iterators()) - if (face->at_boundary()) - { - if ((face->center()[0] == 0.) || (face->center()[0] == 5.)) - face->set_boundary_id(1); - else - face->set_boundary_id(0); - } - - // Next, we set up the linear systems and fill them with content so that - // they can be used throughout the time stepping process: - setup_system(); - - assemble_system(); - - // Finally, we solve the diffusion problem using several of the - // Runge-Kutta methods implemented in namespace TimeStepping, each time - // outputting the error at the end time. (As explained in the - // introduction, since the exact solution is zero at the final time, the - // error equals the numerical solution and can be computed by just taking - // the $l_2$ norm of the solution vector.) - unsigned int n_steps = 0; - const unsigned int n_time_steps = 200; - const double initial_time = 0.; - const double final_time = 10.; - - std::cout << "Explicit methods:" << std::endl; - explicit_method(TimeStepping::FORWARD_EULER, - n_time_steps, - initial_time, - final_time); - std::cout << " Forward Euler: error=" << solution.l2_norm() - << std::endl; - - explicit_method(TimeStepping::RK_THIRD_ORDER, - n_time_steps, - initial_time, - final_time); - std::cout << " Third order Runge-Kutta: error=" << solution.l2_norm() - << std::endl; - - explicit_method(TimeStepping::RK_CLASSIC_FOURTH_ORDER, - n_time_steps, - initial_time, - final_time); - std::cout << " Fourth order Runge-Kutta: error=" << solution.l2_norm() - << std::endl; - std::cout << std::endl; - - - std::cout << "Implicit methods:" << std::endl; - implicit_method(TimeStepping::BACKWARD_EULER, - n_time_steps, - initial_time, - final_time); - std::cout << " Backward Euler: error=" << solution.l2_norm() - << std::endl; - - implicit_method(TimeStepping::IMPLICIT_MIDPOINT, - n_time_steps, - initial_time, - final_time); - std::cout << " Implicit Midpoint: error=" << solution.l2_norm() - << std::endl; - - implicit_method(TimeStepping::CRANK_NICOLSON, - n_time_steps, - initial_time, - final_time); - std::cout << " Crank-Nicolson: error=" << solution.l2_norm() - << std::endl; - - implicit_method(TimeStepping::SDIRK_TWO_STAGES, - n_time_steps, - initial_time, - final_time); - std::cout << " SDIRK: error=" << solution.l2_norm() - << std::endl; - std::cout << std::endl; - - - std::cout << "Embedded explicit methods:" << std::endl; - n_steps = embedded_explicit_method(TimeStepping::HEUN_EULER, - n_time_steps, - initial_time, - final_time); - std::cout << " Heun-Euler: error=" << solution.l2_norm() - << std::endl; - std::cout << " steps performed=" << n_steps << std::endl; - - n_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE, - n_time_steps, - initial_time, - final_time); - std::cout << " Bogacki-Shampine: error=" << solution.l2_norm() - << std::endl; - std::cout << " steps performed=" << n_steps << std::endl; - - n_steps = embedded_explicit_method(TimeStepping::DOPRI, - n_time_steps, - initial_time, - final_time); - std::cout << " Dopri: error=" << solution.l2_norm() - << std::endl; - std::cout << " steps performed=" << n_steps << std::endl; - - n_steps = embedded_explicit_method(TimeStepping::FEHLBERG, - n_time_steps, - initial_time, - final_time); - std::cout << " Fehlberg: error=" << solution.l2_norm() - << std::endl; - std::cout << " steps performed=" << n_steps << std::endl; - - n_steps = embedded_explicit_method(TimeStepping::CASH_KARP, - n_time_steps, - initial_time, - final_time); - std::cout << " Cash-Karp: error=" << solution.l2_norm() - << std::endl; - std::cout << " steps performed=" << n_steps << std::endl; - } -} // namespace Step52 - - - -// @sect3{The main() function} -// -// The following main function is similar to previous examples -// and need not be commented on. -int main() -{ - try - { - Step52::Diffusion diffusion; - diffusion.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/examples/step-52.diff b/tests/examples/step-52.diff deleted file mode 100644 index cd5758a4e2..0000000000 --- a/tests/examples/step-52.diff +++ /dev/null @@ -1,74 +0,0 @@ -669c669 -< triangulation.refine_global(4); ---- -> triangulation.refine_global(3); -693c693 -< unsigned int n_steps = 0; ---- -> unsigned int n_steps [[maybe_unused]] = 0; -703,704c703 -< std::cout << " Forward Euler: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " Forward Euler: [no error reported]" << std::endl; -710,711c709 -< std::cout << " Third order Runge-Kutta: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " Third order Runge-Kutta: [no error reported]" << std::endl; -717,719c715 -< std::cout << " Fourth order Runge-Kutta: error=" << solution.l2_norm() -< << std::endl; -< std::cout << std::endl; ---- -> std::cout << " Fourth order Runge-Kutta: [no error reported]" << std::endl; -727,728c723 -< std::cout << " Backward Euler: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " Backward Euler: [no error reported]" << std::endl; -734,735c729 -< std::cout << " Implicit Midpoint: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " Implicit Midpoint: [no error reported]" << std::endl; -741,742c735 -< std::cout << " Crank-Nicolson: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " Crank-Nicolson: [no error reported]" << std::endl; -748,749c741 -< std::cout << " SDIRK: error=" << solution.l2_norm() -< << std::endl; ---- -> std::cout << " SDIRK: [no error reported]" << std::endl; -758,760c750 -< std::cout << " Heun-Euler: error=" << solution.l2_norm() -< << std::endl; -< std::cout << " steps performed=" << n_steps << std::endl; ---- -> std::cout << " Heun-Euler: [no error reported]" << std::endl; -766,768c756 -< std::cout << " Bogacki-Shampine: error=" << solution.l2_norm() -< << std::endl; -< std::cout << " steps performed=" << n_steps << std::endl; ---- -> std::cout << " Bogacki-Shampine: [no error reported]" << std::endl; -774,776c762 -< std::cout << " Dopri: error=" << solution.l2_norm() -< << std::endl; -< std::cout << " steps performed=" << n_steps << std::endl; ---- -> std::cout << " Dopri: [no error reported]" << std::endl; -782,784c768 -< std::cout << " Fehlberg: error=" << solution.l2_norm() -< << std::endl; -< std::cout << " steps performed=" << n_steps << std::endl; ---- -> std::cout << " Fehlberg: [no error reported]" << std::endl; -790,792c774 -< std::cout << " Cash-Karp: error=" << solution.l2_norm() -< << std::endl; -< std::cout << " steps performed=" << n_steps << std::endl; ---- -> std::cout << " Cash-Karp: [no error reported]" << std::endl; diff --git a/tests/examples/step-52.with_umfpack=true.output b/tests/examples/step-52.with_umfpack=true.output deleted file mode 100644 index 92c3651b6e..0000000000 --- a/tests/examples/step-52.with_umfpack=true.output +++ /dev/null @@ -1,16 +0,0 @@ -Explicit methods: - Forward Euler: [no error reported] - Third order Runge-Kutta: [no error reported] - Fourth order Runge-Kutta: [no error reported] -Implicit methods: - Backward Euler: [no error reported] - Implicit Midpoint: [no error reported] - Crank-Nicolson: [no error reported] - SDIRK: [no error reported] - -Embedded explicit methods: - Heun-Euler: [no error reported] - Bogacki-Shampine: [no error reported] - Dopri: [no error reported] - Fehlberg: [no error reported] - Cash-Karp: [no error reported] -- 2.39.5