.. _acoustic2d:
----------------------------------------
The acoustic equations in two dimensions
----------------------------------------
This example illustrates how to implement a two dimensional version of
the Hermite-Taylor method (see :ref:`quickstart` for an example in
1D) for solving a linear first order hyperbolic system. The source files reside in ``chides/acoustic/2D``
We will consider the equations of acoustics which can be derived from the compressible
Navier-Stokes equations. Assuming an isentropic, inviscid, quiscent flow and small
perturbations the conservation of mass and momentum are reduced to
.. math::
:nowrap:
\begin{eqnarray}
\frac{\partial \rho}{\partial t} + \rho_0 \nabla \cdot \vec{u} = 0,\\
\frac{\partial \vec{u}}{\partial t} + \frac{c_0^2}{\rho_0} \nabla \rho = 0.
\end{eqnarray}
Here we will take :math:`\rho_0 = c_0 = 1`.
To make the equations look more like the implementation we write the
above system in terms of the density and the velocities :math:`u` and
:math:`v`
.. math::
:nowrap:
\begin{eqnarray}
&&\frac{\partial \rho}{\partial t} + \frac{\partial u}{\partial x}+ \frac{\partial v}{\partial y} = 0,\\
&&\frac{\partial u}{\partial t} + \frac{\partial \rho}{\partial x}
= 0, \\
&&\frac{\partial v}{\partial t} + \frac{\partial \rho}{\partial y} = 0.
\end{eqnarray}
A first example in two dimensions with periodic boundary conditions
-------------------------------------------------------------------
As a first example of how CHIDES can be used to construct multi-dimensional solvers we consider the above equations in a periodic box
:math:`(x,y) \in [-\pi,\pi]^2`. Then an exact solution is
.. math::
:nowrap:
\begin{eqnarray}
p(x,y,t) &=& \sin(x) \sin(y) \sin(\sqrt{2} t), \\
u(x,y,t) &=& \frac{1}{\sqrt{2}} \cos(x) \sin(y) \cos(\sqrt{2} t),\\
v(x,y,t) &=& \frac{1}{\sqrt{2}} \sin(x) \cos(y) \cos(\sqrt{2} t).
\end{eqnarray}
The code structure is similar to the one-dimensional case. Here the
main program is contained in the file
`chides/acoustic/2D/herm2d_per.f90`_ and the module describing the
PDE, initial data and discretization related parameters resides in the
file `chides/acoustic/2D/acoustic_per.f90`_.
.. _chides/acoustic/2D/herm2d_per.f90: https://bitbucket.org/appelo/chides/src
.. _chides/acoustic/2D/acoustic_per.f90: https://bitbucket.org/appelo/chides/src
Grids
+++++
As in the one-dimensional case the equations will be advanced one half time step :math:`\frac{\Delta t}{2}` from a primal grid
.. math::
:nowrap:
\begin{eqnarray}
&& x_i = x_L + i h_x, \ \ i = 0, \ldots, n_x, \ \ h_x = (x_R-x_L)/n_x,\\
&& x_j = y_B + j h_y, \ \ j = 0, \ldots, n_y, \ \ h_y = (y_T-y_B)/n_y,
\end{eqnarray}
to a dual grid
.. math::
:nowrap:
\begin{eqnarray}
x_i = x_L + (i+1/2) h_x, \ \ i = 1, \ldots, n_x, \ \ h_x = (x_R-x_L)/n_x,\\
y_j = y_B + (j+1/2) h_y, \ \ j = 1, \ldots, n_y, \ \ h_y = (y_T-y_B)/n_y.
\end{eqnarray}
In ``herm2d_per.f90`` this corresponds to the code block:
.. code-block:: fortran
hx=(xr-xl)/DBLE(nx)
DO i=0,nx
x(i)=xl+hx*DBLE(i)
IF (i > 0) THEN
xh(i)=xl+hx*(DBLE(i)-.5d0)
END IF
END DO
hy=(yt-yb)/DBLE(ny)
DO i=0,ny
y(i)=yb+hy*DBLE(i)
IF (i > 0) THEN
yh(i)=yb+hy*(DBLE(i)-.5d0)
END IF
END DO
Approximating the solution by a piecewise polynomial
++++++++++++++++++++++++++++++++++++++++++++++++++++
In two dimensions we represent the solution of one of the fields
(the density and the velocities in the x and y dierction) by a
polynomial. For example for the velocity in the x direction we have :math:`u(x,y,t_n)` at a
fixed time level :math:`t_n` in terms of degree :math:`m+1` nodal
based scaled tensor product polynomials centered around each grid point. For example on the primal grid we have that
.. math::
:nowrap:
\begin{equation}
u(x,y,t_n) \approx p_{ij}(x,y) = \sum_{k = 0}^{m} \sum_{l = 0}^{m}
c_{kl} \left( \frac{x - x_i }{h_x}\right)^k
\left( \frac{y - y_j }{h_y}\right)^l.
\end{equation}
Again, at a grid point :math:`(x_i,y_j)` the coefficients of the
polynomial can be expressed in terms of derivatives
.. math::
:nowrap:
\begin{equation}
c_{kl} = \frac{h_x^k}{k!} \frac{h_y^l}{l!}\left. \frac{d^{k+l} p}{dx^k dy^l} \right|_{x = x_i,y = y_j}.
\end{equation}
In CHIDES the coefficients :math:`c_{kl}` on each node on the primal and dual grid are stored in two arrays, see ``herm2d_per.f90``:
.. code-block:: fortran
DOUBLE PRECISION, DIMENSION(0:m,0:m,nvar,0:nx,0:ny) :: u
DOUBLE PRECISION, DIMENSION(0:m,0:m,nvar,0:nx+1,0:ny+1) :: uh
The two leading dimension refers to the number of coefficients (or
derivatives) at each node, `nvar` is the number of variables, e.g. 3
for the acoustic system. The last two dimensions refers to the indicies along the grids.
Initial data
++++++++++++
At the beginning of the computation the initial data is used to
initialize the coefficients :math:`c_{kl}` by the call ``call
indat(u,x,y,t,hx,hy,nx,ny,m,nvar)``, with ``indat`` being specified in
``acoustic2d_per.f90``. Here, as the initial data are composed of
trigonometric functions we use the same approach as described in
:ref:`quickstart`.
Forming the Hermite interpolating polynomial
++++++++++++++++++++++++++++++++++++++++++++
The first step in evolving the approximate solution to the next half
timestep at a dual gridpoint
:math:`(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})` is to form a degree :math:`2m+1` tensor product polynomial
.. math::
:nowrap:
\begin{equation}
p_{i-\frac{1}{2},j-\frac{1}{2}}(x) = \sum_{k = 0}^{2m+1} \sum_{l = 0}^{2m+1}
c_{kl} \left( \frac{x - x_{i-\frac{1}{2} } }{h_x}\right)^k
\left( \frac{y - y_{j-\frac{1}{2} }}{h_y}\right)^l.
\end{equation}
from the polynomials representing the initial data at the four nearby
primal gridpoints. That is, the coefficients of the new polynomial are
determined by insisting that it interpolates the nearby
polynomials. For each variable we form the Hermite interpolant at the dual gridpoint ``(xh(i),yh(j))``, i.e. at the
center of the cell with corners ``(x(i-1),y(j))``, ``(x(i),y(j))``, ``(x(i-1),y(j-1))``, ``(x(i),y(j-1))``.
As the approximation is a tensor product polynomial we first interpolate on the top and bottom of the cell
to the two face midpoints. Then we interpolate from the to the cell
center. In the code this is done for one cell at a time by a simple
call to ``Taylor_to_Hermite2d``
.. code-block:: fortran
DO j=1,ny
DO i=1,nx
tcofs = 0.d0
call Taylor_to_Hermite2d(uhermite,u(:,:,:,i-1:i,j-1:j),cintx,cinty,m,nvar)
! We store the Hermite interpolant in the array tcofs. The third index is for the
! time derivatives of the variables. The Hermite interpolant is the 0th time-
! derivative.
do ivar = 1,nvar
tcofs(0:2*m+1,0:2*m+1,0,ivar) = uhermite(:,:,ivar)
end do
....
Evolving the Hermite interpolating polynomial
++++++++++++++++++++++++++++++++++++++++++++++
To evolve in time we first expand the solution in a Taylor series in
time:
.. math::
:nowrap:
\begin{equation}
p^n_{i-\frac{1}{2},j-\frac{1}{2}}(x) = \sum_{k = 0}^{2m+1} \sum_{l =
0}^{2m+1} \sum_{s=0}^{q}
c_{kls} \left( \frac{x - x_{i-\frac{1}{2} } }{h_x}\right)^k
\left( \frac{y - y_{j-\frac{1}{2} }}{h_y}\right)^l \left( \frac{t - t_n }{\Delta t}\right)^s.
\end{equation}
The evolution in time requires that we know all coefficients
:math:`c_{kls}`. In the code these are computed by a call to
the routine ``pde``. This routine fills in the time
derivatives :math:`c_{kls}, \, s =1,\ldots,q` and returns them in the array ``tcofs(:,:,1:q,:)``.
.. code-block:: fortran
....
! The rest of the time-derivatives are filled in using the PDE.
! Here the PDE is autonomous so we do not need to pass in coordinates.
CALL pde(tcofs,2*m+1,q,nvar,cflx,cfly)
Once the time derivatives are known we compute the Taylor series at
the next half timestep and store it in ``uh``.
.. code-block:: fortran
....
! With all the time derivatives we finally add up the Taylor series
! for our degrees of freedom, i.e. derivatives up to degree m.
uh(0:m,0:m,1:nvar,i,j) = tcofs(0:m,0:m,0,1:nvar)
dfact=1.d0
DO l=1,q
dfact=dfact*.5d0
uh(0:m,0:m,1:nvar,i,j) = uh(0:m,0:m,1:nvar,i,j) &
+dfact*tcofs(0:m,0:m,l,1:nvar)
END DO
END DO
END DO
The ``pde`` routine residing in ``acoustic2d_per.f90`` has the same structure as in the one dimensional case
.. code-block:: fortran
SUBROUTINE pde(tcofs,mp,q,nvar,cflx,cfly)
!
! time recursion coefficients as determined by the pde
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nvar,mp,q
DOUBLE PRECISION, DIMENSION(0:mp,0:mp,0:q,nvar), INTENT(INOUT) :: tcofs
DOUBLE PRECISION, INTENT(IN) :: cflx,cfly
INTEGER :: idx,idy,idt
DOUBLE PRECISION, DIMENSION(0:mp,0:mp) :: ux,px,vy,py
!
DO idt=1,q
! Compute the "x-derivatives"
DO idx=0,mp-1
px(idx,:) = cflx*DBLE(idx+1)*tcofs(idx+1,:,idt-1,1)
ux(idx,:) = cflx*DBLE(idx+1)*tcofs(idx+1,:,idt-1,2)
END DO
px(mp,:) = 0.d0
ux(mp,:) = 0.d0
! Compute the "y-derivatives"
DO idy=0,mp-1
py(:,idy) = cfly*DBLE(idy+1)*tcofs(:,idy+1,idt-1,1)
vy(:,idy) = cfly*DBLE(idy+1)*tcofs(:,idy+1,idt-1,3)
END DO
py(:,mp) = 0.d0
vy(:,mp) = 0.d0
!
tcofs(:,:,idt,1) = -ux - vy
tcofs(:,:,idt,2) = -px
tcofs(:,:,idt,3) = -py
tcofs(:,:,idt,:) = tcofs(:,:,idt,:)/DBLE(idt)
END DO
END SUBROUTINE pde
The main difference here is that the equations are coupled and the
recursion relations needs to be evaluated one time derivative at a
time. The derivation of the recursion is analogous to the one
dimensional case (described in detail in :ref:`quickstart`).
Second half timestep and boundary conditions
++++++++++++++++++++++++++++++++++++++++++++
The second half timestep first enforces the periodic boundary
conditions and then evolves the solution in the same way as above:
1. Enforce periodic boundary conditions
.. code-block:: fortran
! second half step
t = t + 0.5d0*dt
! We first enforce periodic boundary conditions by copying
! the periodic values to the "ghost points"
uh(:,:,:,0,:) = uh(:,:,:,nx,:)
uh(:,:,:,nx+1,:) = uh(:,:,:,1,:)
uh(:,:,:,:,0) = uh(:,:,:,:,ny)
uh(:,:,:,:,ny+1) = uh(:,:,:,:,1)
2. Form an interpolating polynomial
.. code-block:: fortran
! Then we proceed as in the first half step.
DO j=0,ny
DO i=0,nx
tcofs=0.d0
! Form the interpolant at a primal gridpoint
! (x(i),y(j))
call Taylor_to_Hermite2d(uhermite,uh(0:m,0:m,1:nvar,i:i+1,j:j+1),cintx,cinty,m,nvar)
do ivar = 1,nvar
tcofs(0:2*m+1,0:2*m+1,0,ivar) = uhermite(:,:,ivar)
end do
...
3. Then use the PDE and update the approximation at the full timestep
.. code-block:: fortran
...
! Compute time derivatives
CALL pde(tcofs,2*m+1,q,nvar,cflx,cfly)
! Evaluate the Taylor series and update the solution
! at a primal gridpoint.
u(0:m,0:m,1:nvar,i,j)=tcofs(0:m,0:m,0,1:nvar)
dfact=1.d0
DO l=1,q
dfact=dfact*.5d0
u(0:m,0:m,1:nvar,i,j)=u(0:m,0:m,1:nvar,i,j) &
+dfact*tcofs(0:m,0:m,l,1:nvar)
END DO
END DO
END DO
Compiling and Running the program
+++++++++++++++++++++++++++++++++
In the directory ``chides/acoustic/2D`` there is a simple ``makefile`` that can be used to compile the program by simply typing ``make`` in a terminal window. This assumies you are standing in the directory ``chides/acoustic/2D`` and that your system has ``gfortran`` installed. If you are using another compiler simply change the makefile.
If you are successful in building the solver you should have an executable named ``periodic.x`` that you can run.
Displaying the results
++++++++++++++++++++++
After running the program you should see something like:
.. code-block:: none
bash-3.2$ ./periodic.x
Computation finished:
-------------------------------------------------------------------------
m dt hx maxerr l2err-p l2err-u l2err-v
-------------------------------------------------------------------------
05 0.1667E+01 0.2094E+01 0.4164E-07 0.1215E-06 0.1471E-07 0.1471E-07
-------------------------------------------------------------------------
The program should also have generated a number of snapshots of the solution, ``u0000000.txt``, ``u0000001.txt``, etc. These can be displayed in Matlab using the script ``plot_sol2D.m``.
Changing the order and number of gridpoints
+++++++++++++++++++++++++++++++++++++++++++
All the parameters associated with a simulation are contained in the
top block of the ``pdemodule`` (residing in ``acoustid2d_per.f90``),
this is where you would change the order and the number of gridpoints
etc.
**NOTE** that as the makefile structure is as simple as possible you
should do ``make clean; make`` when you change the parameters. This
ensures that all dependencies between modules are up to date.
.. code-block:: fortran
MODULE pdemodule
!
! This module solves the pde
! rho_t + u_x + v_y = 0,
! u_t + p_x = 0,
! v_t + v_y = 0,
!
! On a domain (x,y) \in [-\pi,pi]^2,
! with periodic boundary conditions.
!
! A solution is, p(x,y,t) = \sin(x) \sin(y) \sin(\omega t), with \omega = \sqrt(2),
! u(x,y,t) = 1/\omega \cos(x) \sin(y) \cos(\omega t),
! v(x,y,t) = 1/\omega \sin(x) \cos(y) \cos(\omega t).
!
!
IMPLICIT NONE
INTEGER, PARAMETER :: m = 5 ! m is the number of derivatives on each node
! the spatial order of the method is 2m+1
INTEGER, PARAMETER :: q = 4*m+2 ! terms in the Taylor series
!
DOUBLE PRECISION, PARAMETER :: pi = 3.1415926535897932385d0
DOUBLE PRECISION, PARAMETER :: xl = -pi
DOUBLE PRECISION, PARAMETER :: xr = pi
DOUBLE PRECISION, PARAMETER :: yb = -pi
DOUBLE PRECISION, PARAMETER :: yt = pi
DOUBLE PRECISION, PARAMETER :: ttot = 10.d0
INTEGER, PARAMETER :: nx = 3 ! Discretization uses nx+1 points from xl to xr
INTEGER, PARAMETER :: ny = nx
DOUBLE PRECISION, PARAMETER :: CFL = 0.95d0 ! Ratio dt/h
!
INTEGER, PARAMETER :: nref = 50 ! refinement factor for i-o and error computation.
!
INTEGER, PARAMETER :: nout = 1 ! output the solution every nout times
INTEGER, PARAMETER :: nvar = 3 ! Number of variables in solution
An example with boundary conditions
-------------------------------------------------------------------
We next consider a case were we impose boundary conditions. This example
is described in the files ``acoustic2d_mirror.f90`` and ``herm2d_mirror.f90``. The domain
is still a rectangle :math:`(x,y) \in [x_L,x_R] \times [y_B,y_T]` but
now we impose zero Dirichlet boundary conditions for the
density. Differentiating the boundary condition in time and in the
tangential directions together with the governing equations yields the
boundary conditions:
.. math::
:nowrap:
\begin{eqnarray}
&& \rho = 0, \ \ x={x_L,x_R}, y={y_L,y_T},\\
&& u = 0, \ \ v_y = 0, \ \ y={y_L,y_T}, \\
&& v = 0, \ \ u_x = 0, \ \ x={x_L,x_R}.
\end{eqnarray}
The boundary conditions are imposed at the half level by filling in
the ghost point values by mirroring the polynomials so that they
are either odd around the boundary and thus agree with the homogenous
Dirichlet boundary conditions or so that they are even and satisfy
homogenous Neumann conditions. The update of the boundary conditions
is done in ``herm2d_mirror.f90`` by calling the subroutine below:
.. code-block:: fortran
SUBROUTINE update_bc(uh,nx,ny,m,nvar)
implicit none
INTEGER, INTENT(IN) :: nx,ny,m,nvar
DOUBLE PRECISION, DIMENSION(0:m,0:m,nvar,0:nx+1,0:ny+1), INTENT(INOUT) :: uh
DOUBLE PRECISION, DIMENSION(0:m) :: DIR,NEU
INTEGER :: i,j,idx,idy
do i = 0,m
NEU(i) = (-1.d0)**i
DIR(i) = (-1.d0)**(i+1)
end do
NEU(0) = 1.d0
! The pde is
! (1) rho_t = -u_x - v_y,
! (2) u_t = -rho_x,
! (3) v_t = -rho_y.
! We implement rho = 0 on all boundaries.
! By (2) this implies u = 0 on top and bottom.
! By (1)+(2) we have v_y = 0 on top and bottom.
! Also:
! By (3) this implies v = 0 on left and right.
! By (1)+(3) we have u_x = 0 on left and right.
! Left and right side:
do j = 1,ny
do idy = 0,m
uh(0:m,idy,1,0,j) = DIR*uh(0:m,idy,1,1,j)
uh(0:m,idy,1,nx+1,j) = DIR*uh(0:m,idy,1,nx,j)
uh(0:m,idy,2,0,j) = NEU*uh(0:m,idy,2,1,j)
uh(0:m,idy,2,nx+1,j) = NEU*uh(0:m,idy,2,nx,j)
uh(0:m,idy,3,0,j) = DIR*uh(0:m,idy,3,1,j)
uh(0:m,idy,3,nx+1,j) = DIR*uh(0:m,idy,3,nx,j)
end do
end do
! Top and bottom
do i = 0,nx+1
do idx = 0,m
uh(idx,0:m,1,i,0) = DIR*uh(idx,0:m,1,i,1)
uh(idx,0:m,1,i,ny+1) = DIR*uh(idx,0:m,1,i,ny)
uh(idx,0:m,2,i,0) = DIR*uh(idx,0:m,2,i,1)
uh(idx,0:m,2,i,ny+1) = DIR*uh(idx,0:m,2,i,ny)
uh(idx,0:m,3,i,0) = NEU*uh(idx,0:m,3,i,1)
uh(idx,0:m,3,i,ny+1) = NEU*uh(idx,0:m,3,i,ny)
end do
end do
An example in a duct with forcing and Taylor series time stepping
-----------------------------------------------------------------
This example is contained in the files ``herm2d_duct.f90`` and ``acoustic2d_duct.f90``.
.. code-block:: fortran
! This module solves the pde
! p_t = u_x + v_y + f,
! u_t = p_x,
! v_t = p_y,
!
! On a domain (x,y) \in [-10,10] X [0,1]
!
! With wall boundary conditions on y=0,1 and periodic
! boundary conditions in x
!