Silo now available on MacPorts

Scientific computing developers rejoice. If you want your code to output files readable by LLNL's VisIT, there's a good chance you use the Silo library to do it. Well, if you use MacPorts, you can now download and compile silo without hassle on a Mac. I submitted a port file for Silo this morning, and the great volunteers at MacPorts had it cleaned up and in their repository by the end of the day. Now all you need to do is sudo port selfupdate && sudo port install silo, and you're good to go.

Research code is now open sourced

Well, after talking with a group of nuclear engineers who are interested in open sourcing their codes, and after being persuaded by my friend Mike that it was justifiable to return some of my tools to the taxpayers who've sponsored my education, I have decided to open up my research code. The project is called PyTRT, and it's available at GitHub under a Simplified BSD license. Documentation for PyTRT is available, but it is probably overwhelming. We'll see if this ever amounts to anything.

Time-dependent discrete ordinates transport

The discrete ordinates approximation for transport lumps "particles" into certain directions. In a 1-D time-dependent transport problem, the effect is the same as having several waves propagating at different speeds from the source region.

Regression

A study titled The largest Last Supper: depictions of food portions and plate size increased over the millennium has attempted to extrapolate nominal meal portion sizes from paintings of The Last Supper created over an 800-year span. Now, tell me if this is not the worst regression analysis you've ever seen:

Really, a parabolic fit? One that extends 300 years beyond the range of the data? I'm no expert in regression analysis, but even an r-value of 0.5 does not seem worthy of the strong conclusions which the article draws.

Aside from the poor correlation coefficient, and assuming the very method of linking plate-to-head ratio to portion size is valid, it looks like their extrapolation hinges entirely on one outlier (the most recent painting) and a century of plenty (in the 1500's, where there were not only small portions but large portions). This is a crock.

It's science, ok

Today I had a little fun with my code. As I've mentioned too often, it's written with a Python front end. That means I can very easily interact with the hundreds of Python modules out there. Well, this being a Friday afternoon, I decided to tie it into an image processing library, which will let me turn an image into a series of values, which I can then define as a field on my mesh (e.g., the total cross section). So I did that with a test image, and ran it through my transport code.

The code to load an image as grayscale, create the mesh, and set the cross sections:

    im = Image.open(fname).convert("L").transpose(Image.FLIP_TOP_BOTTOM )
    values = ( (255.0 - x) / 255.0 for x in im.getdata() )
    args = list(im.size) + list( 10.0 * float(i) / max(im.size)
                                  for i in im.size )
    mesh = meshlib.Mesh(*args)
    sigma_t = meshlib.CellFieldFloat(mesh, meshlib.VectorFloat(values) )

I also did this with a picture of my advisor.

MCNP criticality search Python script

One of the occasional tasks of a nuclear engineering student is to perform a criticality search using some parameter. That is, by varying (for example) the radius r of a uranium sphere, at some point the neutron production (from fission) and loss (from absorption and leakage) terms exactly balance, and the reactor is critical. The issue is complicated a little when using a Monte Carlo simulation, since the k-eigenvalue (unity when critical) is only reported with a certain probability or tolerance; running more particles will decrease the uncertainty.

Undergraduates, being unthinking animals, usually try this arduous task by hand, laboriously changing the parameter in an input file, running it, and repeating. As a graduate student, such tediousness is beneath me, as I would rather spend an equal amount of time writing a script to do the task for me. (Experience dictates that most of the time, some error or revision will require re-doing the task at least once, which means that a script will often save time by a factor of two or more. Besides, scripts are reusable and an educational experience.)

These Python scripts (gzipped archive) will do just that. It's a fairly simple binary search (no Newton–Krylov iterations) partly because of the complicating factor of the uncertainty in each answer. When criticality is reported to within the experimental standard deviation, the script increases the number of particles until such a parameter is determined that results in criticality within the user's specified tolerance. The script will run MCNP, process the output to determine k-effective, and iterate to your heart's content. If you have matplotlib installed, it will also produce a plot of how k varies with the parameter:

The MCNP search script is by no means perfect, and it assumes you have mcnp5 in your path. You have to know some Python to create a subclass of McnpSolver for your particular problem.

Included are scripts for finding the critical parameters of HEU cubes, cylinders, hemispheres, etc. There's also a script for generating a spherical shell volume tally to calculate the radial flux of the converged critical sphere.

Future updates of the MCNP criticality script will be posted on my projects page, although you should not expect there to be any. If you use this script, please cite its use if you report or present the results.

Insolation pt. 2

Incidentally, the plot I posted yesterday was only of the intensity averaged over the daylight hours. Multiplying by the hours of daylight (generated with the same program) should give the total amount of solar energy deposition per unit area in a day.

UPDATE This actually seems to give the wrong values, according to here and here. Maybe the NASA calculations are weighted differently than they say, or I'm misunderstanding what they mean. Either way, ignore this plot and use the values shown by those links.

Insolation

Today's topic is not "insulation," as in "my apartment doesn't have enough insulation to keep me warm in this frigid weather." Rather, insolation is how much sunlight hits the earth's surface.
Ever since moving up north and noting how low the sun is in the sky and how short are the days during winter, I've wondered about how much actual energy gets dumped onto a unit area on the surface of the earth.

After googling all sorts of permutations on "day integrated solar flux" and similar terms, I stumbled upon the more accurate term of insolation (which, after the fact, I remembered having heard through science bowl). Searching around some more resulted in this factual page, and further searching led me to a NASA page with all sorts of useful information. Based on values from one of NASA's programs, which takes into account most of the important stuff, here is a plot of the solar intensity as a function of the days of the year:

It's no wonder it feels dark and cold here in the winter. A similar plot of the average solar angle shows how the sun is always lower in the sky here, too.

Of course, the plot does not include the extra atmospheric attenuation from the sun coming in at lower angles. (The same thickness of clouds will require sunlight to travel further as a result.) That gives Ann Arbor and the other northern places an even more severe disadvantage.

Incidentally, the MATLAB code to display the months (where the x axis data is stored as days) is:
monthlengths = [0 31 28 31 30 31 30 31 31 30 31 30];
months = {'J' 'F' 'M' 'A' 'M' 'J' 'J' 'A' 'S' 'O' 'N' 'D'};
set(gca,'XTick',cumsum(monthlengths));
set(gca,'XTickLabel',months);
and I created the graphic with the same script I posted about earlier.

Trilinos on the mac

My configuration script for Trilinos on Mac OS X 10.5.5, for building without MPI:

./configure \
--cache-file=config.cache \
--prefix=/usr/local/trilinos-8.0.8-serial \
CC=/usr/bin/gcc CXX=/usr/bin/g++ F77=/usr/local/bin/gfortran \
--with-libs="-framework vecLib" \
--with-ldflags="-Wl,-multiply_defined -Wl,suppress" \
--enable-amesos \
--enable-epetra \
--enable-anasazi \
--enable-aztecoo \
--enable-examples \
--enable-didasko \
--enable-teuchos \
--enable-triutils \
--enable-galeri

I believe that the only software needed to get this to work is the gfortran compiler available at the mac High Performance Computing site.

I found that the with-ldflags line was necessary because otherwise the "pow" check would die from "Undefined symbols: "___gxx_personality_v0", referenced from: __gxx_personality_v0$non_lazy_ptr in ccJ9fzb2.o".

To get a trilinos-base program to compile, run g++ with the Apple compiler with options "-framework vecLib -I/usr/local/trilinos-8.0.8-serial/include -L/usr/local/trilinos-8.0.8-serial/" with the libraries (e.g. "-lepetra").

Publication-quality plots in MATLAB

I've done something similar (but much more simple), but this article has some interesting points. Now if only MathWorks would fix that stupid bug that showed up with the latest Java update.

Here's what I do: it scales the fonts and figures (using hgexport) and markers so that they show up as not to small a size when included as 3-inch figures, changes the font to Times, makes the minimum line width 0.5 points, exports to EPS, and uses the epstopdf tool to convert it to a pdf (since Matlab's pdf export is buggy and always puts the image on an 8.5 x 11 page). The results look consistent and very clean.
exportNiceFig.m

UPDATE 11/30/2008: I fixed the behavior of some of the resizing when there are multiple axes on a single figure.

UPDATE 01/21/2009: For future updates, see the projects page

Evolution personified

It strikes me as ironic when people personify evolution as some sort of designer. I was reading an article in SIAM News describing the complexity and mathematical modeling of finger movements, when I stumbled across this:

But nature, for no obvious reason, took a different route. Opting to locate the muscles in the palm and forearm, at a distance from the fingers themselves, evolution tinkered its way to a quasi-two-dimensional system of interconnected tendons, more like fishing net than fishing line. This decidedly nonlinear and complex solution, anathema to "reasonable" engineering design, has paid off handsomely.

In the 19th century or even early 20th century, I could easily imagine a scientific article marveling about the incredible design of the human body in the context of a divine Creator. But to use that kind of language when talking about random processes and natural selection seems kind of silly.

National Science Bowl 2007

Well, I've been in Washington, D.C., for two days now. Science bowl is a science competition for high schoolers in which I participated back in high school; it was probably my favorite part of all 12 grades.

This year I asked to come back as a moderator -- basically like Alex Trebek, except cooler, and there are a lot of us and we don't have moustaches and aren't on syndicated TV. But it's the same general idea. The actual competition starts tomorrow, although we moderators have spent most of today taking turns as moderators and competitors among ourselves (tons of fun, but still kinda exhausting over an 8-hour time period).

Yesterday and the previous evening we've been around the downtown area. I broke off from the rest of the group on Friday and walked around the museums on the mall. The Hirshorn modern art museum didn't have very good exhibits this time around (when I went two novembers ago it had some pretty cool stuff). I browsed quickly through the new Native American museum; its most interesting blurb was about an American Indian who dressed up in traditional clothes and lay down with his eyes closed in a glass box in a real museum as though he were an exhibit, as a rebellion against the idea that their traditions and heritage are only a thing of the past.

My favorite exhibit of anything in DC is by far the gems and minerals exhibit at the Smithsonian Natural History Museum. It is amazing.

Yesterday evening I met up with my friend who graduated and moved up here to work for the NRC -- we had a good time. The day before that I walked around the chevy chase neighborhood admiring the cherry and dogwood trees in bloom; it was pretty. Also, the house two doors North of the 4H campus has a Lamborghini in its driveway.

It's also weird to be around 350 high schoolers for most of each day. Especially when so many of them are so incredibly pretentious. A lot of them are brilliant, but that's still no excuse for the way they behave. The Secretary of Energy spoke to us on Friday morning -- one kid asked him a question about uranium reprocessing. The Secretary gave a very competent answer, but apparently it wasn't what the kid wanted to hear -- he actually started talking back to the guy about U-238. It must have taken incredible patience for the guy to say, "Yes. I know," and to continue unfazed. When you're speaking to someone that important who knows so much more than you do (not only has he been informed from whatever briefings he gets, but also he has a Ph.D. in chemical engineering or something) you need to ask your question, learn from his answer, and sit down. The kid wasted the time of both the Secretary and the rest of the audience, and he made himself look like an idiot in the process. Ugh.

But besides being around those kind of people part of the time (my roommate is sort of like that -- it's obvious he and his accomplishments are the only thing he cares about) it has been a lot of fun. Plus, my friend from high school and A&M is volunteering as well.

Smoretastic

The latest odd thing I did in my spare time was to create a little marshmallow-cooking furnace out of aluminum foil. The chocolate goes on a separate piece of foil that sits on top of it so that it's not rock-hard whenever it's bitten into. It makes a very even, light toast even without rotating it, since I suppose the incident thermal radiation is so uniform. Near-IR shows up as blueish purple to my digital camera, apparently.

CO2 sequestration

Today I came up with a great idea for burying our carbon dioxide.

Inject it as a high-pressure gas into liquid sugar. And sell it.

Deeelicous. I forgot how awesome pop rocks are.

Radium Ore Revigator

One thing I did while visiting Louisiana with family was go to a rural museum where they had lots of artifacts from the area from the '30s through the mid-19th century. There was one particularly interesting one (from a rad health point of view) was a "Radium ore revigator" which, according to this website, put vigor into "tired or wilted water," which, amusingly enough, isn't far from the modern-day quackery that you'll find in ads in the back of popular science.

Light pollution

Another issue that's been bugging me recently is the relatively great amount of light pollution that such a small town as College Station produces.

What is light pollution?

Light pollution occurs when an outdoor light shines not just at the ground but also into the air. An example would be the glowing spherical lights on a pole; roughly half of the light it emits goes straight up in the air.

What is wrong with light pollution?

Mainly two things.

1. It's a waste of energy. If the light pointed up were reflected down, that would increase the intensity of the light on the ground by roughly a factor of two; or, they could turn the power usage down to half and still get about the same brightness on the ground as the original lamp configuration.
2. It destroys the night sky. I'm not familiar with the exact numbers, but the magnitude of the stars that are visible in the sky drops dramatically. The weekend before last, I went to Enchanted Rock (a state park, plenty of distance away from bright cities), and when I took a look at the sky after it was dark, I was absolutely stunned by the sky (it had been a while since I had last been camping or far away from cities with a clear sky). At A&M, on the other hand, you look up on a clear night and you'll see a few stars and planets; on a cloudy day, you'll see a pink glow. From several miles away from the campus, you will still see a pink glow over where the lights are. This is not good.
   

What can I do to make light pollution go away?

Probably nothing. You could try petitioning the different people in charge, but sadly, chances favor their not caring. If they were to pay attention, you would tell them to modify their lamps so that they didn't shine straight up in the sky.

End of poorly assembled rant.