532 *
static const std::array<Source, 5>
sources{
555 *
const double radius;
563 * <a name=
"step_24-ImplementationofthecodeTATForwardProblemcodeclass"></a>
568 *
Let's start again with the constructor. Setting the member variables is
569 * straightforward. We use the acoustic wave speed of mineral oil (in
570 * millimeters per microsecond, a common unit in experimental biomedical
571 * imaging) since this is where many of the experiments we want to compare
572 * the output with are made in. The Crank-Nicolson scheme is used again,
573 * i.e. theta is set to 0.5. The time step is later selected to satisfy @f$k =
574 * \frac hc@f$: here we initialize it to an invalid number.
578 * TATForwardProblem<dim>::TATForwardProblem()
580 * , dof_handler(triangulation)
581 * , time_step(std::numeric_limits<double>::quiet_NaN())
583 * , timestep_number(1)
585 * , wave_speed(1.437)
589 * The second task in the constructor is to initialize the array that
590 * holds the detector locations. The results of this program were compared
591 * with experiments in which the step size of the detector spacing is 2.25
592 * degree, corresponding to 160 detector locations. The radius of the
593 * scanning circle is selected to be half way between the center and the
594 * boundary to avoid that the remaining reflections from the imperfect
595 * boundary condition spoils our numerical results.
599 * The locations of the detectors are then calculated in clockwise
600 * order. Note that the following of course only works if we are computing
601 * in 2d, a condition that we guard with an assertion. If we later wanted
602 * to run the same program in 3d, we would have to add code here for the
603 * initialization of detector locations in 3d. Due to the assertion, there
604 * is no way we can forget to do this.
607 * Assert(dim == 2, ExcNotImplemented());
609 * const double detector_step_angle = 2.25;
610 * const double detector_radius = 0.5;
612 * for (double detector_angle = 2 * numbers::PI; detector_angle >= 0;
613 * detector_angle -= detector_step_angle / 360 * 2 * numbers::PI)
614 * detector_locations.push_back(
615 * Point<dim>(std::cos(detector_angle), std::sin(detector_angle)) *
624 * <a name="step_24-TATForwardProblemsetup_system"></a>
625 * <h4>TATForwardProblem::setup_system</h4>
629 * The following system is pretty much what we've already done in @
ref step_23
"step-23",
632 *
new:
we've done so before in @ref step_6 "step-6" and @ref step_10 "step-10", where we also explain
633 * how the PolarManifold or SphericalManifold object places new points on
634 * concentric circles when a cell is refined, which we will use here as
639 * One thing we had to make sure is that the time step satisfies the CFL
640 * condition discussed in the introduction of @ref step_23 "step-23". Back in that program,
641 * we ensured this by hand by setting a timestep that matches the mesh
642 * width, but that was error prone because if we refined the mesh once more
643 * we would also have to make sure the time step is changed. Here, we do
644 * that automatically: we ask a library function for the minimal diameter of
645 * any cell. Then we set @f$k=\frac h{c_0}@f$. The only problem is: what exactly
646 * is @f$h@f$? The point is that there is really no good theory on this question
647 * for the wave equation. It is known that for uniformly refined meshes
648 * consisting of rectangles, @f$h@f$ is the minimal edge length. But for meshes
649 * on general quadrilaterals, the exact relationship appears to be unknown,
650 * i.e. it is unknown what properties of cells are relevant for the CFL
651 * condition. The problem is that the CFL condition follows from knowledge
652 * of the smallest eigenvalue of the Laplace matrix, and that can only be
653 * computed analytically for simply structured meshes.
680 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
683 *
dof_handler.distribute_dofs(fe);
685 *
std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
691 *
sparsity_pattern.copy_from(
dsp);
693 *
system_matrix.reinit(sparsity_pattern);
730 * reason
being that we can't add nonzero entries to a matrix after the
731 * sparsity pattern has been created, so we simply require that the two
732 * matrices have the same sparsity pattern).
765 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
770 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
772 *
for (
const auto &cell : dof_handler.active_cell_iterators())
774 *
if (face->at_boundary())
778 *
fe_values.reinit(cell, face);
781 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
782 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
784 *
fe_values.shape_value(
j,
q_point) *
787 *
cell->get_dof_indices(local_dof_indices);
788 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
789 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
791 *
local_dof_indices[
j],
796 *
system_matrix.copy_from(mass_matrix);
811 *
constraints.close();
818 * <a name=
"step_24-TATForwardProblemsolve_pandTATForwardProblemsolve_v"></a>
837 *
std::cout <<
" p-equation: " << solver_control.last_step()
838 *
<<
" CG iterations." << std::endl;
851 *
std::cout <<
" v-equation: " << solver_control.last_step()
852 *
<<
" CG iterations." << std::endl;
860 * <a name=
"step_24-TATForwardProblemoutput_results"></a>
873 *
data_out.attach_dof_handler(dof_handler);
877 *
data_out.build_patches();
884 *
data_out.write_vtu(output);
892 * <a name=
"step_24-TATForwardProblemrun"></a>
893 * <
h4>TATForwardProblem::run</
h4>
901 *
isn't doing much harm.
905 * The only changes to this function are: first, that we do not have to
906 * project initial values for the velocity @f$v@f$, since we know that it is
907 * zero. And second that we evaluate the solution at the detector locations
908 * computed in the constructor. This is done using the
909 * VectorTools::point_value function. These values are then written to a
910 * file that we open at the beginning of the function.
914 * void TATForwardProblem<dim>::run()
918 * VectorTools::project(dof_handler,
920 * QGauss<dim>(fe.degree + 1),
921 * InitialValuesP<dim>(),
923 * old_solution_v = 0;
926 * std::ofstream detector_data("detectors.dat");
928 * Vector<double> tmp(solution_p.size());
929 * Vector<double> G1(solution_p.size());
930 * Vector<double> G2(solution_v.size());
932 * const double end_time = 0.7;
933 * for (time = time_step; time <= end_time;
934 * time += time_step, ++timestep_number)
936 * std::cout << std::endl;
937 * std::cout << "time_step " << timestep_number << " @ t=" << time
940 * mass_matrix.vmult(G1, old_solution_p);
941 * mass_matrix.vmult(tmp, old_solution_v);
942 * G1.add(time_step * (1 - theta), tmp);
944 * mass_matrix.vmult(G2, old_solution_v);
945 * laplace_matrix.vmult(tmp, old_solution_p);
946 * G2.add(-wave_speed * wave_speed * time_step * (1 - theta), tmp);
948 * boundary_matrix.vmult(tmp, old_solution_p);
949 * G2.add(wave_speed, tmp);
952 * system_rhs_p.add(time_step * theta, G2);
957 * laplace_matrix.vmult(tmp, solution_p);
958 * system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
960 * boundary_matrix.vmult(tmp, solution_p);
961 * system_rhs_v.add(-wave_speed, tmp);
967 * detector_data << time;
968 * for (unsigned int i = 0; i < detector_locations.size(); ++i)
969 * detector_data << ' '
970 * << VectorTools::point_value(dof_handler,
972 * detector_locations[i])
974 * detector_data << std::endl;
976 * old_solution_p = solution_p;
977 * old_solution_v = solution_v;
980 * } // namespace Step24
987 * <a name="step_24-Thecodemaincodefunction"></a>
988 * <h3>The <code>main</code> function</h3>
992 * What remains is the main function of the program. There is nothing here
1000 *
using namespace Step24;
1005 *
catch (std::exception &exc)
1007 *
std::cerr << std::endl
1009 *
<<
"----------------------------------------------------"
1011 *
std::cerr <<
"Exception on processing: " << std::endl
1012 *
<< exc.what() << std::endl
1013 *
<<
"Aborting!" << std::endl
1014 *
<<
"----------------------------------------------------"
1021 *
std::cerr << std::endl
1023 *
<<
"----------------------------------------------------"
1025 *
std::cerr <<
"Unknown exception!" << std::endl
1026 *
<<
"Aborting!" << std::endl
1027 *
<<
"----------------------------------------------------"
1044directions won't contribute to the data. Consequently, we only have to compare
1045our experimental data with two dimensional simulated data.
1047<a name="step_24-Oneabsorber"></a><h3> One absorber </h3>
1050This movie shows the thermoacoustic waves generated by a single small absorber
1051propagating in the medium (in our simulation, we assume the medium is mineral
1052oil, which has a acoustic speed of 1.437 @f$\frac{mm}{\mu s}@f$):
1054<img src="https://www.dealii.org/images/steps/developer/step-24.one_movie.gif" alt="">
1056For a single absorber, we of course have to change the
1057<code>InitialValuesP</code> class accordingly.
1059Next, let us compare experimental and computational results. The visualization
1060uses a technique long used in seismology, where the data of each detector is
1061plotted all in one graph. The way this is done is by offsetting each
1068<
img src=
"https://www.dealii.org/images/steps/developer/step-24.traces.png" alt=
"">
1087<
img src=
"https://www.dealii.org/images/steps/developer/step-24.one.png" alt=
"">
1090<
img src=
"https://www.dealii.org/images/steps/developer/step-24.one_s.png" alt=
"">
1108<
img src=
"https://www.dealii.org/images/steps/developer/step-24.one_sf.png" alt=
"">
1114<
img src=
"https://www.dealii.org/images/steps/developer/step-24.one_s2.png" alt=
"">
1126<
img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_movie.gif" alt=
"">
1133<
img src=
"https://www.dealii.org/images/steps/developer/step-24.multi.png" alt=
"">
1136<
img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_s.png" alt=
"">
1156<
img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_sf.png" alt=
"">
1159<
img src=
"https://www.dealii.org/images/steps/developer/step-24.multi_s2.png" alt=
"">
1168travel diagonal to cells move at slightly different speeds to those that move
1169parallel to mesh lines. This anisotropy leads to wave fronts that aren't
1179<a name=
"step_24-PlainProg"></a>
void add(const size_type i, const size_type j, const number value)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
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.
std::vector< index_type > data
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
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