531 *
static const std::array<Source, 5> sources{
538 *
for (
const auto &source : sources)
539 * if (p.distance(source.location) < source.radius)
554 *
const double radius;
562 * <a name=
"step_24-ImplementationofthecodeTATForwardProblemcodeclass"></a>
563 * <h3>Implementation of the <code>TATForwardProblem</code>
class</h3>
567 * Let
's start again with the constructor. Setting the member variables is
568 * straightforward. We use the acoustic wave speed of mineral oil (in
569 * millimeters per microsecond, a common unit in experimental biomedical
570 * imaging) since this is where many of the experiments we want to compare
571 * the output with are made in. The Crank-Nicolson scheme is used again,
572 * i.e. theta is set to 0.5. The time step is later selected to satisfy @f$k =
573 * \frac hc@f$: here we initialize it to an invalid number.
577 * TATForwardProblem<dim>::TATForwardProblem()
579 * , dof_handler(triangulation)
580 * , time_step(std::numeric_limits<double>::quiet_NaN())
582 * , timestep_number(1)
584 * , wave_speed(1.437)
588 * The second task in the constructor is to initialize the array that
589 * holds the detector locations. The results of this program were compared
590 * with experiments in which the step size of the detector spacing is 2.25
591 * degree, corresponding to 160 detector locations. The radius of the
592 * scanning circle is selected to be half way between the center and the
593 * boundary to avoid that the remaining reflections from the imperfect
594 * boundary condition spoils our numerical results.
598 * The locations of the detectors are then calculated in clockwise
599 * order. Note that the following of course only works if we are computing
600 * in 2d, a condition that we guard with an assertion. If we later wanted
601 * to run the same program in 3d, we would have to add code here for the
602 * initialization of detector locations in 3d. Due to the assertion, there
603 * is no way we can forget to do this.
606 * Assert(dim == 2, ExcNotImplemented());
608 * const double detector_step_angle = 2.25;
609 * const double detector_radius = 0.5;
611 * for (double detector_angle = 2 * numbers::PI; detector_angle >= 0;
612 * detector_angle -= detector_step_angle / 360 * 2 * numbers::PI)
613 * detector_locations.push_back(
614 * Point<dim>(std::cos(detector_angle), std::sin(detector_angle)) *
623 * <a name="step_24-TATForwardProblemsetup_system"></a>
624 * <h4>TATForwardProblem::setup_system</h4>
628 * The following system is pretty much what we've already done in @ref step_23
"step-23",
629 * but with two important differences. First, we have to create a circular
630 * (or spherical) mesh around the origin, with a radius of 1. This nothing
631 *
new: we
've done so before in @ref step_6 "step-6" and @ref step_10 "step-10", where we also explain
632 * how the PolarManifold or SphericalManifold object places new points on
633 * concentric circles when a cell is refined, which we will use here as
638 * One thing we had to make sure is that the time step satisfies the CFL
639 * condition discussed in the introduction of @ref step_23 "step-23". Back in that program,
640 * we ensured this by hand by setting a timestep that matches the mesh
641 * width, but that was error prone because if we refined the mesh once more
642 * we would also have to make sure the time step is changed. Here, we do
643 * that automatically: we ask a library function for the minimal diameter of
644 * any cell. Then we set @f$k=\frac h{c_0}@f$. The only problem is: what exactly
645 * is @f$h@f$? The point is that there is really no good theory on this question
646 * for the wave equation. It is known that for uniformly refined meshes
647 * consisting of rectangles, @f$h@f$ is the minimal edge length. But for meshes
648 * on general quadrilaterals, the exact relationship appears to be unknown,
649 * i.e. it is unknown what properties of cells are relevant for the CFL
650 * condition. The problem is that the CFL condition follows from knowledge
651 * of the smallest eigenvalue of the Laplace matrix, and that can only be
652 * computed analytically for simply structured meshes.
656 * The upshot of all this is that we're not quite sure what exactly we
658 * computes the minimal
diameter of all cells. If the cells were all squares
659 * or cubes, then the minimal edge length would be the minimal
diameter
660 * divided by <code>
std::
sqrt(dim)</code>. We simply generalize this,
661 * without theoretical justification, to the case of non-uniform meshes.
665 * The only other significant change is that we need to build the boundary
666 * mass
matrix. We will comment on this further down below.
670 * void TATForwardProblem<dim>::setup_system()
679 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
682 * dof_handler.distribute_dofs(fe);
684 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
690 * sparsity_pattern.copy_from(dsp);
692 * system_matrix.reinit(sparsity_pattern);
694 * laplace_matrix.reinit(sparsity_pattern);
705 * The
second difference, as mentioned, to @ref step_23
"step-23" is that we need to
706 * build the boundary mass
matrix that grew out of the absorbing boundary
711 * A
first observation would be that
this matrix is much sparser than the
713 * interior support contribute to
this matrix. We could therefore
714 * optimize the storage pattern to
this situation and build up a
second
715 * sparsity pattern that only contains the
nonzero entries that we
716 * need. There is a trade-off to make here:
first, we would have to have a
717 *
second sparsity pattern object, so that costs memory. Secondly, the
718 *
matrix attached to
this sparsity pattern is going to be smaller and
719 * therefore
requires less memory; it would also be faster to perform
720 *
matrix-vector multiplications with it. The
final argument, however, is
721 * the one that tips the
scale: we are not primarily interested in
722 * performing
matrix-vector with the boundary
matrix alone (though we need
723 * to
do that
for the right hand side vector once per time step), but
724 * mostly wish to add it up to the other matrices used in the
first of the
725 * two equations since
this is the one that is going to be multiplied with
726 * once per iteration of the CG method, i.e. significantly more often. It
728 *
matrix to another, but only
if they use the same sparsity pattern (the
729 * reason being that we can
't add nonzero entries to a matrix after the
730 * sparsity pattern has been created, so we simply require that the two
731 * matrices have the same sparsity pattern).
735 * So let's go with that:
738 * boundary_matrix.
reinit(sparsity_pattern);
742 * The
second thing to
do is to actually build the
matrix. Here, we need
743 * to integrate over faces of cells, so
first we need a quadrature
object
744 * that works on <code>dim-1</code> dimensional objects. Secondly, the
746 * suggest. And
finally, the other variables that are part of the assembly
747 * machinery. All of
this we put between curly braces to limit the scope
748 * of these variables to where we actually need them.
752 * The actual act of assembling the
matrix is then fairly straightforward:
753 * we
loop over all cells, over all faces of each of these cells, and then
754 *
do something only
if that particular face is at the boundary of the
759 *
const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
761 * quadrature_formula,
764 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
765 *
const unsigned int n_q_points = quadrature_formula.size();
769 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
771 *
for (
const auto &cell : dof_handler.active_cell_iterators())
772 * for (const auto &face : cell->face_iterators())
773 * if (face->at_boundary())
777 * fe_values.reinit(cell, face);
779 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
780 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
781 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
782 *
cell_matrix(i, j) += (fe_values.shape_value(i, q_point) *
783 * fe_values.shape_value(j, q_point) *
784 * fe_values.JxW(q_point));
786 * cell->get_dof_indices(local_dof_indices);
787 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
788 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
789 * boundary_matrix.add(local_dof_indices[i],
790 * local_dof_indices[j],
795 * system_matrix.copy_from(mass_matrix);
796 * system_matrix.add(time_step * time_step * theta * theta * wave_speed *
799 * system_matrix.add(wave_speed * theta * time_step, boundary_matrix);
802 * solution_p.reinit(dof_handler.n_dofs());
803 * old_solution_p.reinit(dof_handler.n_dofs());
804 * system_rhs_p.reinit(dof_handler.n_dofs());
806 * solution_v.reinit(dof_handler.n_dofs());
807 * old_solution_v.reinit(dof_handler.n_dofs());
808 * system_rhs_v.reinit(dof_handler.n_dofs());
810 * constraints.close();
817 * <a name=
"step_24-TATForwardProblemsolve_pandTATForwardProblemsolve_v"></a>
818 * <h4>TATForwardProblem::solve_p and TATForwardProblem::solve_v</h4>
822 * The following two
functions, solving the linear systems
for the pressure
823 * and the velocity variable, are taken pretty much verbatim (with the
824 * exception of the change of name from @f$u@f$ to @f$p@f$ of the primary variable)
825 * from @ref step_23
"step-23":
829 *
void TATForwardProblem<dim>::solve_p()
831 *
SolverControl solver_control(1000, 1e-8 * system_rhs_p.l2_norm());
836 * std::cout <<
" p-equation: " << solver_control.last_step()
837 * <<
" CG iterations." << std::endl;
843 *
void TATForwardProblem<dim>::solve_v()
845 *
SolverControl solver_control(1000, 1e-8 * system_rhs_v.l2_norm());
850 * std::cout <<
" v-equation: " << solver_control.last_step()
851 * <<
" CG iterations." << std::endl;
859 * <a name=
"step_24-TATForwardProblemoutput_results"></a>
860 * <h4>TATForwardProblem::output_results</h4>
864 * The same holds here: the function is from @ref step_23
"step-23".
868 *
void TATForwardProblem<dim>::output_results() const
873 * data_out.add_data_vector(solution_p,
"P");
874 * data_out.add_data_vector(solution_v,
"V");
876 * data_out.build_patches();
878 *
const std::string filename =
882 * std::ofstream output(filename);
883 * data_out.write_vtu(output);
891 * <a name=
"step_24-TATForwardProblemrun"></a>
892 * <h4>TATForwardProblem::run</h4>
896 * This function that does most of the work is pretty much again like in
897 * @ref step_23
"step-23", though we make things a bit clearer by
using the vectors G1 and
898 * G2 mentioned in the introduction. Compared to the overall memory
899 * consumption of the program, the introduction of a few temporary vectors
900 * isn
't doing much harm.
904 * The only changes to this function are: first, that we do not have to
905 * project initial values for the velocity @f$v@f$, since we know that it is
906 * zero. And second that we evaluate the solution at the detector locations
907 * computed in the constructor. This is done using the
908 * VectorTools::point_value function. These values are then written to a
909 * file that we open at the beginning of the function.
913 * void TATForwardProblem<dim>::run()
917 * VectorTools::project(dof_handler,
919 * QGauss<dim>(fe.degree + 1),
920 * InitialValuesP<dim>(),
922 * old_solution_v = 0;
925 * std::ofstream detector_data("detectors.dat");
927 * Vector<double> tmp(solution_p.size());
928 * Vector<double> G1(solution_p.size());
929 * Vector<double> G2(solution_v.size());
931 * const double end_time = 0.7;
932 * for (time = time_step; time <= end_time;
933 * time += time_step, ++timestep_number)
935 * std::cout << std::endl;
936 * std::cout << "time_step " << timestep_number << " @ t=" << time
939 * mass_matrix.vmult(G1, old_solution_p);
940 * mass_matrix.vmult(tmp, old_solution_v);
941 * G1.add(time_step * (1 - theta), tmp);
943 * mass_matrix.vmult(G2, old_solution_v);
944 * laplace_matrix.vmult(tmp, old_solution_p);
945 * G2.add(-wave_speed * wave_speed * time_step * (1 - theta), tmp);
947 * boundary_matrix.vmult(tmp, old_solution_p);
948 * G2.add(wave_speed, tmp);
951 * system_rhs_p.add(time_step * theta, G2);
956 * laplace_matrix.vmult(tmp, solution_p);
957 * system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
959 * boundary_matrix.vmult(tmp, solution_p);
960 * system_rhs_v.add(-wave_speed, tmp);
966 * detector_data << time;
967 * for (unsigned int i = 0; i < detector_locations.size(); ++i)
968 * detector_data << ' '
969 * << VectorTools::point_value(dof_handler,
971 * detector_locations[i])
973 * detector_data << std::endl;
975 * old_solution_p = solution_p;
976 * old_solution_v = solution_v;
979 * } // namespace Step24
986 * <a name="step_24-Thecodemaincodefunction"></a>
987 * <h3>The <code>main</code> function</h3>
991 * What remains is the main function of the program. There is nothing here
992 * that hasn't been shown in several of the previous programs:
999 *
using namespace Step24;
1001 * TATForwardProblem<2> forward_problem_solver;
1002 * forward_problem_solver.run();
1004 *
catch (std::exception &exc)
1006 * std::cerr << std::endl
1008 * <<
"----------------------------------------------------"
1010 * std::cerr <<
"Exception on processing: " << std::endl
1011 * << exc.what() << std::endl
1012 * <<
"Aborting!" << std::endl
1013 * <<
"----------------------------------------------------"
1020 * std::cerr << std::endl
1022 * <<
"----------------------------------------------------"
1024 * std::cerr <<
"Unknown exception!" << std::endl
1025 * <<
"Aborting!" << std::endl
1026 * <<
"----------------------------------------------------"
1034<a name=
"step_24-Results"></a><h1>Results</h1>
1037The program writes both graphical data
for each time step as well as the
1038values evaluated at each detector location to disk. We then
1039draw them in plots. Experimental data were also collected
for comparison.
1040Currently our experiments have only been done in two dimensions by
1041circularly scanning a single detector. The tissue sample here is a thin slice
1042in the @f$X-Y@f$ plane (@f$Z=0@f$), and we assume that signals from other @f$Z@f$
1043directions won
't contribute to the data. Consequently, we only have to compare
1044our experimental data with two dimensional simulated data.
1046<a name="step_24-Oneabsorber"></a><h3> One absorber </h3>
1049This movie shows the thermoacoustic waves generated by a single small absorber
1050propagating in the medium (in our simulation, we assume the medium is mineral
1051oil, which has a acoustic speed of 1.437 @f$\frac{mm}{\mu s}@f$):
1053<img src="https://www.dealii.org/images/steps/developer/step-24.one_movie.gif" alt="">
1055For a single absorber, we of course have to change the
1056<code>InitialValuesP</code> class accordingly.
1058Next, let us compare experimental and computational results. The visualization
1059uses a technique long used in seismology, where the data of each detector is
1060plotted all in one graph. The way this is done is by offsetting each
1061detector's signal a bit compared to the previous one. For example, here is a
1062plot of the
first four detectors (from bottom to top, with time in
1063microseconds running from left to right)
using the source setup used in the
1064program, to make things a bit more interesting compared to the present
case of
1065only a single source:
1067<img src=
"https://www.dealii.org/images/steps/developer/step-24.traces.png" alt=
"">
1069One thing that can be seen,
for example, is that the arrival of the
second and
1070fourth signals shifts to earlier times
for greater detector
numbers (i.e. the
1071topmost ones), but not the
first and the third;
this can be interpreted to
1072mean that the origin of these signals must be closer to the latter detectors
1073than to the former ones.
1075If we stack not only 4, but all 160 detectors in one graph, the individual
1076lines blur, but where they
run together they create a pattern of darker or
1077lighter grayscales. The following two figures show the results obtained at
1078the detector locations stacked in that way. The left figure is obtained from
1079experiments, and the right is the simulated data.
1080In the experiment, a single small strong absorber was embedded in
1081weaker absorbing tissue:
1086<img src=
"https://www.dealii.org/images/steps/developer/step-24.one.png" alt=
"">
1089<img src=
"https://www.dealii.org/images/steps/developer/step-24.one_s.png" alt=
"">
1094It is obvious that the source location is closer to the detectors at
angle
1095@f$180^\circ@f$. All the other signals that can be seen in the experimental data
1096result from the fact that there are weak absorbers also in the rest of the
1097tissue, which surrounds the signals generated by the small strong absorber in
1098the
center. On the other hand, in the simulated data, we only simulate the
1099small strong absorber.
1101In reality, detectors have limited bandwidth. The thermoacoustic waves passing
1102through the detector will therefore be filtered. By
using a high-pass filter
1103(implemented in MATLAB and run against the data file produced by
this program),
1104the simulated results can be made to look closer to the experimental
1107<img src=
"https://www.dealii.org/images/steps/developer/step-24.one_sf.png" alt=
"">
1109In our simulations, we see spurious signals behind the main wave that
1110result from numerical artifacts. This problem can be alleviated by
using finer
1111mesh, resulting in the following plot:
1113<img src=
"https://www.dealii.org/images/steps/developer/step-24.one_s2.png" alt=
"">
1117<a name=
"step_24-Multipleabsorbers"></a><h3>Multiple absorbers</h3>
1120To further verify the program, we will also show simulation results
for
1121multiple absorbers. This corresponds to the
case that is actually implemented
1122in the program. The following movie shows the propagation of the generated
1123thermoacoustic waves in the medium by multiple absorbers:
1125<img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_movie.gif" alt=
"">
1127Experimental data and our simulated data are compared in the following two
1132<img src=
"https://www.dealii.org/images/steps/developer/step-24.multi.png" alt=
"">
1135<img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_s.png" alt=
"">
1140Note that in the experimental data, the
first signal (i.e. the left-most dark
1141line) results from absorption at the tissue boundary, and therefore reaches
1142the detectors
first and before any of the signals from the interior. This
1143signal is also faintly visible at the
end of the traces, around 30 @f$\mu s@f$,
1144which indicates that the signal traveled through the entire tissue to reach
1145detectors at the other side, after all the signals originating from the
1146interior have reached them.
1148As before, the numerical result better matches experimental ones by applying a
1149bandwidth filter that matches the actual behavior of detectors (left) and by
1150choosing a finer mesh (right):
1155<img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_sf.png" alt=
"">
1158<img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_s2.png" alt=
"">
1163One of the important differences between the left and the right figure is that
1164the curves look much less
"angular" at the right. The angularity comes from
1165the fact that while waves in the continuous equation travel equally fast in
1166all directions, this isn
't the case after discretization: there, waves that
1167travel diagonal to cells move at slightly different speeds to those that move
1168parallel to mesh lines. This anisotropy leads to wave fronts that aren't
1169perfectly circular (and would produce sinusoidal signals in the stacked
1170plots), but are bulged out in certain directions. To make things worse, the
1171circular mesh we use (see for example @ref step_6
"step-6" for a view of the
1172coarse mesh) is not isotropic either. The net result is that the signal fronts
1173are not sinusoidal unless the mesh is sufficiently fine. The right image is a
1174lot better in this respect, though artifacts in the form of trailing spurious
1175waves can still be seen.
1178<a name=
"step_24-PlainProg"></a>
1179<h1> The plain program</h1>
1180@include
"step-24.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void add(const size_type i, const size_type j, const number value)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
void hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DataOutBase::CompressionLevel compression_level