From b57b7adf0d12337ea2f934f9fa74baf2b164fbd8 Mon Sep 17 00:00:00 2001 From: allmaras Date: Mon, 10 Sep 2007 03:32:23 +0000 Subject: [PATCH] Updated intro for step-29 git-svn-id: https://svn.dealii.org/trunk@15185 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-29/doc/intro.dox | 191 ++++++++++++++++++++++++- 1 file changed, 185 insertions(+), 6 deletions(-) diff --git a/deal.II/examples/step-29/doc/intro.dox b/deal.II/examples/step-29/doc/intro.dox index 1a5015de61..0360423334 100644 --- a/deal.II/examples/step-29/doc/intro.dox +++ b/deal.II/examples/step-29/doc/intro.dox @@ -1,6 +1,185 @@ - -

Introduction

- -This program was contributed by Moritz Allmaras at Texas A&M -University. Some of the work on this tutorial program has been funded -by NSF under grant DMS-0604778. + +

Introduction

+ + +This program was contributed by Moritz Allmaras at Texas A&M +University. Some of the work on this tutorial program has been funded +by NSF under grant DMS-0604778. + + +A question that comes up frequently is how to solve problems involving complex +valued functions with deal.II. For many problems, instead of working with +complex valued finite elements directly, which are not readily available in +the library, it is often much more convenient to split complex valued +functions into their real and imaginary parts and use separate scalar finite +element fields for discretizing each one of them. Basically this amounts to +viewing a single complex valued equation as a system of two real valued +equations. This short example demonstrates how this can be implemented in +deal.II by using an FE_system object to stack two finite element +fields representing real and imaginary parts. We will also revisit the +ParameterHandler class introduced in @ref step_19 "step-19", which provides a +convenient way for reading parameters from a configuration file at runtime +without the need to recompile the program code. + +

Problem setting

+ +The original purpose of this program is to simulate the focussing properties +of an ultrasound wave generated by a transducer lens with variable +geometry. Recent applications in medical imaging use ultrasound waves not only +for imaging porposes, but also to excite certain local effects in a +material, like changes in optical properties, that can then be measured by +other imaging techniques. A vital ingredient for these methods is the ability +to focus the intensity of the ultrasound wave in a particular part of the +material, ideally in a point, to be able to examine the properties of the +material at that particular location. + +To derive a model for this problem, we think of ultrasound as a pressure wave +governed by the wave equation: +@f[ + \frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0 +@f] +where $c$ is the wave speed (that for simplicity we assume to be constant), $U += U(x,t),\;x \in \Omega,\;t\in\mathrm{R}$. The boundary +$\Gamma=\partial\Omega$ is divided into two parts $\Gamma_1$ and +$\Gamma_2=\Gamma\setminus\Gamma_1$, with $\Gamma_1$ representing the +transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to +choose boundary conditions on $\Gamma_2$ in such a way that they imitate a +larger domain). On $\Gamma_1$, the transducer generates a wave of constant +frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): +@f[ +U(x,t) = \sin{\omega t}, \qquad x\in \Gamma_1 +@f] + +If there are no other (interior or boundary) sources, and since the only +source has frequency $\omega$, then the solution admits a separation of +variables of the form $U(x,t) = \textrm{Re}\left(u(x)\,e^{i\omega +t})\right)$. The complex-valued function $u(x)$ describes the spatial +dependency of amplitude and phase (relative to the source) of the waves of +frequency ${\omega}$, with the amplitude being the quantity that we are +interested in. By plugging this form of the solution into the wave equation, +we see that for $u$ we have +@f{eqnarray*} +-\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\ +u(x) &=& 1, \qquad x\in\Gamma_1. +@f} + +For finding suitable conditions on $\Gamma_2$ that model an absorbing +boundary, consider a wave of the form $V(x,t)=e^{i(k\cdot x -\omega t)}$ with +frequency ${\omega}$ traveling in direction $k\in {\mathrm{R}^2}$. In order +for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must +hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right +angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of +$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation +@f[ +c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0. +@f] +Hence, by enforcing the boundary condition +@f[ +c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2, +@f] +waves that hit the boundary $\Gamma_2$ at a right angle will be perfectly +absorbed. On the other hand, those parts of the wave field that do not hit a +boundary at a right angle do not satisfy this condition and enforcing it as a +boundary condition will yield partial reflections, i.e. only parts of the wave +will pass through the boundary as if it wasn't here whereas the remaining +fraction of the wave will be reflected back into the domain. + +If we are willing to accept this as a suffient approximation to an absorbing boundary we finally arrive at the following problem for $u$: +@f{eqnarray*} +-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\ +c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\ +u &=& 1, \qquad x\in\Gamma_1. +@f} +This is a Helmholtz equation (similar to the one in @ref step_7 "step-7", but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system: +@f{eqnarray*} + \left.\begin{array}{ccc} + -\omega^2 v - c^2\Delta v &=& 0 \quad\\ + -\omega^2 w - c^2\Delta w &=& 0 \quad + \end{array}\right\} &\;& x\in\Omega, + \\ + \left.\begin{array}{ccc} + c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\ + c (n\cdot\nabla w) + \omega\,v &=& 0 \quad + \end{array}\right\} &\;& x\in\Gamma_2, + \\ + \left.\begin{array}{ccc} + v &=& 1 \quad\\ + w &=& 0 \quad + \end{array}\right\} &\;& x\in\Gamma_1. +@f} + +For test functions $\phi,\psi$ with $\phi|_{\Gamma_1}=\psi|_{\Gamma_1}=0$, after the usual multiplication, integration over $\Omega$ and applying integration by parts, we get the weak formulation +@f{eqnarray*} +-\omega^2 \langle \phi, v \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \phi, \nabla v \rangle_{\mathrm{L}^2(\Omega)} - c \omega \langle \phi, w \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0, \\ +-\omega^2 \langle \psi, w \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi, \nabla w \rangle_{\mathrm{L}^2(\Omega)} + c \omega \langle \psi, v \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0. +@f} + +We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions +@f[ +v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j. +@f] +Plugging into the variational form yields the equation system +@f[ +\renewcommand{\arraystretch}{2.0} +\left.\begin{array}{ccc} +\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\ +\sum_{j=1}^n \left(-\omega^2 \langle \psi_i, \psi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \psi_i, \nabla \psi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\beta_j + \left(c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\alpha_j &=& 0 +\end{array}\right\}\;\;\forall\; i =1,\ldots,n. +@f] +In matrix notation: +@f[ +\renewcommand{\arraystretch}{2.0} +\left( +\begin{array}{cc} +-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)} & -c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)} \\ +c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)} & -\omega^2 \langle \psi_{i}, \psi_j \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi_{i}, \nabla \psi_j \rangle_{\mathrm{L}^2(\Omega)} +\end{array} +\right) +\left( +\begin{array}{c} +\alpha \\ \beta +\end{array} +\right) += +\left( +\begin{array}{c} +0 \\ 0 +\end{array} +\right) +@f] +(One should not be fooled by the right hand side being zero here, that is because we haven't included the Dirichlet boundary data yet.) +Because of the alternating sign in the off-diagonal blocks, we can already see that this system is non-symmetric, in fact it is even indefinite. +Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be +the same. However, we expect the real and imaginary parts of the solution to +have similar properties and will therefore indeed choose $V_h=W_h$. We choose +the current notation using different symbols merely to be able to distinguish +between test functions for $v$ and $w$, as this distinction plays an important +role in the implementation. + + +

The test case

+ +For the computations, we will consider wave propagation in a rectangular +area. Sound is generated by a transducer that we model using a part of the +boundary shaped like a segment of the circle with center at $(0.5, d)$ and a +radius slightly greater than $d$; this shape leads to a focusing of the sound +wave to the center of the circle. Consequently, we take $\Omega$ to be the +unit square $[0,1]^2$ with the piece $[0.4,0.6]\times\{0\}$ removed and +replaced by the transducer lens (see the Results section below for a graphical +depiction of this domain). We let $d$ vary to change the "focus" of the lens +and see how that affects the spatial distribution of the amplitude of $u$, +i.e. we will consider how well focussed $|u|=\sqrt{v^2+w^2}$ is. + +In the program below, we will implement the complex-valued Helmholtz equations +using the formulation with split real and imaginary parts. We will also +discuss how to generate a domain that looks like a square with a slight bulge +simulating the transducer (in the +UltrasoundProblem::make_grid() function), and how to +generate graphical output that not only contains the solution components $v,w$ +but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in +UltrasoundProblem::output_results()). Finally, we will +explain how to use the ParameterHandler class to easily handle the situation +where we want to prescribe the focal distance $d$, wave speed $c$, frequency +$\omega$, and a number of other parameters, in an input file that is read at +run-time, rather than in the source code where we would have to re-compile the +file every time we want to change parameters. -- 2.39.5