.. !split .. _wave:2D3D: Finite difference methods for 2D and 3D wave equations ====================================================== A natural next step is to consider extensions of the methods for various variants of the one-dimensional wave equation to two-dimensional (2D) and three-dimensional (3D) versions of the wave equation. .. _wave:2D3D:models: Multi-dimensional wave equations -------------------------------- The general wave equation in :math:`d` space dimensions, with constant wave velocity :math:`c`, can be written in the compact form .. _Eq:wave:2D3D:model1: .. math:: \tag{263} \frac{\partial^2 u}{\partial t^2} = c^2\nabla^2 u\hbox{ for }\boldsymbol{x}\in\Omega\subset\mathbb{R}^d,\ t\in (0,T] , where .. math:: \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} , in a 2D problem (:math:`d=2`) and .. math:: \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}, in three space dimensions (:math:`d=3`). Many applications involve variable coefficients, and the general wave equation in :math:`d` dimensions is in this case written as .. _Eq:wave:2D3D:model2: .. math:: \tag{264} \varrho\frac{\partial^2 u}{\partial t^2} = \nabla\cdot (q\nabla u) + f\hbox{ for }\boldsymbol{x}\in\Omega\subset\mathbb{R}^d,\ t\in (0,T], which in, e.g., 2D becomes .. _Eq:_auto114: .. math:: \tag{265} \varrho(x,y) \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x}\left( q(x,y) \frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left( q(x,y) \frac{\partial u}{\partial y}\right) + f(x,y,t) {\thinspace .} To save some writing and space we may use the index notation, where subscript :math:`t`, :math:`x`, or :math:`y` means differentiation with respect to that coordinate. For example, .. math:: \begin{align*} \frac{\partial^2 u}{\partial t^2} &= u_{tt},\\ \frac{\partial}{\partial y}\left( q(x,y) \frac{\partial u}{\partial y}\right) &= (q u_y)_y {\thinspace .} \end{align*} These comments extend straightforwardly to 3D, which means that the 3D versions of the two wave PDEs, with and without variable coefficients, can be stated as .. _Eq:wave:2D3D:model1:v2: .. math:: \tag{266} u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f, .. _Eq:wave:2D3D:model2:v2: .. math:: \tag{267} \varrho u_{tt} = (q u_x)_x + (q u_z)_z + (q u_z)_z + f{\thinspace .} At *each point* of the boundary :math:`\partial\Omega` (of :math:`\Omega`) we need *one* boundary condition involving the unknown :math:`u`. The boundary conditions are of three principal types: 1. :math:`u` is prescribed (:math:`u=0` or a known time variation of :math:`u` at the boundary points, e.g., modeling an incoming wave), 2. :math:`\partial u/\partial n = \boldsymbol{n}\cdot\nabla u` is prescribed (zero for reflecting boundaries), 3. an open boundary condition (also called radiation condition) is specified to let waves travel undisturbed out of the domain, see :ref:`wave:app:exer:radiationBC` for details. All the listed wave equations with *second-order* derivatives in time need *two* initial conditions: 1. :math:`u=I`, 2. :math:`u_t = V`. .. _wave:2D3D:mesh: Mesh ---- We introduce a mesh in time and in space. The mesh in time consists of time points .. math:: t_0=0 < t_1 < \cdots < t_{N_t}, normally, for wave equation problems, with a constant spacing :math:`\Delta t= t_{n+1}-t_{n}`, :math:`n\in{{\mathcal{I^-}_t}}`. Finite difference methods are easy to implement on simple rectangle- or box-shaped *spatial domains*. More complicated shapes of the spatial domain require substantially more advanced techniques and implementational efforts (and a finite element method is usually a more convenient approach). On a rectangle- or box-shaped domain, mesh points are introduced separately in the various space directions: .. math:: \begin{align*} &x_0 < x_1 < \cdots < x_{N_x} \hbox{ in the }x \hbox{ direction},\\ &y_0 < y_1 < \cdots < y_{N_y} \hbox{ in the }y \hbox{ direction},\\ &z_0 < z_1 < \cdots < z_{N_z} \hbox{ in the }z \hbox{ direction}{\thinspace .} \end{align*} We can write a general mesh point as :math:`(x_i,y_j,z_k,t_n)`, with :math:`i\in{\mathcal{I}_x}`, :math:`j\in{\mathcal{I}_y}`, :math:`k\in{\mathcal{I}_z}`, and :math:`n\in{\mathcal{I}_t}`. It is a very common choice to use constant mesh spacings: :math:`\Delta x = x_{i+1}-x_{i}`, :math:`i\in{{\mathcal{I^-}_x}}`, :math:`\Delta y = y_{j+1}-y_{j}`, :math:`j\in{{\mathcal{I^-}_y}}`, and :math:`\Delta z = z_{k+1}-z_{k}`, :math:`k\in{{\mathcal{I^-}_z}}`. With equal mesh spacings one often introduces :math:`h = \Delta x = \Delta y =\Delta z`. The unknown :math:`u` at mesh point :math:`(x_i,y_j,z_k,t_n)` is denoted by :math:`u^{n}_{i,j,k}`. In 2D problems we just skip the :math:`z` coordinate (by assuming no variation in that direction: :math:`\partial/\partial z=0`) and write :math:`u^n_{i,j}`. .. _wave:2D3D:discretization: Discretization (1) --------------------------- Two- and three-dimensional wave equations are easily discretized by assembling building blocks for discretization of 1D wave equations, because the multi-dimensional versions just contain terms of the same type as those in 1D. Discretizing the PDEs ~~~~~~~~~~~~~~~~~~~~~ Equation :ref:`(266) ` can be discretized as .. _Eq:_auto115: .. math:: \tag{268} [D_tD_t u = c^2(D_xD_x u + D_yD_yu + D_zD_z u) + f]^n_{i,j,k} {\thinspace .} A 2D version might be instructive to write out in detail: .. math:: [D_tD_t u = c^2(D_xD_x u + D_yD_yu) + f]^n_{i,j,k}, which becomes .. math:: \frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2} = c^2 \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2} + c^2 \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2} + f^n_{i,j}, Assuming, as usual, that all values at time levels :math:`n` and :math:`n-1` are known, we can solve for the only unknown :math:`u^{n+1}_{i,j}`. The result can be compactly written as .. _Eq:wave:2D3D:models:unp1: .. math:: \tag{269} u^{n+1}_{i,j} = 2u^n_{i,j} + u^{n-1}_{i,j} + c^2\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}{\thinspace .} As in the 1D case, we need to develop a special formula for :math:`u^1_{i,j}` where we combine the general scheme for :math:`u^{n+1}_{i,j}`, when :math:`n=0`, with the discretization of the initial condition: .. math:: [D_{2t}u = V]^0_{i,j}\quad\Rightarrow\quad u^{-1}_{i,j} = u^1_{i,j} - 2\Delta t V_{i,j} {\thinspace .} The result becomes, in compact form, .. _Eq:wave:2D3D:models:u1: .. math:: \tag{270} u^{1}_{i,j} = u^0_{i,j} -2\Delta V_{i,j} + {\frac{1}{2}} c^2\Delta t^2[D_xD_x u + D_yD_y u]^0_{i,j}{\thinspace .} The PDE :ref:`(267) ` with variable coefficients is discretized term by term using the corresponding elements from the 1D case: .. _Eq:_auto116: .. math:: \tag{271} [\varrho D_tD_t u = (D_x\overline{q}^x D_x u + D_y\overline{q}^y D_yu + D_z\overline{q}^z D_z u) + f]^n_{i,j,k} {\thinspace .} When written out and solved for the unknown :math:`u^{n+1}_{i,j,k}`, one gets the scheme .. math:: \begin{align*} u^{n+1}_{i,j,k} &= - u^{n-1}_{i,j,k} + 2u^{n}_{i,j,k} + \\ &\quad \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i+1,j,k})(u^{n}_{i+1,j,k} - u^{n}_{i,j,k}) - \\ &\qquad\qquad \frac{1}{2}(q_{i-1,j,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i-1,j,k})) + \\ &\quad \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i,j+1,k})(u^{n}_{i,j+1,k} - u^{n}_{i,j,k}) - \\ &\qquad\qquad\frac{1}{2}(q_{i,j-1,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j-1,k})) + \\ &\quad \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \frac{1}{2}(q_{i,j,k} + q_{i,j,k+1})(u^{n}_{i,j,k+1} - u^{n}_{i,j,k}) -\\ &\qquad\qquad \frac{1}{2}(q_{i,j,k-1} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j,k-1})) + \\ &\quad \Delta t^2 f^n_{i,j,k} {\thinspace .} \end{align*} Also here we need to develop a special formula for :math:`u^1_{i,j,k}` by combining the scheme for :math:`n=0` with the discrete initial condition, which is just a matter of inserting :math:`u^{-1}_{i,j,k}=u^1_{i,j,k} - 2\Delta tV_{i,j,k}` in the scheme and solving for :math:`u^1_{i,j,k}`. Handling boundary conditions where :math:`u` is known ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The schemes listed above are valid for the internal points in the mesh. After updating these, we need to visit all the mesh points at the boundaries and set the prescribed :math:`u` value. Discretizing the Neumann condition ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The condition :math:`\partial u/\partial n = 0` was implemented in 1D by discretizing it with a :math:`D_{2x}u` centered difference, followed by eliminating the fictitious :math:`u` point outside the mesh by using the general scheme at the boundary point. Alternatively, one can introduce ghost cells and update a ghost value for use in the Neumann condition. Exactly the same ideas are reused in multiple dimensions. Consider the condition :math:`\partial u/\partial n = 0` at a boundary :math:`y=0` of a rectangular domain :math:`[0, L_x]\times [0,L_y]` in 2D. The normal direction is then in :math:`-y` direction, so .. math:: \frac{\partial u}{\partial n} = -\frac{\partial u}{\partial y}, and we set .. math:: [-D_{2y} u = 0]^n_{i,0}\quad\Rightarrow\quad \frac{u^n_{i,1}-u^n_{i,-1}}{2\Delta y} = 0 {\thinspace .} From this it follows that :math:`u^n_{i,-1}=u^n_{i,1}`. The discretized PDE at the boundary point :math:`(i,0)` reads .. math:: \frac{u^{n+1}_{i,0} - 2u^{n}_{i,0} + u^{n-1}_{i,0}}{\Delta t^2} = c^2 \frac{u^{n}_{i+1,0} - 2u^{n}_{i,0} + u^{n}_{i-1,0}}{\Delta x^2} + c^2 \frac{u^{n}_{i,1} - 2u^{n}_{i,0} + u^{n}_{i,-1}}{\Delta y^2} + f^n_{i,j}, We can then just insert :math:`u^n_{i,1}` for :math:`u^n_{i,-1}` in this equation and solve for the boundary value :math:`u^{n+1}_{i,0}`, just as was done in 1D. From these calculations, we see a pattern: the general scheme applies at the boundary :math:`j=0` too if we just replace :math:`j-1` by :math:`j+1`. Such a pattern is particularly useful for implementations. The details follow from the explained 1D case in the section :ref:`wave:pde2:Neumann:impl`. The alternative approach to eliminating fictitious values outside the mesh is to have :math:`u^n_{i,-1}` available as a ghost value. The mesh is extended with one extra line (2D) or plane (3D) of ghost cells at a Neumann boundary. In the present example it means that we need a line with ghost cells below the :math:`y` axis. The ghost values must be updated according to :math:`u^{n+1}_{i,-1}=u^{n+1}_{i,1}`. .. _wave:2D3D:impl: Implementation (6) =========================== .. index:: single: wave equation; 2D, implementation We shall now describe in detail various Python implementations for solving a standard 2D, linear wave equation with constant wave velocity and :math:`u=0` on the boundary. The wave equation is to be solved in the space-time domain :math:`\Omega\times (0,T]`, where :math:`\Omega = (0,L_x)\times (0,L_y)` is a rectangular spatial domain. More precisely, the complete initial-boundary value problem is defined by .. _Eq:_auto117: .. math:: \tag{272} u_{tt} = c^2(u_{xx} + u_{yy}) + f(x,y,t),\quad (x,y)\in \Omega,\ t\in (0,T], .. _Eq:_auto118: .. math:: \tag{273} u(x,y,0) = I(x,y),\quad (x,y)\in\Omega, .. _Eq:_auto119: .. math:: \tag{274} u_t(x,y,0) = V(x,y),\quad (x,y)\in\Omega, .. _Eq:_auto120: .. math:: \tag{275} u = 0,\quad (x,y)\in\partial\Omega,\ t\in (0,T], where :math:`\partial\Omega` is the boundary of :math:`\Omega`, in this case the four sides of the rectangle :math:`\Omega = [0,L_x]\times [0,L_y]`: :math:`x=0`, :math:`x=L_x`, :math:`y=0`, and :math:`y=L_y`. The PDE is discretized as .. math:: [D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j}, which leads to an explicit updating formula to be implemented in a program: .. math:: u^{n+1}_{i,j} = -u^{n-1}_{i,j} + 2u^n_{i,j} + \nonumber .. _Eq:wave:2D3D:impl1:2Du0:ueq:discrete: .. math:: \tag{276} \quad C_x^2( u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}) + C_y^2 (u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}) + \Delta t^2 f_{i,j}^n, for all interior mesh points :math:`i\in{{\mathcal{I^i}_x}}` and :math:`j\in{{\mathcal{I^i}_y}}`, for :math:`n\in{{\mathcal{I^+}_t}}`. The constants :math:`C_x` and :math:`C_y` are defined as .. math:: C_x = c\frac{\Delta t}{\Delta x},\quad C_x = c\frac{\Delta t}{\Delta y} {\thinspace .} At the boundary, we simply set :math:`u^{n+1}_{i,j}=0` for :math:`i=0`, :math:`j=0,\ldots,N_y`; :math:`i=N_x`, :math:`j=0,\ldots,N_y`; :math:`j=0`, :math:`i=0,\ldots,N_x`; and :math:`j=N_y`, :math:`i=0,\ldots,N_x`. For the first step, :math:`n=0`, :ref:`(276) ` is combined with the discretization of the initial condition :math:`u_t=V`, :math:`[D_{2t} u = V]^0_{i,j}` to obtain a special formula for :math:`u^1_{i,j}` at the interior mesh points: .. math:: u^{1}_{i,j} = u^0_{i,j} + \Delta t V_{i,j} + \nonumber .. math:: \quad {\frac{1}{2}}C_x^2( u^{0}_{i+1,j} - 2u^{0}_{i,j} + u^{0}_{i-1,j}) + {\frac{1}{2}}C_y^2 (u^{0}_{i,j+1} - 2u^{0}_{i,j} + u^{0}_{i,j-1}) +\nonumber .. _Eq:wave:2D3D:impl1:2Du0:ueq:n0:discrete: .. math:: \tag{277} \quad \frac{1}{2}\Delta t^2f_{i,j}^n, The algorithm is very similar to the one in 1D: 1. Set initial condition :math:`u^0_{i,j}=I(x_i,y_j)` 2. Compute :math:`u^1_{i,j}` from :ref:`(276) ` 3. Set :math:`u^1_{i,j}=0` for the boundaries :math:`i=0,N_x`, :math:`j=0,N_y` 4. For :math:`n=1,2,\ldots,N_t`: a. Find :math:`u^{n+1}_{i,j}` from :ref:`(276) ` for all internal mesh points, :math:`i\in{{\mathcal{I^i}_x}}`, :math:`j\in{{\mathcal{I^i}_y}}` b. Set :math:`u^{n+1}_{i,j}=0` for the boundaries :math:`i=0,N_x`, :math:`j=0,N_y` .. _wave2D3D:impl:scalar: Scalar computations ------------------- The ``solver`` function for a 2D case with constant wave velocity and boundary condition :math:`u=0` is analogous to the 1D case with similar parameter values (see ``wave1D_u0.py``), apart from a few necessary extensions. The code is found in the program `wave2D_u0.py `__. Domain and mesh ~~~~~~~~~~~~~~~ The spatial domain is now :math:`[0,L_x]\times [0,L_y]`, specified by the arguments ``Lx`` and ``Ly``. Similarly, the number of mesh points in the :math:`x` and :math:`y` directions, :math:`N_x` and :math:`N_y`, become the arguments ``Nx`` and ``Ny``. In multi-dimensional problems it makes less sense to specify a Courant number since the wave velocity is a vector and mesh spacings may differ in the various spatial directions. We therefore give :math:`\Delta t` explicitly. The signature of the ``solver`` function is then .. code-block:: python def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T, user_action=None, version='scalar'): Key parameters used in the calculations are created as .. code-block:: python x = linspace(0, Lx, Nx+1) # mesh points in x dir y = linspace(0, Ly, Ny+1) # mesh points in y dir dx = x[1] - x[0] dy = y[1] - y[0] Nt = int(round(T/float(dt))) t = linspace(0, N*dt, N+1) # mesh points in time Cx2 = (c*dt/dx)**2; Cy2 = (c*dt/dy)**2 # help variables dt2 = dt**2 Solution arrays ~~~~~~~~~~~~~~~ We store :math:`u^{n+1}_{i,j}`, :math:`u^{n}_{i,j}`, and :math:`u^{n-1}_{i,j}` in three two-dimensional arrays, .. code-block:: python u = zeros((Nx+1,Ny+1)) # solution array u_n = [zeros((Nx+1,Ny+1)), zeros((Nx+1,Ny+1))] # t-dt, t-2*dt where :math:`u^{n+1}_{i,j}` corresponds to ``u[i,j]``, :math:`u^{n}_{i,j}` to ``u_n[i,j]``, and :math:`u^{n-1}_{i,j}` to ``u_nm1[i,j]``. .. index:: index set notation Index sets ~~~~~~~~~~ It is also convenient to introduce the index sets (cf. the section :ref:`wave:indexset`) .. code-block:: python Ix = range(0, u.shape[0]) Iy = range(0, u.shape[1]) It = range(0, t.shape[0]) Computing the solution ~~~~~~~~~~~~~~~~~~~~~~ Inserting the initial condition ``I`` in ``u_n`` and making a callback to the user in terms of the ``user_action`` function is a straightforward generalization of the 1D code from the section :ref:`wave:string:impl`: .. code-block:: python for i in Ix: for j in Iy: u_n[i,j] = I(x[i], y[j]) if user_action is not None: user_action(u_n, x, xv, y, yv, t, 0) The ``user_action`` function has additional arguments compared to the 1D case. The arguments ``xv`` and ``yv`` will be commented upon in the section :ref:`wave2D3D:impl:vectorized`. The key finite difference formula :ref:`(269) ` for updating the solution at a time level is implemented in a separate function as .. code-block:: python def advance_scalar(u, u_n, u_nm1, f, x, y, t, n, Cx2, Cy2, dt2, V=None, step1=False): Ix = range(0, u.shape[0]); Iy = range(0, u.shape[1]) if step1: dt = sqrt(dt2) # save Cx2 = 0.5*Cx2; Cy2 = 0.5*Cy2; dt2 = 0.5*dt2 # redefine D1 = 1; D2 = 0 else: D1 = 2; D2 = 1 for i in Ix[1:-1]: for j in Iy[1:-1]: u_xx = u_n[i-1,j] - 2*u_n[i,j] + u_n[i+1,j] u_yy = u_n[i,j-1] - 2*u_n[i,j] + u_n[i,j+1] u[i,j] = D1*u_n[i,j] - D2*u_nm1[i,j] + \ Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n]) if step1: u[i,j] += dt*V(x[i], y[j]) # Boundary condition u=0 j = Iy[0] for i in Ix: u[i,j] = 0 j = Iy[-1] for i in Ix: u[i,j] = 0 i = Ix[0] for j in Iy: u[i,j] = 0 i = Ix[-1] for j in Iy: u[i,j] = 0 return u The ``step1`` variable has been introduced to allow the formula to be reused for the first step, computing :math:`u^1_{i,j}`: .. code-block:: python u = advance_scalar(u, u_n, f, x, y, t, n, Cx2, Cy2, dt, V, step1=True) Below, we will make many alternative implementations of the ``advance_scalar`` function to speed up the code since most of the CPU time in simulations is spent in this function. .. admonition:: Remark: How to use the solution The ``solver`` function in the ``wave2D_u0.py`` code updates arrays for the next time step by switching references as described in the section :ref:`wave:pde1:impl:ref:switch`. Any use of ``u`` on the user's side is assumed to take place in the user action function. However, should the code be changed such that ``u`` is returned and used as solution, have in mind that you must return ``u_n`` after the time lime, otherwise a ``return u`` will actually return ``u_nm1`` (due to the switching of array indices in the loop)! .. _wave2D3D:impl:vectorized: Vectorized computations ----------------------- The scalar code above turns out to be extremely slow for large 2D meshes, and probably useless in 3D beyond debugging of small test cases. Vectorization is therefore a must for multi-dimensional finite difference computations in Python. For example, with a mesh consisting of :math:`30\times 30` cells, vectorization brings down the CPU time by a factor of 70 (!). Equally important, vectorized code can also easily be parallelized to take (usually) optimal advantage of parallel computer platforms. In the vectorized case, we must be able to evaluate user-given functions like :math:`I(x,y)` and :math:`f(x,y,t)` for the entire mesh in one operation (without loops). These user-given functions are provided as Python functions ``I(x,y)`` and ``f(x,y,t)``, respectively. Having the one-dimensional coordinate arrays ``x`` and ``y`` is not sufficient when calling ``I`` and ``f`` in a vectorized way. We must extend ``x`` and ``y`` to their vectorized versions ``xv`` and ``yv``: .. code-block:: python from numpy import newaxis xv = x[:,newaxis] yv = y[newaxis,:] # or xv = x.reshape((x.size, 1)) yv = y.reshape((1, y.size)) This is a standard required technique when evaluating functions over a 2D mesh, say ``sin(xv)*cos(xv)``, which then gives a result with shape ``(Nx+1,Ny+1)``. Calling ``I(xv, yv)`` and ``f(xv, yv, t[n])`` will now return ``I`` and ``f`` values for the entire set of mesh points. With the ``xv`` and ``yv`` arrays for vectorized computing, setting the initial condition is just a matter of .. code-block:: python u_n[:,:] = I(xv, yv) One could also have written ``u_n = I(xv, yv)`` and let ``u_n`` point to a new object, but vectorized operations often make use of direct insertion in the original array through ``u_n[:,:]``, because sometimes not all of the array is to be filled by such a function evaluation. This is the case with the computational scheme for :math:`u^{n+1}_{i,j}`: .. code-block:: python def advance_vectorized(u, u_n, u_nm1, f_a, Cx2, Cy2, dt2, V=None, step1=False): if step1: dt = np.sqrt(dt2) # save Cx2 = 0.5*Cx2; Cy2 = 0.5*Cy2; dt2 = 0.5*dt2 # redefine D1 = 1; D2 = 0 else: D1 = 2; D2 = 1 u_xx = u_n[:-2,1:-1] - 2*u_n[1:-1,1:-1] + u_n[2:,1:-1] u_yy = u_n[1:-1,:-2] - 2*u_n[1:-1,1:-1] + u_n[1:-1,2:] u[1:-1,1:-1] = D1*u_n[1:-1,1:-1] - D2*u_nm1[1:-1,1:-1] + \ Cx2*u_xx + Cy2*u_yy + dt2*f_a[1:-1,1:-1] if step1: u[1:-1,1:-1] += dt*V[1:-1, 1:-1] # Boundary condition u=0 j = 0 u[:,j] = 0 j = u.shape[1]-1 u[:,j] = 0 i = 0 u[i,:] = 0 i = u.shape[0]-1 u[i,:] = 0 return u Array slices in 2D are more complicated to understand than those in 1D, but the logic from 1D applies to each dimension separately. For example, when doing :math:`u^{n}_{i,j} - u^{n}_{i-1,j}` for :math:`i\in{{\mathcal{I^+}_x}}`, we just keep ``j`` constant and make a slice in the first index: ``u_n[1:,j] - u_n[:-1,j]``, exactly as in 1D. The ``1:`` slice specifies all the indices :math:`i=1,2,\ldots,N_x` (up to the last valid index), while ``:-1`` specifies the relevant indices for the second term: :math:`0,1,\ldots,N_x-1` (up to, but not including the last index). In the above code segment, the situation is slightly more complicated, because each displaced slice in one direction is accompanied by a ``1:-1`` slice in the other direction. The reason is that we only work with the internal points for the index that is kept constant in a difference. The boundary conditions along the four sides makes use of a slice consisting of all indices along a boundary: .. code-block:: python u[: ,0] = 0 u[:,Ny] = 0 u[0 ,:] = 0 u[Nx,:] = 0 In the vectorized update of ``u`` (above), the function ``f`` is first computed as an array over all mesh points: .. code-block:: text f_a = f(xv, yv, t[n]) We could, alternatively, have used the call ``f(xv, yv, t[n])[1:-1,1:-1]`` in the last term of the update statement, but other implementations in compiled languages benefit from having ``f`` available in an array rather than calling our Python function ``f(x,y,t)`` for every point. Also in the ``advance_vectorized`` function we have introduced a boolean ``step1`` to reuse the formula for the first time step in the same way as we did with ``advance_scalar``. We refer to the ``solver`` function in ``wave2D_u0.py`` for the details on how the overall algorithm is implemented. The callback function now has the arguments ``u, x, xv, y, yv, t, n``. The inclusion of ``xv`` and ``yv`` makes it easy to, e.g., compute an exact 2D solution in the callback function and compute errors, through an expression like ``u - u_exact(xv, yv, t[n])``. .. _wave2D3D:impl:verify: Verification (6) ------------------------- Testing a quadratic solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The 1D solution from the section :ref:`wave:pde2:fd:verify:quadratic` can be generalized to multi-dimensions and provides a test case where the exact solution also fulfills the discrete equations, such that we know (to machine precision) what numbers the solver function should produce. In 2D we use the following generalization of :ref:`(189) `: .. _Eq:wave2D3D:impl:verify:quadratic: .. math:: \tag{278} {u_{\small\mbox{e}}}(x,y,t) = x(L_x-x)y(L_y-y)(1+{\frac{1}{2}}t) {\thinspace .} This solution fulfills the PDE problem if :math:`I(x,y)={u_{\small\mbox{e}}}(x,y,0)`, :math:`V=\frac{1}{2}{u_{\small\mbox{e}}}(x,y,0)`, and :math:`f=2c^2(1+{\frac{1}{2}}t)(y(L_y-y) + x(L_x-x))`. To show that :math:`{u_{\small\mbox{e}}}` also solves the discrete equations, we start with the general results :math:`[D_t D_t 1]^n=0`, :math:`[D_t D_t t]^n=0`, and :math:`[D_t D_t t^2]=2`, and use these to compute .. math:: \begin{align*} [D_xD_x {u_{\small\mbox{e}}}]^n_{i,j} &= [y(L_y-y)(1+{\frac{1}{2}}t) D_xD_x x(L_x-x)]^n_{i,j}\\ &= y_j(L_y-y_j)(1+{\frac{1}{2}}t_n)(-2){\thinspace .} \end{align*} A similar calculation must be carried out for the :math:`[D_yD_y {u_{\small\mbox{e}}}]^n_{i,j}` and :math:`[D_tD_t {u_{\small\mbox{e}}}]^n_{i,j}` terms. One must also show that the quadratic solution fits the special formula for :math:`u^1_{i,j}`. The details are left as :ref:`wave:exer:quadratic:2D`. The ``test_quadratic`` function in the `wave2D_u0.py `__ program implements this verification as a proper test function for the pytest and nose frameworks. Visualization (2) -------------------------- .. index:: 3D visualization .. index:: visualization of 2D scalar fields Eventually, we are ready with a real application with our code! Look at the ``wave2D_u0.py`` and the ``gaussian`` function. It starts with a Gaussian function to see how it propagates in a square with :math:`u=0` on the boundaries: .. code-block:: python def gaussian(plot_method=2, version='vectorized', save_plot=True): """ Initial Gaussian bell in the middle of the domain. plot_method=1 applies mesh function, =2 means surf, =3 means Matplotlib, =4 means mayavi, =0 means no plot. """ # Clean up plot files for name in glob('tmp_*.png'): os.remove(name) Lx = 10 Ly = 10 c = 1.0 from numpy import exp def I(x, y): """Gaussian peak at (Lx/2, Ly/2).""" return exp(-0.5*(x-Lx/2.0)**2 - 0.5*(y-Ly/2.0)**2) def plot_u(u, x, xv, y, yv, t, n): """User action function for plotting.""" ... Nx = 40; Ny = 40; T = 20 dt, cpu = solver(I, None, None, c, Lx, Ly, Nx, Ny, -1, T, user_action=plot_u, version=version) Matplotlib ~~~~~~~~~~ We want to animate a 3D surface in Matplotlib, but this is a really slow process and not recommended, so we consider Matplotlib not an option as long as on-screen animation is desired. One can use the recipes for single shots of :math:`u`, where it does produce high-quality 3D plots. .. index:: Gnuplot Gnuplot ~~~~~~~ Let us look at different ways for visualization. We import SciTools as ``st`` and can access ``st.mesh`` and ``st.surf`` in Matplotlib or Gnuplot, but this is not supported except for the Gnuplot package, where it works really well (Figure :ref:`wave2D3D:impl:viz:fig:gnuplot1`). Then we choose ``plot_method=2`` (or less relevant ``plot_method=1``) and force the backend for SciTools to be Gnuplot (if you have the C package Gnuplot and the ``Gnuplot.py`` Python interface module installed): .. code-block:: text Terminal> python wave2D_u0.py --SCITOOLS_easyviz_backend gnuplot It gives a nice visualization with lifted surface and contours beneath. Figure :ref:`wave2D3D:impl:viz:fig:gnuplot1` shows four plots of :math:`u`. .. _wave2D3D:impl:viz:fig:gnuplot1: .. figure:: wave2D_u0_gnuplot_gaussian.png :width: 800 *Snapshots of the surface plotted by Gnuplot* Video files can be made of the PNG frames: .. code-block:: text Terminal> ffmpeg -i tmp_%04d.png -r 25 -vcodec flv movie.flv Terminal> ffmpeg -i tmp_%04d.png -r 25 -vcodec linx264 movie.mp4 Terminal> ffmpeg -i tmp_%04d.png -r 25 -vcodec libvpx movie.webm Terminal> ffmpeg -i tmp_%04d.png -r 25 -vcodec libtheora movie.ogg It is wise to use a high frame rate - a low one will just skip many frames. There may also be considerable quality differences between the different formats. .. raw:: html

.. index:: Mayavi Mayavi ~~~~~~ The best option for doing visualization of 2D and 3D scalar and vector fields in Python programs is Mayavi, which is an interface to the high-quality package VTK in C++. There is good online documentation and also an introduction in Chapter 5 of [Ref06]_. To obtain Mayavi on Ubuntu platforms you can write .. code-block:: text pip install mayavi --upgrade For Mac OS X and Windows, we recommend using Anaconda. To obtain Mayavi for Anaconda you can write .. code-block:: text conda install mayavi Mayavi has a MATLAB-like interface called ``mlab``. We can do .. code-block:: python import mayavi.mlab as plt # or from mayavi import mlab and have ``plt`` (as usual) or ``mlab`` as a kind of MATLAB visualization access inside our program (just more powerful and with higher visual quality). The official documentation of the ``mlab`` module is provided in two places, one for the `basic functionality `__ and one for `further functionality `__. Basic `figure handling `__ is very similar to the one we know from Matplotlib. Just as for Matplotlib, all plotting commands you do in ``mlab`` will go into the same figure, until you manually change to a new figure. Back to our application, the following code for the user action function with plotting in Mayavi is relevant to add. .. code-block:: python # Top of the file try: import mayavi.mlab as mlab except: # We don't have mayavi pass def solver(...): ... def gaussian(...): ... if plot_method == 3: from mpl_toolkits.mplot3d import axes3d import matplotlib.pyplot as plt from matplotlib import cm plt.ion() fig = plt.figure() u_surf = None def plot_u(u, x, xv, y, yv, t, n): """User action function for plotting.""" if t[n] == 0: time.sleep(2) if plot_method == 1: # Works well with Gnuplot backend, not with Matplotlib st.mesh(x, y, u, title='t=%g' % t[n], zlim=[-1,1], caxis=[-1,1]) elif plot_method == 2: # Works well with Gnuplot backend, not with Matplotlib st.surfc(xv, yv, u, title='t=%g' % t[n], zlim=[-1, 1], colorbar=True, colormap=st.hot(), caxis=[-1,1], shading='flat') elif plot_method == 3: print 'Experimental 3D matplotlib...not recommended' elif plot_method == 4: # Mayavi visualization mlab.clf() extent1 = (0, 20, 0, 20,-2, 2) s = mlab.surf(x , y, u, colormap='Blues', warp_scale=5,extent=extent1) mlab.axes(s, color=(.7, .7, .7), extent=extent1, ranges=(0, 10, 0, 10, -1, 1), xlabel='', ylabel='', zlabel='', x_axis_visibility=False, z_axis_visibility=False) mlab.outline(s, color=(0.7, .7, .7), extent=extent1) mlab.text(6, -2.5, '', z=-4, width=0.14) mlab.colorbar(object=None, title=None, orientation='horizontal', nb_labels=None, nb_colors=None, label_fmt=None) mlab.title('Gaussian t=%g' % t[n]) mlab.view(142, -72, 50) f = mlab.gcf() camera = f.scene.camera camera.yaw(0) if plot_method > 0: time.sleep(0) # pause between frames if save_plot: filename = 'tmp_%04d.png' % n if plot_method == 4: mlab.savefig(filename) # time consuming! elif plot_method in (1,2): st.savefig(filename) # time consuming! This is a point to get started - visualization is as always a very time-consuming and experimental discipline. With the PNG files we can use ``ffmpeg`` to create videos. .. figure:: mayavi2D_gaussian1.png :width: 600 Plot with Mayavi .. raw:: html

Exercises (4) ====================== .. --- begin exercise --- .. _wave:exer:quadratic:2D: Exercise 2.16: Check that a solution fulfills the discrete model ---------------------------------------------------------------- Carry out all mathematical details to show that :ref:`(278) ` is indeed a solution of the discrete model for a 2D wave equation with :math:`u=0` on the boundary. One must check the boundary conditions, the initial conditions, the general discrete equation at a time level and the special version of this equation for the first time level. Filename: ``check_quadratic_solution``. .. --- end exercise --- .. --- begin exercise --- .. _wave:exer:mesh3D:calculus: Project 2.17: Calculus with 2D mesh functions --------------------------------------------- The goal of this project is to redo :ref:`wave:exer:mesh1D:calculus` with 2D mesh functions (:math:`f_{i,j}`). **Differentiation.** The differentiation results in a discrete gradient function, which in the 2D case can be represented by a three-dimensional array ``df[d,i,j]`` where ``d`` represents the direction of the derivative, and ``i,j`` is a mesh point in 2D. Use centered differences for the derivative at inner points and one-sided forward or backward differences at the boundary points. Construct unit tests and write a corresponding test function. **Integration.** The integral of a 2D mesh function :math:`f_{i,j}` is defined as .. math:: F_{i,j} = \int_{y_0}^{y_j} \int_{x_0}^{x_i} f(x,y)dxdy, where :math:`f(x,y)` is a function that takes on the values of the discrete mesh function :math:`f_{i,j}` at the mesh points, but can also be evaluated in between the mesh points. The particular variation between mesh points can be taken as bilinear, but this is not important as we will use a product Trapezoidal rule to approximate the integral over a cell in the mesh and then we only need to evaluate :math:`f(x,y)` at the mesh points. Suppose :math:`F_{i,j}` is computed. The calculation of :math:`F_{i+1,j}` is then .. math:: \begin{align*} F_{i+1,j} &= F_{i,j} + \int_{x_i}^{x_{i+1}}\int_{y_0}^{y_j} f(x,y)dydx\\ & \approx \Delta x \frac{1}{2}\left( \int_{y_0}^{y_j} f(x_{i},y)dy + \int_{y_0}^{y_j} f(x_{i+1},y)dy\right) \end{align*} The integrals in the :math:`y` direction can be approximated by a Trapezoidal rule. A similar idea can be used to compute :math:`F_{i,j+1}`. Thereafter, :math:`F_{i+1,j+1}` can be computed by adding the integral over the final corner cell to :math:`F_{i+1,j} + F_{i,j+1} - F_{i,j}`. Carry out the details of these computations and implement a function that can return :math:`F_{i,j}` for all mesh indices :math:`i` and :math:`j`. Use the fact that the Trapezoidal rule is exact for linear functions and write a test function. Filename: ``mesh_calculus_2D``. .. --- end exercise --- .. --- begin exercise --- .. _wave:app:exer:wave2D:Neumann: Exercise 2.18: Implement Neumann conditions in 2D ------------------------------------------------- Modify the `wave2D_u0.py `__ program, which solves the 2D wave equation :math:`u_{tt}=c^2(u_{xx}+u_{yy})` with constant wave velocity :math:`c` and :math:`u=0` on the boundary, to have Neumann boundary conditions: :math:`\partial u/\partial n=0`. Include both scalar code (for debugging and reference) and vectorized code (for speed). To test the code, use :math:`u=1.2` as solution (:math:`I(x,y)=1.2`, :math:`V=f=0`, and :math:`c` arbitrary), which should be exactly reproduced with any mesh as long as the stability criterion is satisfied. Another test is to use the plug-shaped pulse in the ``pulse`` function from the section :ref:`wave:pde2:software` and the `wave1D_dn_vc.py `__ program. This pulse is exactly propagated in 1D if :math:`c\Delta t/\Delta x=1`. Check that also the 2D program can propagate this pulse exactly in :math:`x` direction (:math:`c\Delta t/\Delta x=1`, :math:`\Delta y` arbitrary) and :math:`y` direction (:math:`c\Delta t/\Delta y=1`, :math:`\Delta x` arbitrary). Filename: ``wave2D_dn``. .. --- end exercise --- .. --- begin exercise --- .. _wave:exer:3D:f77:cy:efficiency: Exercise 2.19: Test the efficiency of compiled loops in 3D ---------------------------------------------------------- Extend the ``wave2D_u0.py`` code and the Cython, Fortran, and C versions to 3D. Set up an efficiency experiment to determine the relative efficiency of pure scalar Python code, vectorized code, Cython-compiled loops, Fortran-compiled loops, and C-compiled loops. Normalize the CPU time for each mesh by the fastest version. Filename: ``wave3D_u0``. .. --- end exercise ---