From: David Wells Date: Sat, 4 May 2019 16:52:10 +0000 (-0400) Subject: Modernize step-24. X-Git-Tag: v9.1.0-rc1~132^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8006%2Fhead;p=dealii.git Modernize step-24. This PR introduces several small changes: 1. Remove trailing whitespace. 2. Remove unnecessary default constructor definitions. 3. Combine declarations and definitions of some simple functions. 4. Remove some unused headers. 5. Clean up the discussion of boundaries. 6. Switch to VTU output. --- diff --git a/examples/step-24/doc/intro.dox b/examples/step-24/doc/intro.dox index cbaf189e27..62dd817ff3 100644 --- a/examples/step-24/doc/intro.dox +++ b/examples/step-24/doc/intro.dox @@ -25,7 +25,7 @@ step_23 "step-23", where the wave equation was first introduced.

The problem

The temperature at a given location, neglecting thermal diffusion, can be -stated as +stated as @f[ \rho C_p \frac{\partial}{\partial t}T(t,\mathbf r) = H(t,\mathbf r) @@ -43,7 +43,7 @@ acoustically homogeneous medium can be described as follows: if $u$ is the vector-valued displacement, then tissue certainly reacts to changes in pressure by acceleration: @f[ -\rho \frac{\partial^2}{\partial t^2}u(t,\mathbf r) = +\rho \frac{\partial^2}{\partial t^2}u(t,\mathbf r) = -\nabla p(t,\mathbf r). @f] Furthermore, it contracts due to excess pressure and expands based on changes in temperature: @@ -59,7 +59,7 @@ time it takes a wave to cross the domain). In that case, the heating rate $H(t,\mathbf r)$ can be written as $H(t,\mathbf r) = a(\mathbf r)\delta(t)$ (where $a(\mathbf r)$ is a map of absorption strengths for microwave energy and $\delta(t)$ is the Dirac delta function), which together -with the first equation above will yield +with the first equation above will yield an instantaneous jump in the temperature $T(\mathbf r)$ at time $t=0$. Using this assumption, and taking all equations together, we can rewrite and combine the above as follows: @@ -67,14 +67,14 @@ rewrite and combine the above as follows: \Delta p-\frac{1}{c_0^2} \frac{\partial^2 p}{\partial t^2} = \lambda a(\mathbf r)\frac{d\delta(t)}{dt} @f] -where $\lambda = - \frac{\beta}{C_p}$. +where $\lambda = - \frac{\beta}{C_p}$. This somewhat strange equation with the derivative of a Dirac delta function on the right hand side can be rewritten as an initial value problem as follows: -@f{eqnarray*} +@f{eqnarray*} \Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial t^2} & = & -0 \\ -\bar{p}(0,\mathbf r) &=& c_0^2 \lambda a(\mathbf r) = b(\mathbf r) \\ +0 \\ +\bar{p}(0,\mathbf r) &=& c_0^2 \lambda a(\mathbf r) = b(\mathbf r) \\ \frac{\partial\bar{p}(0,\mathbf r)}{\partial t} &=& 0. @f} (A derivation of this transformation into an initial value problem is given at @@ -102,7 +102,7 @@ In general, this is a hard problem: Good absorbing boundary conditions are nonlinear and/or numerically very expensive. We therefore opt for a simple first order approximation to absorbing boundary conditions that reads @f[ -\frac{\partial\bar{p}}{\partial\mathbf n} = +\frac{\partial\bar{p}}{\partial\mathbf n} = -\frac{1}{c_0} \frac{\partial\bar{p}}{\partial t} @f] Here, $\frac{\partial\bar{p}}{\partial\mathbf n}$ is the normal derivative at @@ -114,15 +114,15 @@ condition, but it is one of the very few that are reasonably simple to implement As in step-23, one first introduces a second variable, which is defined as the derivative of the pressure potential: -@f[ -v = \frac{\partial\bar{p}}{\partial t} +@f[ +v = \frac{\partial\bar{p}}{\partial t} @f] With the second variable, one then transforms the forward problem into two separate equations: @f{eqnarray*} \bar{p}_{t} - v & = & 0 \\ -\Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f +\Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f @f} with initial conditions: @f{eqnarray*} @@ -139,16 +139,16 @@ introduced in step-23 is then: \left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_\Omega- \left(\theta v^{n}+(1-\theta)v^{n-1},\phi\right)_\Omega & = & 0 \\ -\left(\nabla((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla\phi\right)_\Omega- -\frac{1}{c_0}\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_{\partial\Omega} - -\frac{1}{c_0^2}\left(\frac{v^n-v^{n-1}}{k},\phi\right)_\Omega & = +\frac{1}{c_0}\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_{\partial\Omega} - +\frac{1}{c_0^2}\left(\frac{v^n-v^{n-1}}{k},\phi\right)_\Omega & = & \left(\theta f^{n}+(1-\theta)f^{n-1}, \phi\right)_\Omega, @f} where $\phi$ is an arbitrary test function, and where we have used the absorbing boundary condition to integrate by parts: -absorbing boundary conditions are incorporated into the weak form by using -@f[ +absorbing boundary conditions are incorporated into the weak form by using +@f[ \int_\Omega\varphi \, \Delta p\; dx = --\int_\Omega\nabla \varphi \cdot \nabla p dx + +-\int_\Omega\nabla \varphi \cdot \nabla p dx + \int_{\partial\Omega}\varphi \frac{\partial p}{\partial {\mathbf n}}ds. @f] @@ -157,7 +157,7 @@ functions, and get @f{eqnarray*} M\bar{p}^{n}-k \theta M v^n & = & M\bar{p}^{n-1}+k (1-\theta)Mv^{n-1},\\ -(-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = & +(-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = & (c_0^2k(1-\theta)A-c_0B)\bar{p}^{n-1}-Mv^{n-1}+c_0^2k(\theta F^{n}+(1-\theta)F^{n-1}). @f} The matrices $M$ and $A$ are here as in step-23, and the @@ -234,22 +234,22 @@ the indefinite integral with respect to time of both sides: @f[ \int^t \Delta p\; dt -\int^t \frac{1}{c_0^2} \frac{\partial^2 p}{\partial t^2} \; dt -= += \int^t \lambda a(\mathbf r)\frac{d\delta(t)}{dt} \;dt. @f] This immediately leads to the statement @f[ P(t,\mathbf r) - \frac{1}{c_0^2} \frac{\partial p}{\partial t} -= += \lambda a(\mathbf r) \delta(t), @f] where $P(t,\mathbf r)$ is such that $\frac{dP(t,\mathbf r)}{dt}=\Delta p$. Next, we form the (definite) integral over time from $t=-\epsilon$ to -$t=+\epsilon$ to find +$t=+\epsilon$ to find @f[ -\int_{-\epsilon}^{\epsilon} P(t,\mathbf r)\; dt +\int_{-\epsilon}^{\epsilon} P(t,\mathbf r)\; dt - \frac{1}{c_0^2} \left[ p(\epsilon,\mathbf r) - p(-\epsilon,\mathbf r) \right] -= += \int_{-\epsilon}^{\epsilon} \lambda a(\mathbf r) \delta(t) \; dt. @f] If we use the property of the delta function that $\int_{-\epsilon}^{\epsilon} @@ -257,13 +257,13 @@ If we use the property of the delta function that $\int_{-\epsilon}^{\epsilon} as we let $\epsilon$ go to zero that @f[ - \lim_{\epsilon\rightarrow 0}\frac{1}{c_0^2} \left[ p(\epsilon,\mathbf r) - p(-\epsilon,\mathbf r) \right] -= += \lambda a(\mathbf r). @f] In other words, using that $p(-\epsilon,\mathbf r)=0$, we retrieve the initial -condition +condition @f[ - \frac{1}{c_0^2} p(0,\mathbf r) + \frac{1}{c_0^2} p(0,\mathbf r) = \lambda a(\mathbf r). @f] @@ -303,7 +303,7 @@ see that in fact @f] Now, let $\epsilon\rightarrow 0$. Assuming that $P$ is a continuous function in -time, we see that +time, we see that @f[ P(\epsilon)-P(-\epsilon) \rightarrow 0, @f] @@ -319,4 +319,3 @@ Consequently, we obtain as the second initial condition that \frac{\partial p(0)}{\partial t} = 0, @f] completing the system of equations. - diff --git a/examples/step-24/doc/results.dox b/examples/step-24/doc/results.dox index a5db76bb9c..83742a45a1 100644 --- a/examples/step-24/doc/results.dox +++ b/examples/step-24/doc/results.dox @@ -1,14 +1,13 @@

Results

The program writes both graphical data for each time step as well as the -values evaluated at each detector location to disk. We then -draw them in plots. Experimental data were also collected for comparison. -Currently our experiments have only been done in two dimensions by -circularly scanning a single detector. The tissue sample here is a thin slice -in the X-Y plane (Z=0), and we assume that signals from other Z directions -won't contribute to the data. Consequently, we only have to compare -our experimental data with two dimensional simulated data. - +values evaluated at each detector location to disk. We then +draw them in plots. Experimental data were also collected for comparison. +Currently our experiments have only been done in two dimensions by +circularly scanning a single detector. The tissue sample here is a thin slice +in the $X-Y$ plane ($Z=0$), and we assume that signals from other $Z$ +directions won't contribute to the data. Consequently, we only have to compare +our experimental data with two dimensional simulated data.

One absorber

@@ -75,7 +74,7 @@ data: In our simulations, we see spurious signals behind the main wave that result from numerical artifacts. This problem can be alleviated by using finer mesh, resulting in the following plot: - + @@ -106,7 +105,7 @@ Note that in the experimental data, the first signal (i.e. the left-most dark line) results from absorption at the tissue boundary, and therefore reaches the detectors first and before any of the signals from the interior. This signal is also faintly visible at the end of the traces, around 30 $\mu s$, -which indicates that the signal travelled through the entire tissue to reach +which indicates that the signal traveled through the entire tissue to reach detectors at the other side, after all the signals originating from the interior have reached them. @@ -138,7 +137,3 @@ coarse mesh) is not isotropic either. The net result is that the signal fronts are not sinusoidal unless the mesh is sufficiently fine. The right image is a lot better in this respect, though artifacts in the form of trailing spurious waves can still be seen. - - - - diff --git a/examples/step-24/step-24.cc b/examples/step-24/step-24.cc index ed0622a618..7a29714488 100644 --- a/examples/step-24/step-24.cc +++ b/examples/step-24/step-24.cc @@ -117,10 +117,10 @@ namespace Step24 // @sect3{Equation data} // As usual, we have to define our initial values, boundary conditions, and - // right hand side functions. Except things are a bit simpler this time: we - // are to consider a problem that is driven by initial conditions, so there + // right hand side functions. Things are a bit simpler this time: we + // consider a problem that is driven by initial conditions, so there // is no right hand side function (though you could look up in step-23 to - // see how this can be done. Secondly, there are no boundary conditions: the + // see how this can be done). Secondly, there are no boundary conditions: the // entire boundary of the domain consists of absorbing boundary // conditions. That only leaves initial conditions, and there things are // simple too since for this particular application only nonzero initial @@ -139,12 +139,22 @@ namespace Step24 class InitialValuesP : public Function { public: - InitialValuesP() - : Function() - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + static const std::array sources{ + Source(Point(0, 0), 0.025), + Source(Point(-0.135, 0), 0.05), + Source(Point(0.17, 0), 0.03), + Source(Point(-0.25, 0), 0.02), + Source(Point(-0.05, -0.15), 0.015)}; + + for (const auto &source : sources) + if (p.distance(source.location) < source.radius) + return 1; + + return 0; + } private: struct Source @@ -160,25 +170,6 @@ namespace Step24 }; - template - double InitialValuesP::value(const Point &p, - const unsigned int /*component*/) const - { - static const Source sources[] = {Source(Point(0, 0), 0.025), - Source(Point(-0.135, 0), 0.05), - Source(Point(0.17, 0), 0.03), - Source(Point(-0.25, 0), 0.02), - Source(Point(-0.05, -0.15), 0.015)}; - static const unsigned int n_sources = sizeof(sources) / sizeof(sources[0]); - - for (unsigned int i = 0; i < n_sources; ++i) - if (p.distance(sources[i].location) < sources[i].radius) - return 1; - - return 0; - } - - // @sect3{Implementation of the TATForwardProblem class} // Let's start again with the constructor. Setting the member variables is @@ -231,11 +222,10 @@ namespace Step24 // The following system is pretty much what we've already done in step-23, // but with two important differences. First, we have to create a circular // (or spherical) mesh around the origin, with a radius of 1. This nothing - // new: we've done so before in step-6, step-10, and step-11, where we also - // explain how to attach a boundary object to a triangulation to be used - // whenever the triangulation needs to know where new boundary points lie - // when a cell is refined. Following this, the mesh is refined a number of - // times. + // new: we've done so before in step-6 and step-10, where we also explain + // how the PolarManifold or SphericalManifold object places new points on + // concentric circles when a cell is refined, which we will use here as + // well. // // One thing we had to make sure is that the time step satisfies the CFL // condition discussed in the introduction of step-23. Back in that program, @@ -350,12 +340,7 @@ namespace Step24 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()) for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) if (cell->at_boundary(f)) { @@ -417,6 +402,7 @@ namespace Step24 } + template void TATForwardProblem::solve_v() { @@ -446,9 +432,12 @@ namespace Step24 data_out.build_patches(); const std::string filename = - "solution-" + Utilities::int_to_string(timestep_number, 3) + ".gnuplot"; + "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu"; + DataOutBase::VtkFlags vtk_flags; + vtk_flags.compression_level = + DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; std::ofstream output(filename); - data_out.write_gnuplot(output); + data_out.write_vtu(output); } @@ -510,7 +499,6 @@ namespace Step24 solve_p(); - system_rhs_v = G2; laplace_matrix.vmult(tmp, solution_p); system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp); @@ -522,7 +510,6 @@ namespace Step24 output_results(); - detector_data << time; for (unsigned int i = 0; i < detector_locations.size(); ++i) detector_data << " " @@ -532,7 +519,6 @@ namespace Step24 << " "; detector_data << std::endl; - old_solution_p = solution_p; old_solution_v = solution_v; }