Burgers' Equation on a Moving Grid: Derivative Forcing

Solving Burgers' equation on a moving grid is a way to handle moving fronts or shocks by adapting the grid to the solution at each time-step. One simple monitor function is just the square of the solution derivative. Since we chose the forced diffusion equation as the governing equation for the grid, we just substitute the central difference approximation for the first derivative squared for the forcing function.

We also include two parameters (s1, s2) so we can adjust the strength of the forcing compared to the smoothing (second derivative term). That gives our semi-discrete equation for the rates as

Here's an example of the algorithm applied to an initial condition made-up of a couple exponential functions.

It shows that the delta-x decreases in areas of high gradients and increases in areas of zero derivative. The wave gets diffused pretty significantly because of the backward Euler time integration scheme.

Another illustrative example of this adaptive grid scheme is shown below. The initial condition is a couple wavelengths of cosine. This example shows that we probably want to force grid points into areas of high curvature (second derivative) as well as high first derivative. Notice the large delta-x around the peaks of each wave, zero derivative in this area, but large second derivative.

Central Schemes and Discontinuities

I've been posting lately on solving Burgers' equation on moving grids, but haven't said much about difficulties involved in properly approximating the non-linear term involving the spatial derivative. A naive approach to solving Burgers' equation with central differences leads to problems, termed "blow-up" by early practitioners doing climate modeling. Boyd discusses this problem as it relates to spectral methods.

A backward-time centered-space (BTCS) discretization of Burgers equation leads to a non-linear equation for each grid point that depends on it's nearest neighbors in the spatial coordinate and it's predecessor in the time coordinate. Backward time:

Centered space:

Resulting discrete Burgers' equation:

We can use a Newton method to solve the resulting system, and to do that we'll need the Jacobian:

To solve for the update in each iteration of the Newton's method, we'll generally need an expression for the action of the Jacobian on the update vector:

All of these symbolic manipulations of the governing equation can be done easily in Maxima (see here and here and here for examples).

Here's a graph of several time integration steps of the BTCS Burgers equation starting with part of a cosine:

This shows the characteristic behaviour of non-linear wave equations known as wave steepening, which causes the formation of shock waves (discontinuities). Unfortunately, if we take just a few more time-steps we get the beginnings of "blow-up":


Another interesting case is a small advecting step shown below.

In this case the size of the discontinuity is much smaller compared to the grid resolution, so the solution does not blow-up, but some weird things happen. Burgers' equation has one eigenvalue, or wave speed. All of the waves should travel in the plus x-direction, because u is positive everywhere in this case. That's not what we see though. We see the main step travel in the proper direction (after some smearing due to the first order time integration), but we see wiggly errors going in the opposite direction. Because the central differencing scheme does not respect the eigenvalues of the problem (the way an upwind scheme would), we introduced a new dispersive error wave that travels the wrong way.

Implicit-Time Burgers' Equation on a Moving Grid

In the last post on solving Burgers' equation on a moving grid we ended up with the semi-discrete equation for the rates:

where the spatial derivatives (for the solution and the grid) are approximated by simple central differences. We still haven't committed to a time-integration scheme. This is one of the nice things about the method of lines, we can play with different time integration schemes independent of our spatial discretization.

A good way to start out with implicit integration schemes is with the basic Backward Euler formula:

where the big X represents the unknown vector and the F(X) represents the right-hand side of our semi-discrete equation. The superscripts represent the time level. The reason this is an implicit time integration method is that the spatial terms are evaluated at the next time level rather than the present time level. This substitution is easily made in Maxima (just like the spatial derivatives):

/* substitute in a backward time difference: */
Discrete_Rate : ratsubst((u[i] - u[i,n-1])/dt, 'diff(u,t), Rate_Vector);
Discrete_Rate : ratsubst((x[i] - x[i,n-1])/dt, 'diff(x,t), Discrete_Rate);

Once the choice of time integration scheme is made, we have a fully-discrete problem that we can solve.

The Newton expression, G, in Maxima:

/* expression for Newton's method (zero finding): */
Newton_expression : fullratsimp(Discrete_Rate - Rate_expression);

Since G is non-linear, we'll need to use a Newton's method to solve for the zero (which will give us the solution and grid point locations at the next time-step). This requires the Jacobian:

/* Jacobian of the Newton expression: */
J : diff(Newton_expression, u[i-1]);
J : addcol(J, diff(Newton_expression, u[i]));
J : addcol(J, diff(Newton_expression, u[i+1]));
J : addcol(J, diff(Newton_expression, x[i-1]));
J : addcol(J, diff(Newton_expression, x[i]));
J : addcol(J, diff(Newton_expression, x[i+1]));

The update for each iteration of the Newton's method is given by:


Of course we generally don't invert the Jacobian directly to find the update vector, usually an iterative method such as Successive Over Relaxation or a Krylov subspace method is preferred. Most of these methods require a function which returns the application of the Jacobian on the update vector, so we go ahead and find that in Maxima as well:

/* Jacobian applied on the update vector: */
update_vector : transpose([du[i-1], du[i], du[i+1], dx[i-1], dx[i], dx[i+1]]);
Jdu : factor(J.update_vector);

The great thing about defining the numerical scheme in Maxima this way is that the f90() function will output these (increasingly complex) expressions to Fortran with the proper array indexing. If you've done this upfront part properly they are ready to go into loops as is.

Burgers' Equation on a Moving Grid

The familiar inviscid Burgers' equation on a non-moving, physical-space coordinate system:

Now we want to transform the equation, which is easily done in Maxima:

depends([u, x], [t]);
depends([xi], [x]);
depends([u], [xi]);

burgers : 'diff(u, t) + u * 'diff(u, x);

ev(burgers, nouns);

This gives us Burgers' equation transformed to computational space and a moving grid:

Now we've got some extra terms. The grid velocity:

And the grid transformation (to go from the computational space to the physical space):

The grid velocity is the troublesome term, because while we have an equation governing the time evolution of the unknown, we don't have a governing equation for the grid. This is where moving grids get fun, we get to pick the governing equation. How about a forced, unsteady, linear diffusion equation?

The method of lines is a popular and practical way to formulate PDEs for numerical solution. In Maxima, it is very straightforward to substitute in some finite difference expressions for the spatial derivatives, to get our semi-discrete equation for the time-derivatives.

/* substitute central differences for the spatial derivatives: */
burgers_discrete_space : ratsubst(u[i+1] - u[i-1], 'diff(u, xi), burgers_trans);
burgers_discrete_space : ratsubst(1/(x[i+1] - x[i-1]), 'diff(xi, x), burgers_discrete_space);

/* governing equation for the grid (forced diffusion equation): */
grid : 'diff(x,t) + 'diff(x, xi, 2) = f;
/* substitute central differences for the spatial derivative: */
grid_discrete_space : ratsubst(x[i+1] - 2*x[i] + x[i-1], 'diff(x,xi,2), grid);

The reason for transforming to computational coordinates is because the nodes are (usually) equally spaced in the computational domain. That way we can use the same difference formula for every point no matter what the spacing is in the physical space.

Now we can use the augcoefmatrix() function to get an operator and right-hand-side vector for our rate equation.

Rate_Operator : augcoefmatrix([burgers_discrete_space, grid_discrete_space], ['diff(u,t), 'diff(x,t)]);
Rate_RHS : -submatrix(Rate_Operator, 1, 2);
Rate_RHS : ratsubst(u[i], u, Rate_RHS);
Rate_Operator : submatrix(Rate_Operator, 3);
Rate_Operator_inv : invert(Rate_Operator);
Rate_Vector : transpose(['diff(u,t), 'diff(x,t)]);
Rate_expression : fullratsimp(Rate_Operator_inv.Rate_RHS);

Here's the resulting semi-discrete equations:

Because our choice of governing equation for the grid doesn't depend on the time-derivative of the solution, the operator is easily inverted:

and then applied, to give a vector of nonlinear expressions for our rates:

Note that in the absence of any forcing, the steady state grid is equally distributed grid points. The forcing function is what allows us to cluster grid points where higher resolution is needed. This forcing function should be based on some measure of the local error in the numerical solution. It is convenient just to use magnitude of the solution derivative since we've already calculated it, and areas with high gradients tend to need more points to resolve well.

Equal spacing in computational space is standard for using finite differences, but what about spectral methods where we usually choose nodes according to the Chebyshev polynomial roots? In this case it is convenient to choose our computational spacing to be the roots spacing rather than equal spacing. Here's the derivative of the coordinate with respect to computational coordinate using the DCT derivative routine.

Which shows that the roots nodes give a constant grid transform while the equally spaced nodes have a non-constant transform (opposite the normal situation of equally spaced computational coordinates).