.. _burgers:
------------------------------------------
Viscous Burgers' equation in one dimension
------------------------------------------
Consider the non-linear conservation
law:
.. math::
:nowrap:
u_t+f(u)_x=0,
with the viscous Burgers' equation as a special case with :math:`f(x)
= \frac{u^2}{2}-\nu u_x`. As before we will take derivatives in space
and time of the equation and insist that the resulting relations hold
for the approximation. Typically the flux function will be non-linear
and we must find an polynomial approximation of it using the
polynomial manipulation routines that are part of CHIDES. As the
non-linearity depends on time it becomes less favorable to use a
Taylor series expansion in time and instead we opt to use a one-step
method to evolve the solution and its derivatives. Among other things this example illustrate:
1. How to use a one-step time-integrators, for example a Runge-Kutta
method together with Hermite interpolation.
2. How to improve the accuracy and time-step restrictions
of the method by taking multiple sub-timesteps per half timestep.
3. How to discretize second order derivatives.
4. Semi-conservative and non-conservative discretizations of conservation laws.
The source files reside in ``chides/burgers/1D``
To find an approximation to the solution of the above conservation law
we assume we have interpolated two degree :math:`m` polynomials at two
adjacent gridpoints on the primal grid. That is, at a dual gridpoint
:math:`x_{i+\frac{1}{2}}` is to we start with the degree :math:`2m+1` polynomial
.. math::
:nowrap:
\begin{equation}
p_{i-\frac{1}{2}}(x,t)= \sum_{l = 0}^{2m+1} c_l(t) \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}
Note that we now assume that the coefficients :math:`c_l(t)` depend on
time so that if we compute one time derivative and :math:`l`
derivatives in space and evaluate at
:math:`x = x_{i-\frac{1}{2}}` we get terms of the kind
:math:`c^\prime_l(t)` that we can be identified to be scaled
approximations of :math:`\frac{\partial^l u_t(0,t)}{\partial x^l}`.
Evaluating the non-linearity
++++++++++++++++++++++++++++
As we obtained the time derivatives :math:`c^\prime_l(t)` by
differentiating the approximation it is natural evolve them according
to the system of ODE that is obtained by differentiating the
conservation law in space and insisting that the approximation
satisfies those equations at :math:`x = x_{i-\frac{1}{2}}`. At first
this may appear as a cumbersome proposition as we must compute
multiple derivatives of :math:`f(u)_x`. However, if we consider the
flux to be approximated by a polynomial:
.. math::
:nowrap:
\begin{equation}
f(u) \approx \sum_{l = 0}^{2m+1} a_l(t) \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}
we can immediately see that the system of ODE governing :math:`c_l(t)`
is:
.. math::
:nowrap:
\begin{eqnarray}
c^\prime_l(t) & = & -\frac{(l+1)}{h_x} a_{l+1}(t),\ \ l = 0,\ldots,2m, \\
c^\prime_{2m+1}(t) & = & 0.
\end{eqnarray}
For the particular example we consider here we have that :math:`f(u) =
\frac{u^2}{2}-\nu u_x` which is computed in the routine ``compute_flux(u,ut,mp,nvar,hx)``
.. code-block:: fortran
SUBROUTINE compute_flux(u,flx,mp,nvar,hx)
!
! compute (f(u)-nu*u_x)
!
implicit none
INTEGER, INTENT(IN) :: mp,nvar
DOUBLE PRECISION, DIMENSION(0:mp,nvar), INTENT(IN) :: u
DOUBLE PRECISION, DIMENSION(0:mp,nvar), INTENT(OUT) :: flx
DOUBLE PRECISION, INTENT(IN) :: hx
INTEGER :: idx
DOUBLE PRECISION, DIMENSION(0:mp) :: q1,p1
!
! Compute the hyperbolic part of the flux
CALL tpolymul(u(:,1),u(:,1),p1,mp)
flx(:,1) = 0.5d0*p1
! Then add the viscous part
DO idx=0,mp-1
q1(idx) = dble(idx+1)*u(idx+1,1)
END DO
q1(mp) = 0.d0
flx(:,1) = flx(:,1) - (nu/hx)*q1
!
END SUBROUTINE compute_flux
In the subroutine we see the use of ``tpolymul(u(:,1),u(:,1),p1,mp)``
which is used to compute the square of the polynomial approximation.
After having found the flux it is straightforward to find
(:math:`\Delta t` times) the coefficients :math:`a_l` by
differentiating the flux as is done in the subroutine ``pde(fp,flx,mp,nvar,cflx)``
.. code-block:: fortran
fp = 0.d0
DO idx=0,mp-1
fp(idx,1) = -cflx*DBLE(idx+1)*flx(idx+1,1)
END DO
The output of this subroutine can then be used to evolve the
interpolant :math:`p(x,t)` over half a timestep.
The driver ``herm1d_non_cons``
++++++++++++++++++++++++++++++++
An implementation based on the above description can be found in
``/chides/burgers/1d/herm1d_non_cons.f90``. As before the basic steps
of the algorithm are:
1. Setup the initial data
.. code-block:: fortran
CALL indat(u,x,hx,nx,m,nvar)
2. Form the Hermite interpolant.
.. code-block:: fortran
DO i=1,nx
! Interpolate solution and fluxes to the cell center
call Taylor_to_Hermite1d(uhermite,u(:,:,i-1:i),Aintx,m,nvar)
...
3. And evolution of ``uhermite`` using a Runge-Kutta method
.. code-block:: fortran
call Taylor_to_Hermite1d(uhermite,u(:,:,i-1:i),Aintx,m,nvar)
do irk_sub = 1,rksub
urk = uhermite
do irk=1,rkstages ! do a local RK loop
urk = uhermite
! if irk > 1 we add on previous data
do ii = 1,irk-1
urk = urk + ButA(irk,ii)*k_stages(:,:,ii)
end do
call compute_flux(urk,flxcofs,2*m+1,nvar,hx)
call pde(ut,flxcofs,2*m+1,nvar,cflx)
k_stages(:,:,irk) = ut
end do
! At the end we add up the stages
do ii = 1,rkstages
uhermite = uhermite + ButB(ii)*k_stages(:,:,ii)
end do
end do
uh(:,:,i) = uhermite(0:m,1:nvar)
This section of code evolves the interpolant using a the
Runge-Kutta method defined by the Butcher arrays ``ButA`` and
``ButB`` that have been initialized by the ``call
set_butcher_tables``. The time stepping performs ``rksub`` local
timesteps with the local timestep ``dts = 0.5d0*dt/DBLE(rksub)``
to update ``uh``.
Local timestepping and CFL constraints
++++++++++++++++++++++++++++++++++++++
The purpose of this local timestepping is to maximize the
computation-to-communication ratio. For linear problems with
exact Taylor timestepping we have
already seen that we can run the solvers at the full physically
allowable CFL number. For non-linear problems with approximate
timestepping we can use large timesteps if the local evolution is
accurate enough.
For this example we can choose between three different RK solvers in
the module ``rkmodule.f90``, the classic fourth order Runge-Kutta, the
eight order DORPRI8 and the third order SSP(3,3) by Shu and Osher. To
illustrate the effect of local timestepping we consider the problem
defined by:
.. code-block:: fortran
INTEGER, PARAMETER :: m = 9 !spatial order 2m+1
INTEGER, PARAMETER :: nx = 100
DOUBLE PRECISION, PARAMETER :: xl = -1d0
DOUBLE PRECISION, PARAMETER :: xr = 1d0
DOUBLE PRECISION, PARAMETER :: ttot = 0.2d0
DOUBLE PRECISION, PARAMETER :: nu = 0.0001d0
for different ``CFL`` and ``rksub``. With DORPRI8 we can solve the
problem with one single local timestep at a ``CFL=0.65``, if we use
RK4 we must take ``CFL=0.30`` and if we use SSP(3,3) we must take ``CFL=0.20``.
However if we perform 6 local timesteps for SSP(3,3) or 4 local
timesteps with RK4 or 2 local timesteps with DORPRI8 we can use
``CFL=0.80``.
As this problem also has a viscous term we have choose the timestep
according to
.. code-block:: fortran
dt = min(CFL*hx/maxval(abs(u(0,1,:))),CFL_V*hx**2/(nu*m))
where ``CFL_V`` typically be chosen as 0.1.
The driver ``herm1d_semi_cons``
++++++++++++++++++++++++++++++++
The driver ``herm1d_semi_cons.f90`` illustrates an alternative
approach to discretize conservation laws. The alterative
discretization is motivated by the control volume statement of the
conservation law obtained by integrating over a space-time cell
.. math::
:nowrap:
\int_{0}^{\Delta t } \int_{x_i}^{x_{i+1}} u_t + f(u)_x \, dx dt = 0,
which can be simplified into:
.. math::
:nowrap:
\int_{x_i}^{x_{i+1}} u(\Delta t,x) - u(0,x) dx = - \int_{0}^{\Delta t } \left[ f(u(t,x_{i+1})) -f(u(t,x_{i})) \right] \,dt.
Now for a problem posed on a periodic domain (or with zero flux
boundary conditions) we have that the global mass is conserved
.. math::
:nowrap:
\sum_i \int_{x_i}^{x_{i+1}} u(\Delta t,x) dx = \sum_i \int_{x_i}^{x_{i+1}} u(0,x) dx.
For a numerical method to be able to have a similar conservation
property we must require that the flux from the left and right of an
interface coincide. Unless the flux is a linear function this is not
true in general for the above formulation. In the alternative
formulation in ``herm1d_semi_cons.f90`` we compute the fluxes at the
faces of each cell and interpolate them to the cell center. This will
guarantee that the fluxes are continuous at interfaces. The
interpolated fluxes are then differentiated to compute the first slope
in the Runge-Kutta method.
.. code-block:: fortran
DO i=0,nx
call compute_flux(u(:,:,i),flx(:,:,i),m,nvar,hx)
END DO
! Add the viscous flux
DO i=1,nx
! Semi-conservative
! Interpolate solution and fluxes to the cell center
call Taylor_to_Hermite1d(uhermite,u(:,:,i-1:i),Aintx,m,nvar)
call Taylor_to_Hermite1d(flxcofs,flx(:,:,i-1:i),Aintx,m,nvar)
do irk_sub = 1,rksub
if(irk_sub.gt.1) then
urk = uhermite
call compute_flux(urk,flxcofs,2*m+1,nvar,hx)
end if
! Get the first slope
call pde(ut,flxcofs,2*m+1,nvar,cflx)
k_stages(:,:,1) = ut
DO irk=2,rkstages ! do a local RK loop
urk = uhermite
! if irk > 1 we add on previous data
do ii = 1,irk-1
urk = urk + ButA(irk,ii)*k_stages(:,:,ii)
end do
call compute_flux(urk,flxcofs,2*m+1,nvar,hx)
call pde(ut,flxcofs,2*m+1,nvar,cflx)
k_stages(:,:,irk) = ut
end do
! At the end we add up the stages
do ii = 1,rkstages
uhermite = uhermite + ButB(ii)*k_stages(:,:,ii)
end do
end do
uh(:,:,i) = uhermite(0:m,1:nvar)
END DO
Although this formulation requires one additional interpolation it is
more robust than the non-conservative formulation. For the latter
multi dimensional examples with Euler's equations and compressible
Navier-Stokes equations, this is the formulation we will use.