step-24.h

Go to the documentation of this file.

 = 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
 *       {
 *         Source(const Point<dim> &l, const double r)
 *           : location(l)
 *           , radius(r)
 *         {}
 *   
 *         const Point<dim> location;
 *         const double     radius;
 *       };
 *     };
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeTATForwardProblemcodeclass"></a> 
 * <h3>Implementation of the <code>TATForwardProblem</code> class</h3>
 * 
 
 * 
 * Let's start again with the constructor. Setting the member variables is
 * straightforward. We use the acoustic wave speed of mineral oil (in
 * millimeters per microsecond, a common unit in experimental biomedical
 * imaging) since this is where many of the experiments we want to compare
 * the output with are made in. The Crank-Nicolson scheme is used again,
 * i.e. theta is set to 0.5. The time step is later selected to satisfy @f$k =
 * \frac hc@f$: here we initialize it to an invalid number.
 * 
 * @code
 *     template <int dim>
 *     TATForwardProblem<dim>::TATForwardProblem()
 *       : fe(1)
 *       , dof_handler(triangulation)
 *       , time_step(std::numeric_limits<double>::quiet_NaN())
 *       , time(time_step)
 *       , timestep_number(1)
 *       , theta(0.5)
 *       , wave_speed(1.437)
 *     {
 * @endcode
 * 
 * The second task in the constructor is to initialize the array that
 * holds the detector locations. The results of this program were compared
 * with experiments in which the step size of the detector spacing is 2.25
 * degree, corresponding to 160 detector locations. The radius of the
 * scanning circle is selected to be half way between the center and the
 * boundary to avoid that the remaining reflections from the imperfect
 * boundary condition spoils our numerical results.
 *     
 
 * 
 * The locations of the detectors are then calculated in clockwise
 * order. Note that the following of course only works if we are computing
 * in 2d, a condition that we guard with an assertion. If we later wanted
 * to run the same program in 3d, we would have to add code here for the
 * initialization of detector locations in 3d. Due to the assertion, there
 * is no way we can forget to do this.
 * 
 * @code
 *       Assert(dim == 2, ExcNotImplemented());
 *   
 *       const double detector_step_angle = 2.25;
 *       const double detector_radius     = 0.5;
 *   
 *       for (double detector_angle = 2 * numbers::PI; detector_angle >= 0;
 *            detector_angle -= detector_step_angle / 360 * 2 * numbers::PI)
 *         detector_locations.push_back(
 *           Point<dim>(std::cos(detector_angle), std::sin(detector_angle)) *
 *           detector_radius);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="TATForwardProblemsetup_system"></a> 
 * <h4>TATForwardProblem::setup_system</h4>
 * 
 
 * 
 * The following system is pretty much what we've already done in @ref step_23 "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 @ref step_6 "step-6" and @ref step_10 "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 @ref step_23 "step-23". Back in that program,
 * we ensured this by hand by setting a timestep that matches the mesh
 * width, but that was error prone because if we refined the mesh once more
 * we would also have to make sure the time step is changed. Here, we do
 * that automatically: we ask a library function for the minimal diameter of
 * any cell. Then we set @f$k=\frac h{c_0}@f$. The only problem is: what exactly
 * is @f$h@f$? The point is that there is really no good theory on this question
 * for the wave equation. It is known that for uniformly refined meshes
 * consisting of rectangles, @f$h@f$ is the minimal edge length. But for meshes
 * on general quadrilaterals, the exact relationship appears to be unknown,
 * i.e. it is unknown what properties of cells are relevant for the CFL
 * condition. The problem is that the CFL condition follows from knowledge
 * of the smallest eigenvalue of the Laplace matrix, and that can only be
 * computed analytically for simply structured meshes.
 *   
 
 * 
 * The upshot of all this is that we're not quite sure what exactly we
 * should take for @f$h@f$. The function GridTools::minimal_cell_diameter
 * computes the minimal diameter of all cells. If the cells were all squares
 * or cubes, then the minimal edge length would be the minimal diameter
 * divided by <code>std::sqrt(dim)</code>. We simply generalize this,
 * without theoretical justification, to the case of non-uniform meshes.
 *   
 
 * 
 * The only other significant change is that we need to build the boundary
 * mass matrix. We will comment on this further down below.
 * 
 * @code
 *     template <int dim>
 *     void TATForwardProblem<dim>::setup_system()
 *     {
 *       const Point<dim> center;
 *       GridGenerator::hyper_ball(triangulation, center, 1.);
 *       triangulation.refine_global(7);
 *   
 *       time_step = GridTools::minimal_cell_diameter(triangulation) / wave_speed /
 *                   std::sqrt(1. * dim);
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl;
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl
 *                 << std::endl;
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *       mass_matrix.reinit(sparsity_pattern);
 *       laplace_matrix.reinit(sparsity_pattern);
 *   
 *       MatrixCreator::create_mass_matrix(dof_handler,
 *                                         QGauss<dim>(fe.degree + 1),
 *                                         mass_matrix);
 *       MatrixCreator::create_laplace_matrix(dof_handler,
 *                                            QGauss<dim>(fe.degree + 1),
 *                                            laplace_matrix);
 *   
 * @endcode
 * 
 * The second difference, as mentioned, to @ref step_23 "step-23" is that we need to
 * build the boundary mass matrix that grew out of the absorbing boundary
 * conditions.
 *     
 
 * 
 * A first observation would be that this matrix is much sparser than the
 * regular mass matrix, since none of the shape functions with purely
 * interior support contribute to this matrix. We could therefore
 * optimize the storage pattern to this situation and build up a second
 * sparsity pattern that only contains the nonzero entries that we
 * need. There is a trade-off to make here: first, we would have to have a
 * second sparsity pattern object, so that costs memory. Secondly, the
 * matrix attached to this sparsity pattern is going to be smaller and
 * therefore requires less memory; it would also be faster to perform
 * matrix-vector multiplications with it. The final argument, however, is
 * the one that tips the scale: we are not primarily interested in
 * performing matrix-vector with the boundary matrix alone (though we need
 * to do that for the right hand side vector once per time step), but
 * mostly wish to add it up to the other matrices used in the first of the
 * two equations since this is the one that is going to be multiplied with
 * once per iteration of the CG method, i.e. significantly more often. It
 * is now the case that the SparseMatrix::add class allows to add one
 * matrix to another, but only if they use the same sparsity pattern (the
 * reason being that we can't add nonzero entries to a matrix after the
 * sparsity pattern has been created, so we simply require that the two
 * matrices have the same sparsity pattern).
 *     
 
 * 
 * So let's go with that:
 * 
 * @code
 *       boundary_matrix.reinit(sparsity_pattern);
 *   
 * @endcode
 * 
 * The second thing to do is to actually build the matrix. Here, we need
 * to integrate over faces of cells, so first we need a quadrature object
 * that works on <code>dim-1</code> dimensional objects. Secondly, the
 * FEFaceValues variant of FEValues that works on faces, as its name
 * suggest. And finally, the other variables that are part of the assembly
 * machinery. All of this we put between curly braces to limit the scope
 * of these variables to where we actually need them.
 *     
 
 * 
 * The actual act of assembling the matrix is then fairly straightforward:
 * we loop over all cells, over all faces of each of these cells, and then
 * do something only if that particular face is at the boundary of the
 * domain. Like this:
 * 
 * @code
 *       {
 *         const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
 *         FEFaceValues<dim>     fe_values(fe,
 *                                     quadrature_formula,
 *                                     update_values | update_JxW_values);
 *   
 *         const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *         const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *         FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *   
 *         std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *         for (const auto &cell : dof_handler.active_cell_iterators())
 *           for (const auto &face : cell->face_iterators())
 *             if (face->at_boundary())
 *               {
 *                 cell_matrix = 0;
 *   
 *                 fe_values.reinit(cell, face);
 *   
 *                 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) += (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);
 *                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                     boundary_matrix.add(local_dof_indices[i],
 *                                         local_dof_indices[j],
 *                                         cell_matrix(i, j));
 *               }
 *       }
 *   
 *       system_matrix.copy_from(mass_matrix);
 *       system_matrix.add(time_step * time_step * theta * theta * wave_speed *
 *                           wave_speed,
 *                         laplace_matrix);
 *       system_matrix.add(wave_speed * theta * time_step, boundary_matrix);
 *   
 *   
 *       solution_p.reinit(dof_handler.n_dofs());
 *       old_solution_p.reinit(dof_handler.n_dofs());
 *       system_rhs_p.reinit(dof_handler.n_dofs());
 *   
 *       solution_v.reinit(dof_handler.n_dofs());
 *       old_solution_v.reinit(dof_handler.n_dofs());
 *       system_rhs_v.reinit(dof_handler.n_dofs());
 *   
 *       constraints.close();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="TATForwardProblemsolve_pandTATForwardProblemsolve_v"></a> 
 * <h4>TATForwardProblem::solve_p and TATForwardProblem::solve_v</h4>
 * 
 
 * 
 * The following two functions, solving the linear systems for the pressure
 * and the velocity variable, are taken pretty much verbatim (with the
 * exception of the change of name from @f$u@f$ to @f$p@f$ of the primary variable)
 * from @ref step_23 "step-23":
 * 
 * @code
 *     template <int dim>
 *     void TATForwardProblem<dim>::solve_p()
 *     {
 *       SolverControl solver_control(1000, 1e-8 * system_rhs_p.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       cg.solve(system_matrix, solution_p, system_rhs_p, PreconditionIdentity());
 *   
 *       std::cout << "   p-equation: " << solver_control.last_step()
 *                 << " CG iterations." << std::endl;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void TATForwardProblem<dim>::solve_v()
 *     {
 *       SolverControl solver_control(1000, 1e-8 * system_rhs_v.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       cg.solve(mass_matrix, solution_v, system_rhs_v, PreconditionIdentity());
 *   
 *       std::cout << "   v-equation: " << solver_control.last_step()
 *                 << " CG iterations." << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="TATForwardProblemoutput_results"></a> 
 * <h4>TATForwardProblem::output_results</h4>
 * 
 
 * 
 * The same holds here: the function is from @ref step_23 "step-23".
 * 
 * @code
 *     template <int dim>
 *     void TATForwardProblem<dim>::output_results() const
 *     {
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution_p, "P");
 *       data_out.add_data_vector(solution_v, "V");
 *   
 *       data_out.build_patches();
 *   
 *       const std::string filename =
 *         "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
 *       DataOutBase::VtkFlags vtk_flags;
 *       vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="TATForwardProblemrun"></a> 
 * <h4>TATForwardProblem::run</h4>
 * 
 
 * 
 * This function that does most of the work is pretty much again like in
 * @ref step_23 "step-23", though we make things a bit clearer by using the vectors G1 and
 * G2 mentioned in the introduction. Compared to the overall memory
 * consumption of the program, the introduction of a few temporary vectors
 * isn't doing much harm.
 *   
 
 * 
 * The only changes to this function are: first, that we do not have to
 * project initial values for the velocity @f$v@f$, since we know that it is
 * zero. And second that we evaluate the solution at the detector locations
 * computed in the constructor. This is done using the
 * VectorTools::point_value function. These values are then written to a
 * file that we open at the beginning of the function.
 * 
 * @code
 *     template <int dim>
 *     void TATForwardProblem<dim>::run()
 *     {
 *       setup_system();
 *   
 *       VectorTools::project(dof_handler,
 *                            constraints,
 *                            QGauss<dim>(fe.degree + 1),
 *                            InitialValuesP<dim>(),
 *                            old_solution_p);
 *       old_solution_v = 0;
 *   
 *   
 *       std::ofstream detector_data("detectors.dat");
 *   
 *       Vector<double> tmp(solution_p.size());
 *       Vector<double> G1(solution_p.size());
 *       Vector<double> G2(solution_v.size());
 *   
 *       const double end_time = 0.7;
 *       for (time = time_step; time <= end_time;
 *            time += time_step, ++timestep_number)
 *         {
 *           std::cout << std::endl;
 *           std::cout << "time_step " << timestep_number << " @ t=" << time
 *                     << std::endl;
 *   
 *           mass_matrix.vmult(G1, old_solution_p);
 *           mass_matrix.vmult(tmp, old_solution_v);
 *           G1.add(time_step * (1 - theta), tmp);
 *   
 *           mass_matrix.vmult(G2, old_solution_v);
 *           laplace_matrix.vmult(tmp, old_solution_p);
 *           G2.add(-wave_speed * wave_speed * time_step * (1 - theta), tmp);
 *   
 *           boundary_matrix.vmult(tmp, old_solution_p);
 *           G2.add(wave_speed, tmp);
 *   
 *           system_rhs_p = G1;
 *           system_rhs_p.add(time_step * theta, G2);
 *   
 *           solve_p();
 *   
 *           system_rhs_v = G2;
 *           laplace_matrix.vmult(tmp, solution_p);
 *           system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
 *   
 *           boundary_matrix.vmult(tmp, solution_p);
 *           system_rhs_v.add(-wave_speed, tmp);
 *   
 *           solve_v();
 *   
 *           output_results();
 *   
 *           detector_data << time;
 *           for (unsigned int i = 0; i < detector_locations.size(); ++i)
 *             detector_data << ' '
 *                           << VectorTools::point_value(dof_handler,
 *                                                       solution_p,
 *                                                       detector_locations[i])
 *                           << ' ';
 *           detector_data << std::endl;
 *   
 *           old_solution_p = solution_p;
 *           old_solution_v = solution_v;
 *         }
 *     }
 *   } // namespace Step24
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * What remains is the main function of the program. There is nothing here
 * that hasn't been shown in several of the previous programs:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step24;
 *   
 *         TATForwardProblem<2> forward_problem_solver;
 *         forward_problem_solver.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><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 @f$X-Y@f$ plane (@f$Z=0@f$), and we assume that signals from other @f$Z@f$
directions won't contribute to the data. Consequently, we only have to compare
our experimental data with two dimensional simulated data.
 
<a name="Oneabsorber"></a><h3> One absorber </h3>
 
 
This movie shows the thermoacoustic waves generated by a single small absorber
propagating in the medium (in our simulation, we assume the medium is mineral
oil, which has a acoustic speed of 1.437 @f$\frac{mm}{\mu s}@f$):
 
<img src="https://www.dealii.org/images/steps/developer/step-24.one_movie.gif" alt="">
 
For a single absorber, we of course have to change the
<code>InitialValuesP</code> class accordingly.
 
Next, let us compare experimental and computational results. The visualization
uses a technique long used in seismology, where the data of each detector is
plotted all in one graph. The way this is done is by offsetting each
detector's signal a bit compared to the previous one. For example, here is a
plot of the first four detectors (from bottom to top, with time in
microseconds running from left to right) using the source setup used in the
program, to make things a bit more interesting compared to the present case of
only a single source:
 
<img src="https://www.dealii.org/images/steps/developer/step-24.traces.png" alt="">
 
One thing that can be seen, for example, is that the arrival of the second and
fourth signals shifts to earlier times for greater detector numbers (i.e. the
topmost ones), but not the first and the third; this can be interpreted to
mean that the origin of these signals must be closer to the latter detectors
than to the former ones.
 
If we stack not only 4, but all 160 detectors in one graph, the individual
lines blur, but where they run together they create a pattern of darker or
lighter grayscales.  The following two figures show the results obtained at
the detector locations stacked in that way. The left figure is obtained from
experiments, and the right is the simulated data.
In the experiment, a single small strong absorber was embedded in
weaker absorbing tissue:
 
<table width="100%">
<tr>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.one.png" alt="">
</td>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.one_s.png" alt="">
</td>
</tr>
</table>
 
It is obvious that the source location is closer to the detectors at angle
@f$180^\circ@f$. All the other signals that can be seen in the experimental data
result from the fact that there are weak absorbers also in the rest of the
tissue, which surrounds the signals generated by the small strong absorber in
the center. On the other hand, in the simulated data, we only simulate the
small strong absorber.
 
In reality, detectors have limited bandwidth. The thermoacoustic waves passing
through the detector will therefore be filtered. By using a high-pass filter
(implemented in MATLAB and run against the data file produced by this program),
the simulated results can be made to look closer to the experimental
data:
 
<img src="https://www.dealii.org/images/steps/developer/step-24.one_sf.png" alt="">
 
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="">
 
 
 
<a name="Multipleabsorbers"></a><h3>Multiple absorbers</h3>
 
 
To further verify the program, we will also show simulation results for
multiple absorbers. This corresponds to the case that is actually implemented
in the program. The following movie shows the propagation of the generated
thermoacoustic waves in the medium by multiple absorbers:
 
<img src="https://www.dealii.org/images/steps/developer/step-24.multi_movie.gif" alt="">
 
Experimental data and our simulated data are compared in the following two
figures:
<table width="100%">
<tr>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.multi.png" alt="">
</td>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s.png" alt="">
</td>
</tr>
</table>
 
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 @f$\mu s@f$,
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.
 
As before, the numerical result better matches experimental ones by applying a
bandwidth filter that matches the actual behavior of detectors (left) and by
choosing a finer mesh (right):
 
<table width="100%">
<tr>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.multi_sf.png" alt="">
</td>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s2.png" alt="">
</td>
</tr>
</table>
 
One of the important differences between the left and the right figure is that
the curves look much less "angular" at the right. The angularity comes from
the fact that while waves in the continuous equation travel equally fast in
all directions, this isn't the case after discretization: there, waves that
travel diagonal to cells move at slightly different speeds to those that move
parallel to mesh lines. This anisotropy leads to wave fronts that aren't
perfectly circular (and would produce sinusoidal signals in the stacked
plots), but are bulged out in certain directions. To make things worse, the
circular mesh we use (see for example @ref step_6 "step-6" for a view of the
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.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-24.cc"
*/