<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)
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:
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:
\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
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
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*}
\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]
@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
@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}
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]
@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]
\frac{\partial p(0)}{\partial t} = 0,
@f]
completing the system of equations.
-
<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>
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="">
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.
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.
-
-
-
-
// @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
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
};
- 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
// 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,
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))
{
}
+
template <int dim>
void TATForwardProblem<dim>::solve_v()
{
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);
}
solve_p();
-
system_rhs_v = G2;
laplace_matrix.vmult(tmp, solution_p);
system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
output_results();
-
detector_data << time;
for (unsigned int i = 0; i < detector_locations.size(); ++i)
detector_data << " "
<< " ";
detector_data << std::endl;
-
old_solution_p = solution_p;
old_solution_v = solution_v;
}
projects / dealii.git / commitdiff