Open Source Development of Community Monte Carlo Model

We are in the process of developing a Monte Carlo radiative transfer model that makes state-of-the-art techniques available to a wide range of users. The initial release of this code is available here. Questions or comments about the code and the development effort can be addressed to Robert Pincus.

The approach

Monte Carlo models get used in an enormous range of applications. Because the models are so computationally intensive, investigators build codes to answer the specific question at hand. We do not intent to build a model that to meet everyone's needs. Instead we will build a collection of interlocking code pieces, including a set of Monte Carlo integration routines. The core integrators will do the unpolarized narrowband problem (i.e. one set of scattering properties) and each can be designed to compute certain quantities as efficiently as possible.

We will put together a collection of independent, reusable modules that individuals can use to build a model for a particular application. These software Legos will speed development by accomplishing common tasks in a robust, efficient manner. The integrators are one of the Legos, but we'll need lots more. We hope to develop a large enough collection of packages to be useful in solving a wide range of specific problems (e.g. average fluxes in a 2D domain; heating rates in a 3D domain).

Computational Legos

As a concrete example of a software component, consider the treatment of surface reflectance in a narrow-band simulation. We usually assume that the surface is either black or Lambertian, but more generally we might want to include BRDFs for vegetated surfaces, or allow the possibility of a spatially variable surface. What the Monte Carlo integrator needs, then, is a way to say "a photon strikes the surface at this location and this angle; what is the weight and direction of the reflected photon?"

We will develop one or more package(s) to treat surface reflectance. Each package needs one subroutine in which the surface reflectance properties are specified as a function of position, and a second subroutine in which photons encountering the surface are mapped to reflected photons. Calls to the first subroutine will be simple or complicated, depending on the complexity of the reflectance model. But only the second subroutine gets called by the Monte Carlo integrator, and the calls made by the Monte Carlo integrator are the same for every implementation of the surface reflectance package. This will let us write a Monte Carlo integrator that can handle any kind of surface reflection model transparently. If you're tired of using a uniform surface, replace the uniform surface reflectance Lego with a spatially variable one.

Formally, what we will do is define a standard interface (subroutine calls) for each software component. The implementation underlying the interface will vary depending on the application.

A software shopping list

We require (at least) the following components:

• Status handling - a way to pass information about the success or failure of a task between subroutines. Every package we write can use the same standard way of communicating this information. An example implementation is available at the UW CIMSS web site.
• Random number generation. Fortran 90 provides random number generation, but it's a little clumsy and will take some massaging to work in a parallel environment. We will devise an interface that suits our needs explicitly.
• General input and output routines for fields. It might be useful to have routines that read and/or write gridded 1, 2, or 3D fields, for example. Some folks will want to see ASCII files; others will want netCDF or HDF files. Eugene is working along these lines already.
• Parameter input. This would include things like the number of photons, the names of input and output files, etc. One version of this routine might read from standard input, another might have a clickie-poo graphical interface; if we're careful in our design we'll be able to use whatever pleases us transparently.

To describe the optical properties of the domain we'll need

• A description of the scattering medium. An integrator should be able to ask about the size of the domain, and get the optical properties (extinction, single scattering albedo, phase function, and equal probability angles) at any location. Different implementations will need to be written for 1, 2, and 3-D domains. The integrator would use the package to compute the extinction along a path and to find the new direction of scattered photons.
• A description of the absorbing medium. K-distributions can be done either inside the integrator or post-facto, if one keeps the path lengths. We should design a package that will support either approach.
• Surface reflectance, as described above. We'll need to be able to ensure that the spatial domains are the same in each of these three modules.

It may make more sense to use a single software component to describe scattering, absorption, and surface reflection; this would let us use the idea of a spatial domain without knowing details about which processes are being included. This is open to discussion.

And of course we need

• Integrators. There will be at least 4 integration routines, one for each permutation of 2D/3D with radiance/flux computations. I'm sure there will be others. Among other things, we'll want to tell the integrator at what spatial locations we want the output. We can probably do this generally, so we can specify a single point if we want domain averages or a set of points if we want spatial fields. Integrators should certainly provide uncertainty estimates. Someone probably wants to write an integrator that runs on a massively parallel machine.

Are there other common Monte Carlo tasks that can be broken out like this? Pipe up if you've got ideas. We ultimately want to provide enough components to do the whole job.