From: David Wells Date: Wed, 8 May 2019 04:09:53 +0000 (-0400) Subject: Modernize step-29. X-Git-Tag: v9.1.0-rc1~112^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b46d3d11119f21b04b15d5a18ae21dc3b52742d0;p=dealii.git Modernize step-29. 1. Inline some function definitions 2. Use ranged-for loops 3. Minor typography fixes (e.g., delete trailing whitespace) 4. Update console output --- diff --git a/examples/step-29/doc/intro.dox b/examples/step-29/doc/intro.dox index 452a8c20d3..35944314af 100644 --- a/examples/step-29/doc/intro.dox +++ b/examples/step-29/doc/intro.dox @@ -27,15 +27,15 @@ deal.II by using an FE_system object to stack two finite element fields representing real and imaginary parts. We also revisit the ParameterHandler class first used in step-19, which provides a convenient way for reading parameters from a configuration file at runtime -without the need to recompile the program code. +without the need to recompile the program code. The equations covered here fall into the class of vector-valued problems. A -toplevel overview of this topic can be found in the @ref vector_valued module. +top level overview of this topic can be found in the @ref vector_valued module.

Problem setting

-The original purpose of this program is to simulate the focussing properties +The original purpose of this program is to simulate the focusing properties of an ultrasound wave generated by a transducer lens with variable geometry. Recent applications in medical imaging use ultrasound waves not only for imaging purposes, but also to excite certain local effects in a @@ -43,10 +43,10 @@ material, like changes in optical properties, that can then be measured by other imaging techniques. A vital ingredient for these methods is the ability to focus the intensity of the ultrasound wave in a particular part of the material, ideally in a point, to be able to examine the properties of the -material at that particular location. +material at that particular location. To derive a model for this problem, we think of ultrasound as a pressure wave -governed by the wave equation: +governed by the wave equation: @f[ \frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0 @f] @@ -57,7 +57,7 @@ $\Gamma_2=\Gamma\setminus\Gamma_1$, with $\Gamma_1$ representing the transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to choose boundary conditions on $\Gamma_2$ in such a way that they imitate a larger domain). On $\Gamma_1$, the transducer generates a wave of constant -frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): +frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): @f[ U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1 @f] @@ -72,7 +72,7 @@ interested in. By plugging this form of the solution into the wave equation, we see that for $u$ we have @f{eqnarray*} -\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\ -u(x) &=& 1, \qquad x\in\Gamma_1. +u(x) &=& 1, \qquad x\in\Gamma_1. @f} For finding suitable conditions on $\Gamma_2$ that model an absorbing @@ -81,11 +81,11 @@ frequency ${\omega}$ traveling in direction $k\in {\mathrm{R}^2}$. In order for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of -$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation +$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation @f[ c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0. @f] -Hence, by enforcing the boundary condition +Hence, by enforcing the boundary condition @f[ c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2, @f] @@ -96,26 +96,26 @@ boundary condition will yield partial reflections, i.e. only parts of the wave will pass through the boundary as if it wasn't here whereas the remaining fraction of the wave will be reflected back into the domain. -If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$: +If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$: @f{eqnarray*} -\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\ c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\ -u &=& 1, \qquad x\in\Gamma_1. +u &=& 1, \qquad x\in\Gamma_1. @f} This is a Helmholtz equation (similar to the one in step-7, but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system: @f{eqnarray*} \left.\begin{array}{ccc} - -\omega^2 v - c^2\Delta v &=& 0 \quad\\ + -\omega^2 v - c^2\Delta v &=& 0 \quad\\ -\omega^2 w - c^2\Delta w &=& 0 \quad \end{array}\right\} &\;& x\in\Omega, \\ \left.\begin{array}{ccc} - c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ + c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ c (n\cdot\nabla w) + \omega\,v &=& 0 \quad \end{array}\right\} &\;& x\in\Gamma_2, \\ \left.\begin{array}{ccc} - v &=& 1 \quad\\ + v &=& 1 \quad\\ w &=& 0 \quad \end{array}\right\} &\;& x\in\Gamma_1. @f} @@ -126,19 +126,19 @@ For test functions $\phi,\psi$ with $\phi|_{\Gamma_1}=\psi|_{\Gamma_1}=0$, after -\omega^2 \langle \psi, w \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi, \nabla w \rangle_{\mathrm{L}^2(\Omega)} + c \omega \langle \psi, v \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0. @f} -We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions +We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions @f[ -v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. +v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. @f] Plugging into the variational form yields the equation system @f[ \renewcommand{\arraystretch}{2.0} \left.\begin{array}{ccc} -\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\ +\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\ \sum_{j=1}^n \left(-\omega^2 \langle \psi_i, \psi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \psi_i, \nabla \psi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\beta_j + \left(c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\alpha_j &=& 0 \end{array}\right\}\;\;\forall\; i =1,\ldots,n. @f] -In matrix notation: +In matrix notation: @f[ \renewcommand{\arraystretch}{2.0} \left( @@ -159,38 +159,38 @@ c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)} & -\omega^2 \langl \end{array} \right) @f] -(One should not be fooled by the right hand side being zero here, that is -because we haven't included the Dirichlet boundary data yet.) -Because of the alternating sign in the off-diagonal blocks, we can already +(One should not be fooled by the right hand side being zero here, that is +because we haven't included the Dirichlet boundary data yet.) +Because of the alternating sign in the off-diagonal blocks, we can already see that this system is non-symmetric, in fact it is even indefinite. Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be the same. However, we expect real and imaginary part of the solution to -have similar properties and will therefore indeed take $V_h=W_h$ in the -implementation, and also use the same basis functions $\phi_i = \psi_i$ for -both spaces. The reason for the notation using different symbols is just that -it allows us to distinguish between shape functions for $v$ and $w$, as this -distinction plays an important role in the implementation. +have similar properties and will therefore indeed take $V_h=W_h$ in the +implementation, and also use the same basis functions $\phi_i = \psi_i$ for +both spaces. The reason for the notation using different symbols is just that +it allows us to distinguish between shape functions for $v$ and $w$, as this +distinction plays an important role in the implementation.

The test case

-For the computations, we will consider wave propagation in the unit square, -with ultrasound generated by a transducer lens that is shaped like a segment -of the circle with center at $(0.5, d)$ and a +For the computations, we will consider wave propagation in the unit square, +with ultrasound generated by a transducer lens that is shaped like a segment +of the circle with center at $(0.5, d)$ and a radius slightly greater than $d$; this shape should lead to a focusing of the sound wave at the center of the circle. Varying $d$ changes the "focus" of the lens and affects the spatial distribution of the intensity of $u$, where our main -concern is how well $|u|=\sqrt{v^2+w^2}$ is focused. +concern is how well $|u|=\sqrt{v^2+w^2}$ is focused. In the program below, we will implement the complex-valued Helmholtz equations using the formulation with split real and imaginary parts. We will also discuss how to generate a domain that looks like a square with a slight bulge simulating the transducer (in the UltrasoundProblem::make_grid() function), and how to -generate graphical output that not only contains the solution components $v$ and +generate graphical output that not only contains the solution components $v$ and $w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in -UltrasoundProblem::output_results()). Finally, we use the -ParameterHandler class to easily read parameters like the focal distance $d$, -wave speed $c$, frequency $\omega$, and a number of other parameters from an -input file at run-time, rather than fixing those parameters in the source code +UltrasoundProblem::output_results()). Finally, we use the +ParameterHandler class to easily read parameters like the focal distance $d$, +wave speed $c$, frequency $\omega$, and a number of other parameters from an +input file at run-time, rather than fixing those parameters in the source code where we would have to re-compile every time we want to change parameters. diff --git a/examples/step-29/doc/results.dox b/examples/step-29/doc/results.dox index f463c8e8b9..0d0c7dad43 100644 --- a/examples/step-29/doc/results.dox +++ b/examples/step-29/doc/results.dox @@ -43,33 +43,39 @@ parameters here and therefore stick with their default values. Here's the console output of the program in debug mode: @code -examples/\step-29> make run -============================ Running \step-29 -DEAL::Generating grid... done (1.11607s) +> make run +[ 66%] Built target step-29 +[100%] Run step-29 with Debug configuration +DEAL::Generating grid... done (0.832682s) DEAL:: Number of active cells: 25600 -DEAL::Setting up system... done (1.10807s) +DEAL::Setting up system... done (0.614761s) DEAL:: Number of degrees of freedom: 51842 -DEAL::Assembling system matrix... done (6.50841s) -DEAL::Solving linear system... done (2.89218s) -DEAL::Generating output... done (4.52428s) +DEAL::Assembling system matrix... done (2.82536s) +DEAL::Solving linear system... done (2.27971s) +DEAL::Generating output... done (1.84418s) +[100%] Built target run @endcode (Of course, execution times will differ if you run the program locally.) The fact that most of the time is spent on assembling -the system matrix and generating output is due to the many assertion +the system matrix and generating output is due to the many assertions that need to be checked in debug mode. In release mode these parts of the program run much faster whereas solving the linear system is hardly sped up at all: @code -============================ Running \step-29 -DEAL::Generating grid... done (0.0280020s) +> make run +[ 66%] Built target step-29 +Scanning dependencies of target run +[100%] Run step-29 with Release configuration +DEAL::Generating grid... done (0.0144960s) DEAL:: Number of active cells: 25600 -DEAL::Setting up system... done (0.112007s) +DEAL::Setting up system... done (0.0356880s) DEAL:: Number of degrees of freedom: 51842 -DEAL::Assembling system matrix... done (0.160010s) -DEAL::Solving linear system... done (2.21614s) -DEAL::Generating output... done (0.880055s) +DEAL::Assembling system matrix... done (0.0436570s) +DEAL::Solving linear system... done (1.54733s) +DEAL::Generating output... done (0.720528s) +[100%] Built target run @endcode The graphical output of the program looks as follows: @@ -122,29 +128,41 @@ refined and which parts don't. Here are the run times for 5, 6, and 7 global refinements: @code -DEAL::Generating grid... done (0.0320020s) +> make run +[ 66%] Built target step-29 +[100%] Run step-29 with Release configuration +DEAL::Generating grid... done (0.0135260s) DEAL:: Number of active cells: 25600 -DEAL::Setting up system... done (0.104006s) +DEAL::Setting up system... done (0.0213910s) DEAL:: Number of degrees of freedom: 51842 -DEAL::Assembling system matrix... done (0.164011s) -DEAL::Solving linear system... done (2.23214s) -DEAL::Generating output... done (0.900056s) - -DEAL::Generating grid... done (0.132009s) +DEAL::Assembling system matrix... done (0.0414300s) +DEAL::Solving linear system... done (1.56621s) +DEAL::Generating output... done (0.729605s) +[100%] Built target run + +> make run +[ 66%] Built target step-29 +[100%] Run step-29 with Release configuration +DEAL::Generating grid... done (0.0668490s) DEAL:: Number of active cells: 102400 -DEAL::Setting up system... done (0.408025s) +DEAL::Setting up system... done (0.109694s) DEAL:: Number of degrees of freedom: 206082 -DEAL::Assembling system matrix... done (0.656041s) -DEAL::Solving linear system... done (14.8849s) -DEAL::Generating output... done (3.57222s) - -DEAL::Generating grid... done (0.504031s) +DEAL::Assembling system matrix... done (0.160784s) +DEAL::Solving linear system... done (7.86577s) +DEAL::Generating output... done (2.89320s) +[100%] Built target run + +> make run +[ 66%] Built target step-29 +[100%] Run step-29 with Release configuration +DEAL::Generating grid... done (0.293154s) DEAL:: Number of active cells: 409600 -DEAL::Setting up system... done (1.72011s) +DEAL::Setting up system... done (0.491301s) DEAL:: Number of degrees of freedom: 821762 -DEAL::Assembling system matrix... done (2.63216s) -DEAL::Solving linear system... done (117.811s) -DEAL::Generating output... done (15.1489s) +DEAL::Assembling system matrix... done (0.605386s) +DEAL::Solving linear system... done (45.1989s) +DEAL::Generating output... done (11.2292s) +[100%] Built target run @endcode Each time we refine the mesh once, so the number of cells and degrees diff --git a/examples/step-29/step-29.cc b/examples/step-29/step-29.cc index 805f55801f..db527110bc 100644 --- a/examples/step-29/step-29.cc +++ b/examples/step-29/step-29.cc @@ -20,8 +20,7 @@ // @sect3{Include files} -// The following header files are unchanged from step-7 and have been -// discussed before: +// The following header files have all been discussed before: #include #include @@ -103,38 +102,26 @@ namespace Step29 : Function(2) {} - virtual void vector_value(const Point &p, - Vector & values) const override; + virtual void vector_value(const Point & /*p*/, + Vector &values) const override + { + Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2)); + + values(0) = 1; + values(1) = 0; + } virtual void vector_value_list(const std::vector> &points, - std::vector> & value_list) const override; - }; - - - template - inline void - DirichletBoundaryValues::vector_value(const Point & /*p*/, - Vector &values) const - { - Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2)); - - values(0) = 1; - values(1) = 0; - } - - - template - void DirichletBoundaryValues::vector_value_list( - const std::vector> &points, - std::vector> & value_list) const - { - Assert(value_list.size() == points.size(), - ExcDimensionMismatch(value_list.size(), points.size())); + std::vector> & value_list) const override + { + Assert(value_list.size() == points.size(), + ExcDimensionMismatch(value_list.size(), points.size())); - for (unsigned int p = 0; p < points.size(); ++p) - DirichletBoundaryValues::vector_value(points[p], value_list[p]); - } + for (unsigned int p = 0; p < points.size(); ++p) + DirichletBoundaryValues::vector_value(points[p], value_list[p]); + } + }; // @sect3{The ParameterReader class} @@ -377,7 +364,6 @@ namespace Step29 { public: UltrasoundProblem(ParameterHandler &); - ~UltrasoundProblem(); void run(); private: @@ -411,13 +397,6 @@ namespace Step29 , fe(FE_Q(1), 2) {} - - template - UltrasoundProblem::~UltrasoundProblem() - { - dof_handler.clear(); - } - // @sect4{UltrasoundProblem::make_grid} // Here we setup the grid for our domain. As mentioned in the exposition, @@ -467,10 +446,7 @@ namespace Step29 // boundary indicator. GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1); - typename Triangulation::cell_iterator cell = triangulation.begin(), - endc = triangulation.end(); - - for (; cell != endc; ++cell) + for (auto &cell : triangulation.cell_iterators()) for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) if (cell->face(face)->at_boundary() && @@ -582,11 +558,7 @@ namespace Step29 FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); std::vector local_dof_indices(dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { // On each cell, we first need to reset the local contribution matrix // and request the FEValues object to compute the shape functions for