.. _quickstart: ========== Quickstart ========== This section describes the plain vanilla Hermite-Taylor method for solving the scalar transport equation with periodic boundary conditions .. math:: :nowrap: \begin{eqnarray} && u_t = a u_x, \ \ x \in [x_L,x_R], \, t>0, \\ && u(x,0) = f(x), \ \ u(x_L,t) = u(x_R,t). \end{eqnarray} For this particular example we choose :math:`x_L = 0, \, x_R = 1, \, f(x) = \sin(2 \pi x)`. The implementation of this basic example has the same structure as the more involved examples discussed later. The main program is contained in the file `chides/transport/1d/herm1d.f90`_. The main program will call subroutines that are specific to the above problem and contained in the module ``pdemodule`` implemented in the file `chides/transport/1d/transport1d.f90`_. All examples will have these two elements, a main program (``hermXd.f90``, X being the number of dimensions) and a PDE module (``pdename.f90``). In addition, a typical implementation will call various ``service`` routines for polynomial manipulation, input-output, etc. For this first example we try to keep the complexity low and only use two service routines `chides/service/Hermite.f90`_ related to the central task of forming a Hermite interpolating polynomial, and `chides/service/get_Taylor1d.f90`_ to turn point values into Taylor polynomials. Grids +++++ The equations will be advanced one half time step :math:`\frac{\Delta t}{2}` from a primal grid .. math:: :nowrap: \begin{equation} x_i = x_L + i h_x, \ \ i = 0, \ldots, n_x, \ \ h_x = (x_R-x_L)/n_x, \end{equation} to a dual grid .. math:: :nowrap: \begin{equation} x_i = x_L + (i+1/2) h_x, \ \ i = 1, \ldots, n_x, \ \ h_x = (x_R-x_L)/n_x. \end{equation} In ``herm1d.f90`` this corresponds to the code block: .. code-block:: fortran hx=(xr-xl)/DBLE(nx) ! Primal grid DO i=0,nx x(i)=xl+hx*DBLE(i) END DO ! Dual grid DO i=1,nx xh(i)=xl+hx*(DBLE(i)-0.5d0) END DO Approximating the solution by a piecewise polynomial ++++++++++++++++++++++++++++++++++++++++++++++++++++ Hermite methods are based on the representation of the solution :math:`u(x,t_n)` at some fixed time level :math:`t_n` in terms of degree :math:`m+1` nodal based scaled polynomials centered around each grid point. For example on the primal grid we have that .. math:: :nowrap: \begin{equation} u(x,t_n) \approx p_i(x) = \sum_{l = 0}^{m} c_l \left( \frac{x - x_i }{h_x}\right)^l, x \in [x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}], \ \ i = 0,\ldots,n_x. \end{equation} Note that on the grid point :math:`x_i` the coefficients of the polynomial are related to the derivatives of the polynomial in the following way .. math:: :nowrap: \begin{equation} c_l = \frac{h_x^l}{l!} \left. \frac{d^l p}{dx^l} \right|_{x = x_i}. \end{equation} Thus, the degrees of freedom in the method can (interchangeably) be thought of as derivatives of the solution or coefficients of a local polynomial. In CHIDES the coefficients :math:`c_l` on each node on the primal and dual grid are stored in two arrays, again see ``herm1d.f90``: .. code-block:: fortran DOUBLE PRECISION, DIMENSION(0:m,nvar,0:nx) :: u DOUBLE PRECISION, DIMENSION(0:m,nvar,1:nx) :: uh The leading dimension refers to the number of coefficients (or derivatives) at each node, `nvar` is the number of variables, e.g. 1 for a scalar problem and >1 for a system of equations. The third dimension refers to the index along the grids. Initial data +++++++++++++ At the beginning of the computation the initial data is used to initialize the coefficients :math:`c_l` by call ``CALL indat(u,x,hx,nx,m,nvar)``, with ``indat`` being specified in ``transport1d.f90``. Here, as the initial data is a sin function computing derivatives is rather easy, but instead we use the service routine ``get_Taylor1d.f90`` which computes the coefficients using polynomial interpolation: .. code-block:: fortran u = 0.d0 DO j=0,nx ! CALL get_Taylor1d(u(:,1,j),x(j),m,2*m+1,hx,infun) ! END DO Here the inputs to get_Taylor1d are the local coordinate ``x(j)``, the degree of the Taylor polynomial ``m``, the degree of the interpolant ``2*m+1`` chosen to guarantee convergence at the design order, the cell width ``hx``, and finally the function evaluating the initial data ``infun``. In this case .. code-block:: fortran ! infun=SIN(freq*pi*x) ! The output is the array of scaled Taylor coefficients at the jth node. The interpolation ``get_Taylor1d.f90`` is carried out on a Chebyshev grid of width ``hx`` centered at ``x(j)``. See the description in Chapter 4 for details. 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}}` is to form a degree :math:`2m+1` polynomial .. math:: :nowrap: \begin{equation} p_{i-\frac{1}{2}}(x) = \sum_{l = 0}^{2m+1} c_l \left( \frac{x - x_{i-\frac{1}{2} }}{h_x}\right)^l,\ \ x \in [x_{i-1},x_{i}], \ \ i = 1,\ldots,n_x, \end{equation} from the polynomials representing the initial data at the two nearby primal gridpoints. That is, the coefficients of the new polynomial are determined by insisting that it interpolates the nearby polynomials. Functionally, this amounts to applying a the linear map corresponding to Hermite interpolation, say :math:`A`, to the nearby polynomials to find :math:`p_{i-\frac{1}{2}} = A(p_{i-1},p_{i})`. In the actual code this just corresponds to assembling the nearby coefficients into a :math:`2m+2 \times 1` vector which is multiplied by a matrix ``Aintx`` corresponding to :math:`A`. For example in the section of the code where the initial data is output to file we have: .. code-block:: fortran DOUBLE PRECISION, DIMENSION(0:2*m+1,0:2*m+1) :: Aintx DOUBLE PRECISION, DIMENSION(0:2*m+1) :: udata, pplot . . udata(0:m) = u(0:m,1,i) udata(m+1:2*m+1) = u(0:m,1,i+1) pplot = MATMUL(Aintx,udata) The matrix ``Aintx`` is computed in the program ``herm1d.f90`` .. code-block:: fortran CALL Hermite(m,node_left,node_right,node_interp,Aintx,0) using the subroutine ``Hermite`` from in `chides/service/Hermite.f90`_. Evolving the Hermite interpolating polynomial +++++++++++++++++++++++++++++++++++++++++++++ To evolve the above polynomial we consider the space-time approximation .. math:: :nowrap: \begin{equation} p^n_{i-\frac{1}{2}}(x,t) = \sum_{l=0}^{2m+1} \sum_{s=0}^{q} d_{ls} \left( \frac{x - x_{i+\frac{1}{2}} }{h_x}\right)^l \left( \frac{t - t_n }{\Delta t}\right)^s, \end{equation} which when evaluated at time :math:`t = t_n` reduces to .. math:: :nowrap: \begin{equation} p_{i-\frac{1}{2}}(x,t_n) = \sum_{l = 0}^{2m+1} d_{l0} \left( \frac{x - x_{i-\frac{1}{2} }}{h_x}\right)^l,\ \ x \in [x_{i-1},x_{i}], \ \ i = 1,\ldots,n_x, \end{equation} and we thus have :math:`d_{l0} = c_l`. To find the remaining coefficients :math:`d_{ls}, s > 0` we will use the PDE. In particular we will insist that the approximation :math:`p_{i-\frac{1}{2}}(x,t)` satisfies the PDE at :math:`(x^\ast,t_n)`. For notational convenience consider the approximation at one dual gridpoint, say :math:`x^\ast` and denote it by :math:`p(x,t)`. Then we have that .. math:: :nowrap: \begin{eqnarray} \frac{\partial p(x^\ast,t_n)}{\partial x} = \frac{d_{10}}{h_x},\\ \frac{\partial p(x^\ast,t_n)}{\partial t} = \frac{d_{01}}{\Delta t}, \end{eqnarray} thus from the PDE we find that we may express the first time-derivative in the Taylor series .. math:: :nowrap: \begin{equation} \frac{\partial p(x^\ast,t_n)}{\partial t} = a \frac{\partial p(x^\ast,t_n)}{\partial x} \ \ \Rightarrow \ \ d_{01} = a \frac{{\Delta t}}{h_x} d_{10}. \end{equation} By repeated differentiation in space and time of the PDE .. math:: :nowrap: \begin{equation} \frac{\partial^{r+k} u}{\partial t^r \partial^k} = a \frac{\partial^{r-1+k+1} u}{\partial t^{r-1} \partial x^{k+1}}, \end{equation} which upon replacing :math:`u` with :math:`p` yields the recursion relation: .. math:: :nowrap: \begin{equation} d_{rk} = a \frac{{\Delta t}}{h_x} \frac{k+1}{r} d_{r-1,k+1},, \ \ r = 1,\ldots, \ \ k = 0,\ldots, \end{equation} specifying all time derivatives in the Taylor series in terms of spatial derivatives. In the code this update is done in the main time-stepping loop .. code-block:: fortran DO it=1,nsteps ! first half step t = DBLE(it-1)*dt ! At the end of the first half step all entries of ! uh will have been updated uh = 0.d0 ! Loop over all dual gridpoints DO i=1,nx ! tcofs contains the coefficients c_{l,s} ! l = 0,2*m+1, s = 0,q tcofs=0.d0 ! c_{l,0}, l = 0,2*m+1 are found by interpolating data ! from adjacent primal nodes udata(0:m) = u(:,1,i-1) udata(m+1:2*m+1) = u(:,1,i) tcofs(0:2*m+1,0,1)=MATMUL(Aintx,udata) ! The rest of the coefficients are found via the ! recursion: ! c_{l,s} = dt/h*(l+1)/s * c_{l+1,s-1} ! This is done in the subroutine pde CALL pde(tcofs,2*m+1,q,nvar,cflx,etax,xh(i),t) Filling in all the ''time coefficients'' in ``tcofs`` is done inside the subroutine ``pde`` from `chides/transport/1d/herm1d.f90`_: .. code-block:: fortran SUBROUTINE pde(tcofs,mp,q,nvar,cflx,etax,x,t) ! ! This pde implements the recursion relation ! associated with the linear transport equation ! u_t = a*u_x ! implicit none INTEGER, INTENT(IN) :: mp,q,nvar DOUBLE PRECISION, DIMENSION(0:mp,0:q,nvar), INTENT(INOUT) :: tcofs DOUBLE PRECISION, INTENT(IN) :: cflx,etax,x,t INTEGER :: idx,idt DOUBLE PRECISION :: a_dt_h ! tcofs(:,1,:)=0.d0 a_dt_h = a*cflx DO idt=1,q DO idx=0,mp-idt tcofs(idx,idt,1)=a_dt_h*DBLE(idx+1)/DBLE(idt)*tcofs(idx+1,idt-1,1) END DO END DO END SUBROUTINE pde With all the time-coefficients available we can evaluate the approximation :math:`p(x,t)` at a later time :math:`(x^\ast,t_n+\frac{\Delta t}{2})` for the first :math:`m` spatial derivatives yielding the new approximation of :math:`u(x,t_{n+\frac{1}{2}})`: .. math:: :nowrap: \begin{equation} u(x,t_{n+\frac{1}{2}}) \approx p^n_{i-\frac{1}{2}}(x) = \sum_{l=0}^{m} \sum_{s=0}^{q} \tilde{c}_{l} \left( \frac{x - x_{i+\frac{1}{2}} }{h_x}\right)^l, \ \ x \in [x_{i-1},x_{i}], \ \ i = 1,\ldots,n_x, \end{equation} where the new coefficients :math:`\tilde{c}_l` are .. math:: :nowrap: \begin{equation} \tilde{c}_l = \sum_{s = 0}^{q} d_{ls} \left( \frac{1}{2} \right)^s, \ \ l =0,\ldots,m. \end{equation} In the code this looks like: .. code-block:: fortran ! When all coefficients are known ! the Taylor series can be evaluated ! and uh(i) updated DO idx=0,m uh(idx,1,i) = tcofs(idx,0,1) dfact = 1.d0 DO l = 1,q dfact = dfact*0.5d0 uh(idx,1,i) = uh(idx,1,i) & + dfact*tcofs(idx,l,1) END DO END DO Second half timestep and boundary conditions +++++++++++++++++++++++++++++++++++++++++++++ The second half timestep follows the same procedure as described above: 1. Form an interpolating polynomial .. code-block:: fortran ! second half step ! Same thing as above but we now update ! the interor of u t = t + 0.5d0*dt DO i = 1,nx-1 tcofs = 0.d0 udata(0:m) = uh(:,1,i) udata(m+1:2*m+1) = uh(:,1,i+1) tcofs(0:2*m+1,0,1) = MATMUL(Aintx,udata) 2. Then use the PDE and update the approximation at the full timestep .. code-block:: fortran CALL pde(tcofs,2*m+1,q,nvar,cflx,etax,x(i),t) ! Add up the Taylor series DO idx=0,m u(idx,1,i)=tcofs(idx,0,1) dfact=1.d0 DO l=1,q dfact=dfact*.5d0 u(idx,1,i)=u(idx,1,i) & +dfact*tcofs(idx,l,1) END DO END DO END DO For this example the boundary conditions are taken to be periodic. This means that the interior update procedure can be applied on the boundary as well, with the slight difference that when the interpolating polynomial is formed, data ''from outside the domain'' is replaced by the periodic copy: .. code-block:: fortran ! now the left and right bcs ! left tcofs=0.d0 udata(0:m)=uh(:,1,nx) udata(m+1:2*m+1)=uh(:,1,1) tcofs(0:2*m+1,0,1)=MATMUL(Aintx,udata) CALL pde(tcofs,2*m+1,q,nvar,cflx,etax,x(0),t) ! Add up the Taylor series DO idx=0,m u(idx,1,0)=tcofs(idx,0,1) dfact=1.d0 DO l=1,q dfact=dfact*.5d0 u(idx,1,0)=u(idx,1,0) & +dfact*tcofs(idx,l,1) END DO END DO ! Right u(:,:,nx)=u(:,:,0) Compiling and Running the program +++++++++++++++++++++++++++++++++ In the directory ``chides/transport/1D`` 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/transport/1D`` 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 ``herm1d.x`` that you can run. Displaying the results ++++++++++++++++++++++ After running the program you should see something like: .. code-block:: none bash-3.2$ ./herm1d.x Space step, Max-error 5.00000E-02 6.60343E-10 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_sol.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 ``transport1d.f90``): .. code-block:: fortran MODULE pdemodule ! ! ! Problem setup for solving 1d transport equation: ! u_t = a u_x ! With periodic boundary conditions ! ! IMPLICIT NONE INTEGER, PARAMETER :: m = 3 ! spatial order 2m+1 INTEGER, PARAMETER :: q = 2*m+2 ! terms in the Taylor series INTEGER, PARAMETER :: nx = 20 ! Discretization uses nx+1 points from xl to xr INTEGER, PARAMETER :: nxp = 500 ! Number of points in the output DOUBLE PRECISION, PARAMETER :: CFL = 0.95d0 ! Ratio dt/h INTEGER, PARAMETER :: nout = 1 ! output the solution every nout times INTEGER, PARAMETER :: nvar = 1 ! DOUBLE PRECISION, PARAMETER :: a = 1.0d0 ! Wave speed DOUBLE PRECISION, PARAMETER :: ttot = 10.0d0 ! Final time DOUBLE PRECISION, PARAMETER :: xl = 0.0d0 ! Domain boundaries DOUBLE PRECISION, PARAMETER :: xr = 1.0d0 ! ! Other constants DOUBLE PRECISION, PARAMETER :: pi = 3.1415926535897932385d0 DOUBLE PRECISION, PARAMETER :: freq = 2.d0 ! CONTAINS This is where you would change the order and the number of gridpoints etc. .. _chides/transport/1d/herm1d.f90: https://bitbucket.org/appelo/chides .. _chides/transport/1d/transport1d.f90: https://bitbucket.org/appelo/chides .. _chides/service/Hermite.f90: https://bitbucket.org/appelo/chides .. _chides/service/get_Taylor1d.f90: https://bitbucket.org/appelo/chides