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.

Micro Air Vehicle Aeroelastic Analysis

VDM-Verlag publishes specialist monographs, and they decided that my old masters thesis fit the bill. I imagine that it's probably because of the very rapid growth in the use of unmanned aerial vehicles (UAVs) for surgical strikes, close air support and long-duration reconnaissance missions in Iraq and Afgahnistan. The full title is Aeroelastic Analysis of Micro Air Vehicles: Coupled Computational and Experimental Techniques.

My book is about a coupled aeroelastic method that uses mode-shapes for the flexible wing Micro Air Vehicle (MAV) that were determined with laser-vibrometry. These mode-shapes are loaded and integrated in time with a computational fluid dynamics (CFD) code. It's a unique way to model the behaviour of these flexible wing vehicles. The reason for including a flexible wing is for portability (roll it up and stick it in a tube), and for damping gust response. This second item can provide for a more stable sensor/shooter platform.

The "about this book blurb" from the back cover:

A computational aeroelastic analysis of a micro air vehicle (MAV) is conducted. This MAV has a 24-inch wing span, and is designed for local area reconnaissance. Wind tunnel data for the MAV with a rigid carbon fiber wing and a flexible carbon fiber-ribbed nylon wing are compared to CFD results incorporating static and dynamic deformations. We use laser vibrometry to determine the mode shapes of the flexible wing. From these shapes, CFD grid deformations are calculated as part of a closely-coupled aeroelastic solution method. The accuracy of MAV performance predictions using CFD with and without aeroelastic modeling is evaluated against previous wind-tunnel experiments. The performance benefits of the flexible wing, and the applicability and limitations of the model are evaluated in the present research effort. Some suggestions are made as to improvements that can be made to increase the range of applicability and the accuracy of the model.

My experience with VDM publishing was a positive one. They were quick, efficient and professional. The proof copy of the book turned out well (just make sure all your figures look decent in black and white).

Further Reading

• A related AIAA paper I presented in 2005 in Toronto that's focused more on the computational aspects.
  
  
• more UAVs in the news about Afghanistan.
  
  
• A slightly more critical view of UAVs
  
  
• A pretty good New York Times write-up on UAVs
  

Treating Wounds from PowerPoint Bullets

Or "how to retain your critical thinking skills while being PowerPointed to death".

I recently wrote about PowerPoint Engineering from the perspective of improving technical communications towards decision makers. What about when you are on the receiving end of a volley of PP bullets?

Tufte has a nice critical write-up with "methods for how not to get fooled" while enduring a presentation. The gist (in "bullet" form of course):

• Be skeptical of analytical claims of "conservatism", these are often rhetorical tricks which hide a lack of quantitative rigor: "we're so uncertain we can't put a number on it, so we were conservative"
  
• Close reading pays, beware of contradictory sub-sub-bullets tucked in under vague "big" bullets or headings
  
• Be skeptical of claims of "significance" without supporting quantitative statistical analysis, another rhetorical slight of hand
  
• Be aware of the scope of an analysis, true statements about a particular case are easily misunderstood to apply generally
  
• Disambiguate pronouns, often the importance of the subject (life or death) is hidden by innocuous looking pronouns
  
• Be aware of your organization's biases as well as your own (echos of Aristotle's Nicomachean Ethics)
  

Of course, in the end it boils down to being an engaged, critical thinker. PowerPoint has been shown to be ill-suited for technical communication, so the information you need is probably not on the slides! To quote from the Tao Te Ching:

We put thirty spokes together and call it a wheel;
But it is on the space where there is nothing that the usefulness of the wheel depends.
We turn clay to make a vessel;
But it is on the space where there is nothing that the usefulness of the vessel depends.
We pierce doors and windows to make a house;
And it is on these spaces where there is nothing that the usefulness of the house depends.
Therefore just as we take advantage of what is, we should recognize the usefulness of what is not.

Understanding their omissions is as important as understanding their presentations. In modern, complex systems we often suffer from information and data overload, so there's likely to be lots omitted. Hopefully not the critical few pieces of information, but failure analysis after failure analysis shows it often is.

As analysis becomes more causal, multivariate, comparative, evidence-based, and resolution-intense, the more damaging the bullet list becomes. -- E. Tufte, on bullets

A better list of how to handle the PP barrage than my paraphrase above (emphasis mine):

One's inability to produce the devastating comeback live during the course of discussion-- l'esprit de l'escalier--is why it is helpful, in deciding serious matters (1) to have the material under discussion distributed in advance, (2) to rehearse the possible exchange in advance, (3) to take a time-out, leave the meeting, escape the groupthink and bullying, go for a long walk down the hall (or up the stairs in the French version of after-the-fact-wit), and ask yourself "What would Richard Feynman do?"-- E. Tufte, on analytical reasoning

Beware technocrats bearing bullet points.

NASA keeps coming up in these sorts of discussions, not because they are particularly bad, but because they are a large, technical organization that is public. They don't hide behind the veil of secrecy (proprietary or classified). A recent discussion on Slashdot highlights some interesting stuff. One being that someone who by all rights should be technically competent can behave like a foolish bureaucrat when addled by the wrong sort of system.

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).

PowerPoint Engineering Considered Harmful

I have always been inspired by Edward Tufte's work. As an engineer in the aerospace industry, successfully working in large organizations (read bureaucracies) is a necessary skill. Unfortunately, bureaucracy has its own particular pathology. The trick is to not catch the diseases yourself while at the same time thriving in close proximity. In that vein, I have lately been considering the Cognitive Style of PowerPoint, which is a popular tool in the industry.

Everybody likes to beat up on NASA, and there's plenty to find fault with (as there is in any bureaucracy). Honest past criticism from a much admired physicist:

It appears that there are enormous differences of opinion as to the probability of a failure with loss of vehicle and of human life. The estimates range from roughly 1 in 100 to 1 in 100,000. The higher figures come from the working engineers, and the very low figures from management. What are the causes and consequences of this lack of agreement?
...
Let us make recommendations to ensure that NASA officials deal in a world of reality in understanding technological weaknesses and imperfections well enough to be actively trying to eliminate them. They must live in reality in comparing the costs and utility of the Shuttle to other methods of entering space. And they must be realistic in making contracts, in estimating costs, and the difficulty of the projects. Only realistic flight schedules should be proposed, schedules that have a reasonable chance of being met. If in this way the government would not support them, then so be it. NASA owes it to the citizens from whom it asks support to be frank, honest, and informative, so that these citizens can make the wisest decisions for the use of their limited resources.

For a successful technology, reality must take precedence over public relations, for nature cannot be fooled.
-- Rogers Commission, Appendix F, by R. P. Feynman

The same "debt of honesty" that Feynman claims NASA owes the public can be reasonably attributed to any organization towards their stakeholders. Pretending the problems experienced by this particular organization are isolated would be foolish. NASA does provide an accessible case study of an organization being firmly infected by the pitch culture though. Others might put it as the triumph of the personality-ethic over the character-ethic, style over substance.

From a more recent report precipitated by disaster:

The Mission Management Team Chairʼs position in the hierarchy governed what information she would or would not receive. Information was lost as it traveled up the hierarchy. A demoralized Debris Assessment Team did not include a slide about the need for better imagery in their presentation to the Mission Evaluation Room. Their presentation included the Crater analysis, which they reported as incomplete and uncertain. However, the Mission Evaluation Room manager perceived the Boeing analysis as rigorous and quantitative. The choice of headings, arrangement of information, and size of bullets on the key chart served to highlight what manage- ment already believed. The uncertainties and assumptions that signaled danger dropped out of the information chain when the Mission Evaluation Room manager condensed the Debris Assessment Teamʼs formal presentation to an informal verbal brief at the Mission Management Team meeting.
-- CAIB Report, Vol. 1, pp 201

Any technical organization will have trouble clearly and accurately communicating technical information to decision makers, who most likely are not experts in the particular domain being addressed, and often are not even technically competent in any speciality. Even when they are experts, using an inappropriate tool never helps.

From Thomas Ricks' book Fiasco:

[Army Lt. General David] McKiernan had another, smaller but nagging issue: He couldn't get Franks to issue clear orders that stated explicitly what he wanted done, how he wanted to do it, and why. Rather, Franks passed along PowerPoint briefing slides that he had shown to Rumsfeld: "It's quite frustrating the way this works, but the way we do things nowadays is combatant commanders brief their products in PowerPoint up in Washington to OSD and Secretary of Defense...In lieu of an order, or a frag [fragmentary order], or plan, you get a bunch of PowerPoint slides...[T]hat is frustrating, because nobody wants to plan against PowerPoint slides." That reliance on slides rather than formal written orders seemed to some military professionals to capture the essence of Rumsfeld's amateurish approach to war planning. "Here may be the clearest manifestation of OSD's contempt for the accumulated wisdom of the military profession and of the assumption among forward thinkers that technology--above all information technology--has rendered obsolete the conventions traditionall governing the preparation and conduct of war," commented retired Army Col. Andrew Bacevich, a former commander of an armored cavalry regiment. "To imagine that PowerPoint slides can substitute for such means is really the height of recklessness." It was like telling an automobile mechanic to use a manufacturer's glossy sales brochure to figure out how to repair an engine.

This raises some of the same issues discussed in the report by members of the NASA Return to Flight Task Force in the contribution at the top of this thread.
-- Edward Tufte, August 11, 2006

Writing a good technical report or white-paper takes thought and hard work. Clearly communicating intent with effective mission orders takes initiative and competence. Throwing bullets on a slide is the path of less resistance (I say this as someone who has regrettably succumb to this temptation on more than one occasion). It is also the refuge of the bureaucrat "of an entirely different character" planning to profit from ambiguity and disingenuous politicking.

Knowing a thing is bad is nice in a theoretical way, but what about actually doing well because of that knowledge. How to effectively communicate technical results to hurried managers? I think the Tufte-inspired handout format is a nice start. Doing a good job of tightly incorporating relevant graphics into succinct narrative is probably the best an engineer can do. Improving the report to the point that it isn't a long-winded, dry, snooze-inducer means that management may actually read it, and their decisions may benefit from it. Of course, the hand-out can reference more in-depth technical documents for the competent, curious and critical reader.

Chebyshev Spectral Method Shock Filtering

Now that I've got a nice handy spectral derivative routine, time to test it out on some practical problems (ie non-smooth).

A useful scalar test function is a square wave with a little Gaussian blip added to it:

This is good test function because it contains discontinuities as well as small smooth features that we'd like to be able to resolve.

There are tons of options to filter out high frequency "ringing" around discontinuities (Gibbs phenomenon). One nice option is the Butterworth filter, the gain for which is:

Where the cut-off frequency and the order of the filter are free parameters you can "fiddle" with. For this set of derivatives I'm going to choose a cut-off frequency of half the highest frequency in my grid and see what the effect of the order parameter is. Here's a plot of the gain for a couple different order filters:

Notice they all have the same gain at the cut-off frequency, they just pass through it differently.

The un-smoothed spectral derivative compared to the derivative after different order Butterworth filters are applied to the Chebyshev coefficients before the inverse transform are shown below.

The presence of the discontinuities causes "ringing" throughout the entire grid, but the filters do a pretty good job of containing it to the area immediately around the shocks. Here's a zoomed in view of the smooth feature area:

Without filtering the derivatives of the smooth feature are totally swamped by ringing, but you actually get a nice result out after some filtering.