25 November 2009

Who needs five digits of accuracy anyway

At the ANS Mathematics & Computational topical conference this spring, I heard a rumor (from a very well-placed source) about an ironic error in MCNP.

MCNP is considered by many in the labs and industry to be "the" nuclear simulation program of choice. It's a Monte Carlo transport program with continuous energy data and the ability to represent virtually exact geometry (with no truncation error for smooth surface, unlike in deterministic methods). It's had many decades of development, and it's been tested and tweaked considerably to match up with experimental and numerical benchmarks. Naturally, then, a simulated result from MCNP is often used as a reference answer when comparing new codes.

No simulation can exactly represent reality: there are always errors in data, simplifications in the physics, etc. This should always be kept in mind when doing computational anything. Case in point, I heard of a "run-off" benchmark between MCNP and another reactor code, SCALE. The result was that MCNP got the right answer and SCALE was wrong: the only problem was that their definition of the "right" answer was based on an MCNP-run calculation.

The error in MCNP—which you can verify yourself with the card entry print 98—is laughably simple yet indeterminately old. Avogadro's number, which converts between the number of atoms and the mass in a material, is programmed to be 6.022043447E+23 rather than the correct 6.02214179E+23. (In fact, it's written to 16 digits of "accuracy" even though only about ten are scientifically known to any precision.) The irony is that if they change Avogadro's number, their carefully benchmarked answers will come out "wrong."

It's not a huge problem, since I doubt most answers would be sensitive to such a small perturbation, but it serves as a warning to carefully consider the difference between "true" and "best approximation."

20 November 2009

Grammar fact of the day

I am apparently overly prone to hyphenating words. I knew that compound adjectives are supposed to be hyphenated. (Compound adjectives are modifiers composed of two words stuck together.) However, I wasn't sure if predicate adjectives were supposed to be hyphenated. I would certainly have a time-dependent simulation, but could I say that my simulation is time-dependent? As it turns out, I cannot grammatically say so. My problem is time dependent, as is my comprehension of the English language and its quaint usage rules.

SWIG and Traits classes

With all the new code I've been writing—namely, a large number of classes and template parameters and typedefs—I've been using a number of traits classes. This unfortunately posed a problem for SWIG (which provides a Python front-end to the C++ code), because it had trouble deducing the true object types through the maze (albeit short) of inheritance, typedefs, and template parameters. It had trouble despite the %include "Traits.hpp" command, which makes it aware of the traits class with template parameters.

My first solution was to wrap the troublesome traits-derived typenames with #if SWIG [explicitly substitute types] #else [use normal traits] #endif. It nullified a subset of the purpose of the traits class, but it did the job well enough.

The problem cropped up again today with a new class. Frustrated, I cast about for a solution, and tried explicitly instantiating the template in SWIG with %template(TraitsT) Traits<OneD, Gray>;. To my surprise, it worked! So apparently it's necessary to instantiate a traits class (even if it is a struct with no methods or data) in order to get SWIG to properly parse the typedefs in classes that use those traits.

EDIT: Looks like the manual actually has an entry about this. It recommends %template() traits<double,double>; to avoid generating unnecessary wrappers for it.