National Strategy for Advancing Climate Modeling

NAP has a new report out on a national strategy for advancing climate modeling. I've just done a quick skim so far. Some tidbits below that touch on V&V and uncertainty.

Convergence for Falkner-Skan Solutions

About 6 months ago Dan Hughes sent me a link to an interesting paper on "chaotic" behavior in the trajectory of iterates in a numerical Falkner-Skan solution. It struck me that the novel results reported in that paper were an artifact of the numerical method, and had little to do with any "chaotic" physics that might be going on in boundary layers or other systems that might be well described by this equation. This is similar to the point I made in the Fun with Filip post: the choice of numerical method matters. Do not rush to judgment about problems until you have brought the most appropriate methods to bear.

There are some things about the paper that are not novel, and others that seem to be nonsense. It is well-known that there can be multiple solutions at given parameter values (non-uniqueness) for this equation, see White. There is the odd claim that "the flow starts to create shock waves in the medium [above the critical wedge angle], which is a representation of chaotic behavior in the flow field." Weak solutions (solutions with discontinuities/shocks) and chaotic dynamics are two different things. They use the fact that the method they choose does not converge when two solutions are possible as evidence of chaotic dynamics. Perhaps the iterates really do exhibit chaos, but this is purely an artifact of the method (i.e. there is no physical time in this problem, only the pseudo-time of the iterative scheme). By using a different approach you will get different "dynamics", and with proper choice of method, can get convergence (spectral even!) to any of the multiple solutions depending on what initial condition you give your iterative scheme. They introduce a parameter, \(\eta_{\infty}\), for the finite value of the independent variable at "infinity" (i.e. the domain is truncated). There is nothing wrong with this (actually it's a commonly used approach for this problem), but it is not a good idea to solve for this parameter as well as the shear at the wall in your Newton iteration. A more careful approach of mapping the boundary point "to infinity" as the grid resolution is increased (following one of Boyd's suggested mappings) removes the need to solve for this parameter, and gives spectral convergence for this problem even in the presence of non-uniqueness and the not uncommon vexation of a boundary condition defined at infinity (all of external aerodynamics has this helpful feature).

Now you have N problems

(10:40) It's not that regulators don't understand information technology, because it should be possible to be a non-expert and still make a good law. MPs and Congressmen and so on are elected to represent districts and people, not disciplines and issues.

That couple of sentences reminded me of a recent post on single-issue advocacy by Roger Pielke Jr. So, in that spirit, here's a fun word game: How applicable is Doctorow's criticism if you substitute "climate" for "copyright" below?

(22:20) But the reality is, copyright legislation gets as far as it does precisely because it's not taken seriously, which is why on one hand, Canada has had Parliament after Parliament introduce one stupid copyright bill after another, but on the other hand, Parliament after Parliament has failed to actually vote on the bill. [...] It's why the World Intellectual Property Organization is gulled time and again into enacting crazed, pig-ignorant copyright proposals because when the nations of the world send their U.N. missions to Geneva, they send water experts, not copyright experts; they send health experts, not copyright experts; they send agriculture experts, not copyright experts, because copyright is just not important to pretty much everyone!
Canada's Parliament didn't vote on its copyright bills because, of all the things that Canada needs to do, fixing copyright ranks well below health emergencies on first nations reservations, exploiting the oil patch in Alberta, interceding in sectarian resentments among French- and English-speakers, solving resources crises in the nation's fisheries, and thousand other issues! The triviality of copyright tells you that when other sectors of the economy start to evince concerns about the internet and the PC, that copyright will be revealed for a minor skirmish, and not a war. Why would other sectors nurse grudges against computers? Well, because the world we live in today is /made/ of computers. We don't have cars anymore, we have computers we ride in; we don't have airplanes anymore, we have flying Solaris boxes with a big bucketful of SCADA controllers [laughter]; a 3D printer is not a device, it's a peripheral, and it only works connected to a computer; a radio is no longer a crystal, it's a general-purpose computer with a fast ADC and a fast DAC and some software.

Full Transcript

Recurrence, Averaging and Predictability

Motivation and Background

Yet another installment in the Lorenz63 series. This time motivated by a commenter on Climate Etc. Tomas Milanovic claims that time averages are chaotic too in response to the oft repeated claim that the predictability limitations of nonlinear dynamical systems are not a problem in the case of climate prediction. Lorenz would seem to agree, “most climatic elements, and certainly climatic means, are not predictable in the first sense at infinite range, since a non-periodic series cannot be made periodic through averaging [1].” We’re not going to just take his word on it. We’ll see if we can demonstrate this with our toy model.

That’s the motivation, but before we get to toy model results a little background discussion is in order. In this previous entry I illustrated the different types of functionals that you might be interested in depending on whether you are doing weather prediction or climate prediction. I also made the remark, “A climate prediction is trying to provide a predictive distribution of a time-averaged atmospheric state which is (hopefully) independent of time far enough into the future.” It was pointed out to me that this is a testable hypothesis [2], and that the empirical evidence doesn’t seem to support the existence of time-averages (or other functionals) describing the Earth’s climate system that are independent of time [3]. In fact, the above assumption was critiqued by none other than Lorenz in 1968 [4]. In that paper he states,

Questions concerning the existence and uniqueness of long-term statistics fall into the realm of ergodic theory. [...] In the case of nonlinear equations, the uniqueness of long-term statistics is not assured. From the way in which the problem is formulated, the system of equations, expressed in deterministic form, together with a specified set of initial conditions, determines a time-dependent solution extending indefinitely into the future, and therefore determines a set of long-term statistics. The question remains as to whether such statistics are independent of the choice of initial conditions.

He goes on to define a system as transitive if the long-term statistics are independent of initial condition, and intransitive if there are “two or more sets of long-term statistics, each of which has a greater-than-zero probability of resulting from randomly chosen initial conditions.” Since the concept of climate change has no meaning for statistics over infinitely long intervals, he then defines a system as almost intransitive if the statistics at infinity are unique, but the statistics over finite intervals depend (perhaps even sensitively) on initial conditions. In the context of policy relevance we are generally interested in behavior over finite time-intervals.

In fact, from what I’ve been able to find, different large-scale spatial averages (or coherent structures, which you could track by suitable projections or filtering) of state for the climate system face similar limits to predictability as un-averaged states. The predictability just decays at a slower rate. So instead of predictive limitations for weather-like functionals on the order of a few weeks, the more climate-like functionals become unpredictable on slower time-scales. There’s no magic here, things don’t suddenly become predictable a couple decades or a century hence because you take an average. It’s just that averaging or filtering may change the rate that errors for that functional grow (because in spatio-temporal chaos different structures, or state vectors, will have different error growth rates and reach saturation at different times). Again Lorenz puts it well, “the theory which assures us of ultimate decay of atmospheric predictability says nothing about the rate of decay” [1]. Recent work shows that initialization matters for decadal prediction, and that the predictability of various functionals decay at different rates [5]. For instance, sea surface temperature anomalies are predictable at longer forecast horizons than surface temperatures over land. Hind-casts of large spatial averages on decadal time-scales have shown skill in the last two decades of the past century (though they had trouble beating a persistence forecast for much of the rest of the century) [6].

I’ve noticed in on-line discussions about climate science that some people think that the problem of establishing long term statistics for nonlinear systems is a solved one. That is not the case for the complex, nonlinear systems we are generally most interested in (there are results for our toy though [7, 8]). I think this snippet sums things up well,

Atmospheric and oceanic forcings are strongest at global equilibrium scales of 10⁷ m and seasons to millennia. Fluid mixing and dissipation occur at micorscales of 10^-3 m and 10^-3s, and cloud particulate transformations happen at 10^-6 m or smaller. Observed intrinsic variability is spectrally broad band across all intermediate scales. A full representation for all dynamical degrees of freedom in different quantities and scales is uncomputable even with optimistically foreseeable computer technology. No fundamentally reliable reduction of the size of the AOS [atmospheric oceanic simulation] dynamical system (i.e., a statistical mechanics analogous to the transition between molecular kinetics and fluid dynamics) is yet envisioned. [9]

Here McWilliams is making a point similar to that made by Lorenz in [4] about establishing a statistical mechanics for climate. This would be great if it happened, because that would mean that the problem of turbulence would be solved for us engineers too. Right now the best we have (engineers interested in turbulent flows and climate scientists too) is empirically adequate models that are calibrated to work well in specific corners of reality.

Lorenz was responsible for another useful concept concerning predictability, that is predictability of the first and second kind [1]. If you care about the time-accurate evolution of the order of states then you are interested in predictability of the first kind. If, however, you do not care about the order, but only the statistics, then you are concerned with predictability of the second kind. Unfortunately, Lorenz’s concepts of first and second kind predictability have been morphed in to a claim that first kind predictability is about solving initial value problem (IVP)s and second kind predictability is about solving boundary value problem (BVP)s. For example, “Predictability of the second kind focuses on the boundary value problem: how predictable changes in the boundary conditions that affect climate can provide predictive power [5].” This is unsound. If you read Lorenz closely, you’ll see that the important open question he was exploring about whether the climate is transitive, intransitive or almost intransitive has been assumed away by the spurious association of kinds of predictability with kinds of problems [1]. Lorenz never made this mistake, he was always clear that the difference in kinds of predictability depends on the functionals you are interested in, not whether it is appropriate to solve an IVP or a BVP (what reason could you have for expecting meaningful frequency statistics from a solution to a BVP?). Those considerations depend on the sort of system you have. In an intransitive or almost intransitive system even climate-like functionals depend on the initial conditions.

A good early paper on applying information theory concepts to climate predictability is by Leung and North [10], and there is a more recent review article that covers the basic concepts by DelSole and Tippett [11].

Recurrence Plots

Recurrence plots are useful for getting a quick qualitative feel for the type of response exhibited by a time-series [12, 13]. First we run a little initial condition (IC) ensemble with our toy model. The computer experiment we’ll run to explore this question will consist of perturbations to the initial conditions (I chose the size of the perturbation so the ensemble would blow-up around t = 12). Rather than sampling from a distribution for the members of the ensemble, I chose them according a stochastic collocation (this helps in getting the same results every time too).

---

(a)EnsembleTrajectories (b)EnsembleMean Figure 1: Initial Condition Ensemble

---

One thing that these two plots makes clear is that it doesn’t make much sense to compare individual trajectories with the ensemble mean. The mean is a parameter of a distribution describing a population of which the trajectories are members. While the trajectories are all orbits on the attractor, the mean is not.

---

(a)SingleTrajectory (b)EnsembleMean Figure 2: Chaotic Recurrence Plots

---

Comparing the chaotic recurrence plots with the plots below of a periodic series and a stochastic series illustrates the qualitative differences in appearance.

---

(a)PeriodicSeries (b)StochasticSeries Figure 3: Non-chaotic Recurrence Plots

---

Clearly, both the ensemble mean and the individual trajectory are chaotic series, sort of “between” periodic and stochastic in their appearance. Ensemble averaging doesn’t make our chaotic series non-chaotic, what about time averaging?

Predictability Decay

How does averaging affect the decay of predictability for the state of the Lorenz63 system, and can we measure this effect? We can track how the predictability of the future state decays given knowledge of the initial state by using the relative entropy. There are other choices for measures such as mutual information [10]. Since we’ve already got our ensemble though, we can just use entropy like we did before. Rather than just a simple moving average, I’ll be calculating an exponentially weighted one using an FFT-based approach, of course (there’s some edge effects we’d need to worry about if this were a serious analysis, but we’ll ignore that for now). The entropy for the ensemble is shown for three different smoothing levels in Figure 4 (the high entropy prior to t = 5 for the smoothed series is spurious because I didn’t pad the series and it’s calculated with the FFT).

---

Figure 4: Entropy of Exponentially Weighted Smoothed Series

---

While smoothing does lower the entropy of the ensemble (lower entropy for more smoothing / smaller λ), it still experiences the same sort of “blow-up” as the unsmoothed trajectory. This indicates problems for predictability even for our time-averaged functionals. Guess what? The recurrence plot indicates that our smoothed trajectory is still chaotic!

---

Figure 5: Smoothed Trajectory Recurrence Plot

---

This result shouldn't be too surprising, moving averages or smoothing (of whatever type you fancy) are linear operations. It would probably take a pretty clever nonlinear transformation to turn a chaotic series into a non-chaotic one (think about how the series in this case is generated in the first place). I wouldn't expect any combination of linear transformations to accomplish that.

Conclusions

I’ll begin the end with another great point from McWilliams (though I’ve not heard of sub-grid fluctuations referred to as “computational noise,” that term makes me think of round-off error) that should serve to temper our demands of predictive capability from climate models[9]:

Among their other roles, parametrizations regularize the solutions on the grid scale by limiting fine-scale variance (also known as computational noise). This practice makes the choices of discrete algorithms quite influential on the results, and it removes the simulation from the mathematically preferable realm of asymptotic convergence with resolution, in which the results are independent of resolution and all well conceived algorithms yield the same answer.

If I had read this earlier, I wouldn’t have spent so much time searching for something that doesn’t exist.

Regardless of my tortured learning process, what do the toy models tell us? Our ability to predict the future is fundamentally limited. Not really an earth-shattering discovery; it seems a whole lot like common sense. Does this have any implication for how we make decisions? I think it does. Our choices should be robust with respect to these inescapable limitations. In engineering we look for broad optimums that are insensitive to design or requirements uncertainties. The same sort of design thinking applies to strategic decision making or policy design. The fundamental truism for us to remember in trying to make good decisions under the uncertainty caused by practical and theoretical constraints is that limits on predictability do not imply impotence.

References

[1]   Lorenz, E. N., The Physical Basis of Climate and Climate Modeling, Vol. 16 of GARP publication series, chap. Climatic Predictability, World Meteorological Organization, 1975, pp. 132–136.

[2]   Pielke Sr, R. A., “your query,” September 2010, electronic mail to the author.

[3]   Rial, J. A., Pielke Sr, R. A., Beniston, M., Claussen, M., Canadell, J., Cox, P., Held, H., Noblet-Ducoudr, N. D., Prinn, R., Reynolds, J. F., and Salas, J. D., “Nonlinearities, Feedbacks And Critical Thresholds Within The EarthS Climate System,” Climatic Change, Vol. 65, No. 1-2, 2004, pp. 11–38.

[4]   Lorenz, E. N., “Climatic Determinism,” Meteorological Monographs, Vol. 8, No. 30, 1968.

[5]   Collins, M. and Allen, M. R., “Assessing The Relative Roles Of Initial And Boundary Conditions In Interannual To Decadal Climate Predictability,” Journal ofClimate, Vol. 15, No. 21, 2002, pp. 3104–3109.

[6]   Lee, T. C., Zwiers, F. W., Zhang, X., and Tsao, M., “Evidence of Decadal Climate Prediction Skill Resulting from Changes in Anthropogenic Forcing,” Journal of Climate, Vol. 19, 2006.

[7]   Tucker, W., The Lorenz Attractor Exists, Ph.D. thesis, Uppsala University, 1998.

[8]   Kehlet, B. and Logg, A., “Long-Time Computability of the Lorenz System,”http://lorenzsystem.net/.

[9]   McWilliams, J. C., “Irreducible Imprecision In Atmospheric And Oceanic Simulations,” Vol. 104 of National Academy of Sciences, National Academy of Sciences, pp. 8709 – 8713.

[10]   Leung, L.-Y. and North, G. R., “Information Theory and Climate Prediction,”Journal of Climate, Vol. 3, 1990, pp. 5–14.

[11]   DelSole, T. and Tippett, M. K., “Predictability: Recent insights from information theory,” Reviews of Geophysics, Vol. 45, 2007.

[12]   Eckmann, J.-P., Kamphorst, S. O., and Ruelle, D., “Recurrence Plots of Dynamical Systems,” EPL (Europhysics Letters), Vol. 4, No. 9, 1987, pp. 973.

[13]   Marwan, N., “A historical review of recurrence plots,” The European PhysicalJournal - Special Topics, Vol. 164, 2008, pp. 3–12, 10.1140/epjst/e2008-00829-1.

Easy Full-Factorial Ensemble Experiments

This is another installment of the Lorenz63 series. In this post I wanted to demonstrate an ’experiment’ on the Lorenz ’63 system. The basics of these sorts of computer experiments as they relate to climate modeling are covered at climateprediction.net.

The main factors underlying our experimental design are initial conditions, parameters and forcing. Since we followed Roache’s advice about making our code Verifiable with the method of manufactured solutions (MMS) the ability to specify forcings is already included. Python has some nice built-in functions in the itertools module that make doing a computational experiment really easy. First all of the arguments for the time integrator are defined as lists (the ones which aren’t part of the experiment only have one element, so they remain constant):

# factors in our experimental design (all of the inputs to the system 
# integrator, some only have one level): 
nt = 2000 
savevery = [10] 
nsaves = nt / savevery[0] + 1 
T = 2.0 
# the first element of the array is the initial condition, the rest 
# are over-written with the solution: 
x = [1.05*sp.ones(nsaves,dtype=float,order=’F’), 
     0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
y = [-1.05*sp.ones(nsaves,dtype=float,order=’F’), 
      -0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
z = [-10.5*sp.ones(nsaves,dtype=float,order=’F’), 
         -9.5*sp.ones(nsaves,dtype=float,order=’F’)] 
t = [sp.zeros(nsaves, dtype=float, order=’F’)] 
Ra = [27.0, 28.0, 29.0] 
Pr = [9.0, 10.0, 11.0] 
b = [8.0/3.0 - 0.1, 8.0/3.0 + 0.1] 
h = [T / float(nt)] 
tol = [1e-14] 
maxits = [20] 
forcing = [0, 1] 
bx = [0.0] 
by = [-1.0] 
bz = [-10.0]

Then a call to itertools.product() gives us the Cartesian Product of all of those lists, which is a full-factorial set of the inputs (2 ⋅ 2 ⋅ 2 ⋅ 3 ⋅ 3 ⋅ 2 ⋅ 2 = 288 treatments in this example).

# this gives us an iterator of every combination of our factor levels: 
full_factorial = itertools.product( 
        x, y, z, t, Ra, Pr, b, h, tol, maxits, forcing, bx, by, bz, savevery)

Once we have the iterator defined, we can calculate all of the trajectories associated with the input treatments in the experiment:

# iterate over all of the ’treatments’ in our design, the iterator 
# yields the arguments directly: 
for args in full_factorial: 
    L63.time_loop(*args) 
    xresult.append(args[0].copy()) 
    yresult.append(args[1].copy()) 
    zresult.append(args[2].copy()) 
    tresult.append(args[3].copy())

The resulting trajectories are shown in Figure 1. It’s pretty clear that realistic uncertainties in parameters and initial conditions lead quickly to pretty diffuse predictive distributions. Part of the challenge of attribution studies is determining the effect of each change in parameter or forcing. In realistic models, the number of ’factors’ to vary is quite large.

---

(a) X component (b) Y component (c) Z component Figure 1: Trajectories for the Full-Factorial Experiment

---

With such a low-dimension system, and few parameters it is easy enough to run the whole full-factorial, smarter experimental designs are required for experiments based on large simulation runs. A good reference on this is Design and Analysis of Computer Experiments.

Chaotic Time-Convergence and MMS

This previous post covered time-convergence of a single trajectory, and the growth of the spread in an ensemble. The next post covered some special considerations that chaos imposes on us when we try to apply the method of manufactured solutions (MMS). I’ll try to tie up some loose ends in this (perhaps final) installment of the Lorenz63 series. Two questions that I think I haven’t adequately addressed yet:

• What about the time-step convergence behavior of an ensemble?
• So can you successfully apply the MMS to verify implementations of numerical approximations to this system or not? (The literature indicates that you can, but I didn’t actually demonstrate the correctness of my implementation.)

To address the first question I wrote a little Python script that runs a range of trajectories with various time-steps and also perturbations to the initial conditions:

nt = [192001, 384001, 768001] # number of time-steps 
step = [1, 2, 4] 
T = 60.0 # integration period 
# parameters for the Lorentz 63 system to get chaotic response: 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
# knobs for the Newton’s method: 
tol = 1e-14 
maxits = 20 
# number of different IC’s in our ensemble: 
nic = 5 
# number of pairs of differences between ensemble members: 
npairs = comb(nic, 2, exact=1) 
h = [] 
x = [] 
xdiffs = [] 
y = [] 
ydiffs = [] 
z = [] 
zdiffs = [] 
t = [] 
for n in nt: 
    h.append(T / float(n)) 
    x.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    xdiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    y.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    ydiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    z.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    zdiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    t.append(sp.linspace(0.0, T, n)) 
 
# some perturbations: 
eps = 1e-16 * sp.array([0, 1, -1, 2, -2]) 
for i in xrange(len(nt)): 
    x[i][0,:] = 1.0 + eps 
    y[i][0,:] = -1.0 + eps 
    z[i][0,:] = 10.0 + eps 
 
# calculate all of the trajectories: 
for i in xrange(len(nt)): # time-step variation loop 
    for j in xrange(nic): # perturbations loop 
        L63.time_loop( 
            x[i][:,j],y[i][:,j],z[i][:,j],t[i],Ra,Pr,b,h[i],tol,maxits)

Here the result for the x-component of the solution:

---

Figure 1: Convergence and Chaotic spread behavior of Lorenz ensemble

---

Figure 1 shows that we have a time-step converged solution out to t ≈ 30. After which, the trajectory between the different time-steps is significantly different, but the ensemble spread is still negligible. At t ≈ 40 the ensemble begins to spread noticeably for all three time-step sizes, and then shortly thereafter we get an abrupt divergence of the trajectories for all of the ensemble members.

The interesting question is whether the divergence is inevitable (that’s my impression from reading the ’chaos’ literature, let me know if you think/know something different), or if getting more and more accurate solutions will delay that abrupt divergence (the visible divergence seems to happen latter in the most resolved ensemble). Unfortunately, I’m already flirting with swap-thrash on my little laptop with these ensembles, so I need to do a little re-factoring before I’ll be able to demonstrate that one definitively (so this probably isn’t the final installment). What’s really required here is a time-step converged solution of the region t > 40 where the ensemble members diverge. These solutions aren’t there yet (smaller time-step or higher-order method is needed).

The second part of the post will deal with demonstrating an MMS-based verification of my (purportedly third-order accurate) implementation. As discussed in this post, since the system is chaotic, choosing a manufactured solution becomes a bit more difficult. Here’s a Python script which calculates trajectories with MMS forcing and calculates error norms:

nt = [501, 1001, 2001, 4001, 8001] # number of time-steps 
step = [1, 2, 4, 8, 16] 
T = 2.0 * sp.pi # integration period 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
tol = 1e-14 
maxits = 20 
mms_switch = 1 
h = [] 
x = [] 
y = [] 
z = [] 
t = [] 
ms = [] 
for n in nt: 
    h.append(T / float(n)) # time-step size 
    x.append(sp.zeros(n, dtype=float, order=’F’)) 
    y.append(sp.zeros(n, dtype=float, order=’F’)) 
    z.append(sp.zeros(n, dtype=float, order=’F’)) 
    t.append(sp.linspace(0.0, T, n)) 
    ms.append(sp.zeros((n,3), dtype=float, order=’F’)) 
 
bx = 0.0 
by = -1.0 
bz = 10.0 
# manufactured solutions: 
for i in xrange(len(nt)): 
    ms[i][:,0] = sp.cos(t[i]) + bx 
    ms[i][:,1] = sp.sin(t[i]) + by 
    ms[i][:,2] = sp.sin(t[i]) + bz 
 
# trajectories: 
for i in xrange(len(nt)): 
    x[i][0] = bx + 1.0 
    y[i][0] = by 
    z[i][0] = bz 
    L63.time_loop( 
        x[i],y[i],z[i],t[i],Ra,Pr,b,h[i],tol,maxits,mms_switch,bx,by,bz) 
 
# error L2 norms: 
err_x = sp.zeros(len(nt), dtype=float) 
err_y = sp.zeros(len(nt), dtype=float) 
err_z = sp.zeros(len(nt), dtype=float) 
for i in xrange(len(nt)): 
    err_x[i] = sp.sqrt(sp.sum((x[i] - ms[i][:,0])**2)/nt[i]) 
    err_y[i] = sp.sqrt(sp.sum((y[i] - ms[i][:,1])**2)/nt[i]) 
    err_z[i] = sp.sqrt(sp.sum((z[i] - ms[i][:,2])**2)/nt[i]) 
 
# find the observed order of convergence: 
px = sp.polyfit(sp.log(h), sp.log(err_x), 1) 
py = sp.polyfit(sp.log(h), sp.log(err_y), 1) 
pz = sp.polyfit(sp.log(h), sp.log(err_z), 1)

Figure 2 shows the behavior of the solution with MMS forcing when the manufactured solution is not in the convergent region. The solution really is not converging.

---

(a) X-component with MMS forcing (b) Error convergence Figure 2: Convergence behavior with MMS forcing, bx = 0.0, by = -1.0, bz = 10.0

---

Contrast this with the behavior in Figure 3, which actually shows at least first order convergence (and the solution looks like the one we chose).

---

(a) X-component with MMS forcing (b) Error convergence Figure 3: Convergence behavior with MMS forcing, bx = 10.0, by = -10.0, bz = 40.0

---

The fact that the slope of the line in Figure 3 is 1 rather than 3 means I still have a lurking ordered-error in my implementation (likely in the error calculations themselves), or I’m not yet in the asymptotic range, though this seems unlikely given the small time-step size and very consistent slope. The good news illustrated by Figure 3 is that we can apply MMS to even chaotic systems. Now where’s that bug…

Chaotic Method of Manufactured Solutions: Lorenz 63

I quickly brushed over the method of manufactured solutions (MMS) aspects in the previous post on Lorenz 63, you might wonder why. Well, it turns out that applying the MMS to chaotic systems includes an additional consideration that is not normally of concern. Here is the ’standard’ sort of checklist for choosing manufactured solutions (taken from Verification of Computer Codes in Computational Science and Engineering, emphasis mine):

1. Manufactured solutions should be sufficiently smooth on the problem domain so that the theoretical order-of-accuracy can be matched by the observed order-of-accuracy obtained from the test. A manufactured solution that is not sufficiently smooth may decrease the observed order-of-accuracy (however, see # 7 in this list).
2. The solution should be general enough that it exercises every term in the governing equation. For example, one should not choose temperature T in the unsteady heat equation to be independent of time. If the governing equation contains spatial cross-derivatives, make sure the manufactured solution has a non-zero cross derivative.
3. The solution should have a sufficient number of nontrivial derivatives. For example, if the code that solves the heat conduction equation is second order in space, picking T in the heat equation to be a linear function of time and space will not provide a sufficient test because the discretization error would be zero (to within round-off) even on coarse grids.
4. Solution derivatives should be bounded by a small constant. This ensures that the solution is not a strongly varying function of space or time or both. If this guideline is not met, then one may not be able to demonstrate the required asymptotic order-of-accuracy using practical grid sizes. Usually, the free constants that are part of the manufactured solution can be selected to meet this guideline.
5. The manufactured solution should not prevent the code from running successfully to completion during testing. Robustness issues are not a part of code order verification. For example, if the code (explicitly or implicitly) assumes the solution is positive, make sure the manufactured solution is positive; or if the heat conduction code expects time units of seconds, do not give it a solution whose units are nanoseconds.
6. Manufactured solutions should be composed of simple analytical functions like polynomials, trigonometric, or exponential functions so that the exact solution can be conveniently and accurately computed. An exact solution composed of infinite series or integrals of functions with singularities is not convenient.
7. The solution should be constructed in such a manner that the differential operators in the governing equations make sense. For example, in the heat conduction equation, the flux is required to be differentiable. Therefore, if one desires to test the code for the case of discontinuous thermal conductivity, then the manufactured solution for temperature must be constructed in such a way that the flux is differentiable. The resultant manufactured solution for temperature is nondifferentiable.
8. Avoid manufactured solutions that grow exponentially in time to avoid confusion with numerical instability.

As emphasized in # 4 we often have to use the free constants in our manufactured solution to bound the derivatives in the domain. For chaotic systems we’ll need to use not only the free constants, but the initial conditions as well, to achieve another end: ensure that the chosen solution is in the convergent region at all times. This is a sufficient but not necessary condition for the source term to be in the basin of entrainment. When the source term is in the basin of entrainment then we can expect to be able to verify the ordered convergence of our implementation as we would in the non-chaotic case. Otherwises we can expect the truncation error to eventually lead to an arbitrarily large departure from the manufactured solution. This terminology comes from the nonlinear dynamics and control literature [1] [3]. A worked example for the Lorenz ’63 system showing the boundary of the convergent region is in Chapter 6 of [2].

The convergent region is defined by the region of phase space for which the eigenvalues of the system Jacobian has negative real parts. As in any other application of the MMS, we start by choosing some convenient solution consistent with the desiderata enumerated above. Trigonometric functions are convenient in this case:

(1)

(2)

(3)

Now we need the Jacobian of the system (and then its eigenvalues), since our governing equations are definied in Maxima, that is easy enough:

/* Need the Jacobian of the system to find the ’convergent region’ so 
   we can select a useful IC/amplitude for our manufactured solutions */ 
sys_eqns : [eqn_1, eqn_2, eqn_3] $ 
J_cr : zeromatrix(3,3) $ 
for j : 1 thru 3 do 
( for i : 1 thru 3 do 
  ( 
    J_cr[i,j] : diff(rhs(sys_eqns[i]), unknowns[j]) 
    ) 
  ) $ 
[J_cr_eigvals, J_cr_eigmult] : eigenvalues(J_cr) $ 
J_cr_eigvals : fullratsimp(J_cr_eigvals) $ $

As before these expressions are output in Fortran90 format for compilation. Here are the eigenvalues in the complex plane for b_x = 1.0, b_y = -1.0, b_z = 10.0.

---

Figure 1: Jacobian eigenvalues for manufactured solution, bx = 1.0, by = -1.0, bz = 10.0

---

We can immediately see the problem from Figure 1: our manufactured solution is in the region of phase space where one of the Jacobian’s eigenvalues has a positive real part. We can adjust the initial conditions of our solution to get all of our eigenvalues on the left half-plane.

---

Figure 2: Jacobian eigenvalues for manufactured solution, bx = 10.0, by = -10.0, bz = 40.0

---

Being limited to only certain portions of the phase space for conducting verification is fine, because a single well-designed manufactured solution verifies the correctness of the entire implementation (as long as all of the terms are activated). This is the same reason that it is often emphasized in the V&V literature that there is no need for the manufactured solution to be physically meaningful.

References

[1]   Jackson, E.A., “Controls of dynamic flows with attractors,” Physics Review A, Vol. 44, Is. 8, 1991.

[2]   Wu, W., “Analytical and Numerical Methods Applied to Nonlinear Vessel Dynamics and Code Verification for Chaotic Systems,” Disertation, Virginia Polytechnic Institute and State University, Dec, 2009.

[3]   W. Wu, L.S. McCue, and C.J. Roy. “The method of manufactured solutions applied to chaotic systems,” Nonlinear Dynamics, 2009, under review.