]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Modernize step-24. 8006/head
authorDavid Wells <drwells@email.unc.edu>
Sat, 4 May 2019 16:52:10 +0000 (12:52 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sun, 5 May 2019 17:23:51 +0000 (13:23 -0400)
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.

examples/step-24/doc/intro.dox
examples/step-24/doc/results.dox
examples/step-24/step-24.cc

index cbaf189e274363182d86de0792785ed75ce6ee6b..62dd817ff3afb0a1aff0dbcb38bbaf6a92a0e9ca 100644 (file)
@@ -25,7 +25,7 @@ step_23 "step-23", where the wave equation was first introduced.
 <h3>The problem</h3>
 
 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.
-
index a5db76bb9cd029fdbcebca52e47370715846f8f9..83742a45a1715fd1126bda9f0f5d8ead4bbd32e7 100644 (file)
@@ -1,14 +1,13 @@
 <h1>Results</h1>
 
 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.
 
 <h3> One absorber </h3>
 
@@ -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:
+
 <img src="https://www.dealii.org/images/steps/developer/step-24.one_s2.png" alt="">
 
 
@@ -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.
-
-
-
-
index ed0622a618ad28b6724ffffa3e3fbda0d2693020..7a29714488c62e684e54729a23e3fd68a0ac1204 100644 (file)
@@ -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<dim>
   {
   public:
-    InitialValuesP()
-      : Function<dim>()
-    {}
-
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      static const std::array<Source, 5> sources{
+        Source(Point<dim>(0, 0), 0.025),
+        Source(Point<dim>(-0.135, 0), 0.05),
+        Source(Point<dim>(0.17, 0), 0.03),
+        Source(Point<dim>(-0.25, 0), 0.02),
+        Source(Point<dim>(-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 <int dim>
-  double InitialValuesP<dim>::value(const Point<dim> &p,
-                                    const unsigned int /*component*/) const
-  {
-    static const Source       sources[] = {Source(Point<dim>(0, 0), 0.025),
-                                     Source(Point<dim>(-0.135, 0), 0.05),
-                                     Source(Point<dim>(0.17, 0), 0.03),
-                                     Source(Point<dim>(-0.25, 0), 0.02),
-                                     Source(Point<dim>(-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 <code>TATForwardProblem</code> 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<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
-
-
-      typename DoFHandler<dim>::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<dim>::faces_per_cell; ++f)
           if (cell->at_boundary(f))
             {
@@ -417,6 +402,7 @@ namespace Step24
   }
 
 
+
   template <int dim>
   void TATForwardProblem<dim>::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;
       }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.