A High Throughput Workflow Environment for Cosmological Simulations

The Blind Cosmology Challenge (BCC) process will require generating multiple galaxy catalogs to the full photometric depth across the full 5000 square degrees of the DES survey. The effort will require roughly 6M SUs and generate 300 TB of output.

Competing effects drive our simulation requirements. On one hand, the dark matter distribution needs to be modeled within a large cosmic volume. On the other hand, galaxy surveys also sample the nearby population of dwarf galaxies, implying a need for high spatial and mass resolution. We have developed an approach that generates a set of discrete N-body sky realizations of dark matter structure spanning a range in resolution and volume. We then dress this dark matter distribution with galaxies brighter than the DES limiting magnitudes in each passband.

To implement the BCC, we plan to produce fifty $10^{10}$-particle N-body runs on XSEDE resources over the next two years. The collaboration has requested that BCC models explore a variety of cosmologies, parameters of which are known only to the Simulation Working Group members. After a blind processing period, constraints on cosmological parameters from the science teams will be compared against their true values, gauging the validity of the processing pipelines. We detail ongoing and proposed simulations along with our strategy for producing synthetic sky surveys of sufficient area and depth in the next section.

The dark matter structure from N-body simulations of a given cosmology forms the basis for galaxy catalog expectations. The N-body models we run on XSEDE resources store particle configurations, $\{x_{i}(t),v_{i}(t)\}$, in two forms: snapshots of the positions and velocities of all particles $i$ in the simulation volume at a fixed time $t$, and lightcones [12] that hold kinematic information for particles lying on the past lightcone of a virtual observer located at a fixed position, $x_{\rm cen}$, in the computational volume.⁹⁹9These particles satisfy the combined space-time requirement, $|x_{i}(t_{i})-x_{\rm cen}|=r_{\rm cosm}(t_{i})$, where $r_{\rm cosm}(t)$ is the cosmic metric distance as a function of time, a known function of the cosmological parameters. Once the N-body steps are completed, we proceed with additional processing using a combination of resources, including XSEDE, SLAC and other collaboration institutions. The post-processing aims to create a synthetic DES catalog of galaxy properties derived from the N-body lightcone output. Figure 1 shows a schematic representation of the post-processing workflow.

Halos, bound systems that host galaxies and clusters of galaxies, are identified in these outputs and their properties can be used to determine their central galaxy characteristics. Halos, as well as local density estimates in under-resolved locations, are used by the ADDGALS¹⁰¹⁰10for Adding Density-Determined GAlaxies to Lightcone Simulations algorithm to assign galaxy properties to suitably selected dark matter particles. The matter along the past lightcone also sets the gravitational lensing shear signal applied to these galaxies, and we are developing a new multi-grid, spherical harmonic algorithm for this computation. Finally, galaxy catalogs that include lensing shear signals are processed with telescope/instrument/noise effects to produce images expected from the back-end electronics of the Dark Energy Camera. From these images, which are generated for only 200 sq deg of sky, we can develop effective transfer functions to create full 5000 sq deg synthetic DES catalogs that contain realistic errors such as star/galaxy mis-classification, blended sources, and appropriate photometric errors.

Generating a cosmological N-body simulation has three main steps, described further below. The first two create an initial particle set consistent with the structure expected at an early time of a chosen cosmological model. The first step samples the linear perturbation spectrum of the cosmological model while the second step realizes that density and velocity field with a set of particles. For the latter, we use a second-order Lagrangian perturbation theory code (2LPTic) code that has been robustly tested in the community [10, 9]. The amount of CPU time required to make these initial conditions is small (typically 300-400 CPU hours for simulations of the scale described here), but the memory requirements are more substantial (slightly less than 2 GB per core), although easily achievable on XSEDE machines.

The final step evolves this particle distribution under its self-gravity within an expanding cosmological background. For this purpose, we are using a lean version of the Gadget cosmological N-body code [22, 23] modified by us to generate output along the past lightcone of synthetic observers in the computational volume. We have worked on XSEDE resources for the past year to modify and optimize the code. The lean version has significantly reduced memory overhead, 44 bytes per particle compared to 84 for standard Gadget. This reduction allows the simulations to fit on a smaller number of processors, affording better scaling. These simulations require 50–100k CPU hours depending on the parameters and resolution of the computation, and generate up to 10 TB of total data output.

Once the simulation codes are deployed on XSEDE computational resources, we register descriptions of these applications with the Apache Airavata registry service. These descriptions are used by the Airavata GFac component to generate the artifacts required to expose the application as a service. The workflow developer can access these wrapped application services and construct workflows and orchestrate executions on target compute resources. The resulting workflow abstractions reduce human inefficiencies by providing a uniform interface for the scientist and hiding unnecessary complexities.

To illustrate the construction of a cosmological workflow, we will describe the process in developing the N-body simulation workflow of the process illustrated in Figure 1. Firstly the nature of the applications, its execution characteristics, and its input and output data are analyzed. The application meta information, including the executable location, its nature like serial or MPI, inputs and outputs, are described and registered with Airavata registry. This process was followed for the following four applications.

BCC Parameter Maker

This initial setup code is written as a python script and prepares necessary configurations and parameter files for the workflow execution. This simple script is forked on the XSEDE Ranger job management nodes.

CAMB

The CAMB (Code for Anisotropies in the Microwave Background) [14] application computes the power spectrum of dark matter, which is necessary for generating the simulation initial conditions. This application is as a serial fortran code. The output files are relatively small ASCII files describing the power spectrum.

2LPTic

The Second-order Lagrangian Perturbation Theory initial conditions code [9] (2LPTic) is programmed using Message Passing MPI C code that computes the initial conditions for the simulation from parameters and an input power spectrum generated by CAMB. The output of this application are a set of binary files that vary in size from $\sim$80–250 GB depending on the simulation resolution.

LGadget

The LGadget simulation code is MPI based C code that uses a TreePM algorithm to evolve a gravitational N-body system [22, 23]. The outputs of this step are system state snapshot files, as well as lightcone files, and some properties of the matter distribution, including diagnostics such as total system energies and momenta. The total output from LGadget depends on resolution and the number of system snapshots stored, and approaches close to 10 TeraBytes for large DES simulation volumes.

After all the above applications are registered, the Blind Cosmology Challenge Workflow is constructed using Airavata XBaya. The resultant workflow graph is shown in Figure 2. The workflow provides capabilities to configure the generation of initial conditions as well as the full N-body simulation components of the BCC process.

Iterative execution support for long running applications The N-Body simulation requires multiple days of execution, but the XSEDE Ranger cluster limits maximum wall time of 48 hours. To mitigate this limitation, the workflow infrastructure has to allow iterative support so the job can be broken down into multiple increments of 48 hour jobs harnessing the check-point restart capabilities within the application. These capabilities required sophistication beyond the blind restarts, in order to account for application execution patterns and exception handling. These capabilities can be potentially matured into a formal Do-While construct semantics of workflow engines.

Output Transfers The workflow executions tend to produce terabytes of data residing on the cluster scratch file systems and have to be persisted for a longer durations. The data movement to archival systems like TACC Ranch for long term storage have to be provided. The large file data movement is non-trivial process. Even though many advancements have been made in this area, the seamless reliable data transfers are still challenging. The emerging solutions like Globus Online [1], GridFTP client API[2] and bbcp [3] are potentially viable options. We leverage the Ranch archival system mounting on Ranger and use bbcp to copy the workflow outputs for test runs. We have yet to explore this as a solution for production executions. To transfer data to the post processing remote locations, Globus Online and GridFTP client are viable options.

The XBaya workflow shown in Figure 2 was tested and refined using smaller simulations over the period Oct 2011 to Mar 2012. We transitioned to production use for our most recent simulations. One workflow-managed simulation has run to completion, while a second job crashed because of a hardware problem on TACC Ranger.

Even with this limited information, we can compare the efficiency of running the required jobs under XBaya to those previously run manually. Jobs were submitted to the long queue at TACC Ranger, which has maximum resource limits of 1024 processors and 48 hour runtime. In Table 2, Total Time is the wallclock time interval for the entire production process, while CPU Time gives the sum of the run times of the required jobs. Times reflect the full N-body production process, from generating initial conditions all the way through to completing of the final N-body timestep. Efficiency is the ratio of CPU to Total times, with $100\%$ representing the ideal scenario of running without interruption.

The first two rows of In Table 2 list different manually-processed simulations. The first is a large simulation that needed four total submissions to the queue: one for initial conditions, and three for N-body computation. The second row was a smaller job requiring three submission, one for initial conditions and two for N-body. These runs are relatively inefficient because, each time the wallclock limit is reached, the user is notified via email, and then must log back into the cluster to submit the next job. If the wallclock limit is reached at an inconvenient time, considerable time can elapse before the next submission. More submissions tends to drive up the inefficiency.

The last row shows the efficiency of running the full production process via the XBaya workflow environment. By enabling immediate submission of jobs when the preceding job finishes, the efficiency improves to 92%, well above the 50–60% found for manual processing. Under the workflow, the only time spent not computing is spent waiting in the queue time on the cluster. That is, our job production becomes limited only by the instantaneous compute resources available on TACC Ranger.

We found that the workflow can also help prevent errors in simulation set-up. Our first production level workflow was designed to have the same parameters as a pair of manually completed simulations. We soon found that the latter, ‘by hand’ simulations were inconsistent with the workflow simulations. Investigating the source of the inconsistency, we found that the workflow was correct, and that a parameter had been mistakenly set to a wrong value in the original simulations. This shows the workflow’s value in reducing the risk of simple human run-time errors.

Synthetic surveys from the first set of simulations are currently in use by DES science groups. We provide here two examples of analyses from cluster and weak lensing science groups.

Members of the Galaxy Cluster Working Group are using the synthetic galaxy catalogs to evaluate different methods for identifying the massive halos that host clusters. Because 5-band optical photometry provides relatively crude distance information for each single galaxy, the ability to identify localized spatial clusters is compromised by poor depth resolution. Different cluster finding algorithms have methods to mitigate this loss of information, but none has been applied to a large survey with the depth of DES.

Model galaxies are assigned to specific, unique dark matter halos, so an ideal cluster finder would return the list of original halos with all members intact. In practice, the distance errors, incorrect choice of cluster location and other effects imply that the set of clusters found does not perfectly match the input set of halos. However, a correspondence map between cluster and halo sets can be derived using joint galaxy membership criteria. The utility of a cluster finder can then be characterized by two parameters: i) purity, the fraction of clusters that correspond to genuinely massive halos, and; ii) completeness, the fraction of halos that are found by the cluster finder. The ideal value for both parameters is one.

Figure 4 shows recent purity and completeness measurements from the redMaPPer cluster finder [19]. The redMaPPer cluster finder works on so-called red-sequence galaxies, galaxies with evolved stellar populations that tend to be found in massive halos. The algorithm is a Bayesian method that assigns galaxies a cluster membership probability based on the galaxy’s color, and the density of nearby galaxies of a similar color. Clusters are identified as clumps of similar high-likelihood member galaxies. Figure 4 shows that redMaPPer is $>80\%$ pure and has similarly high completeness above a redshift-dependent minimum halo mass.

In a related project, members of the Weak Lensing Working Group have been investigating how best to estimate the mass of an observed cluster based on the weak lensing distortion of background galaxies. Light from faint background galaxies is bent by the gravitational field of a massive halo, distorting the shape (shear) and density (magnification) of the galaxies. Using a maximum likelihood estimator [21] to fit the synthetic observations, and comparing to the true halo masses from the simulation, the group can understand how best to estimate masses and also provide feedback on the model used to generate the synthetic observations. A recent calibration plot comparing estimated and true masses is shown in Figure 5.

Implementing the Airavata workflow for this project has entailed some overhead. Scripts that set up the input parameter files needed to be developed, and new features were added to the existing codebase so that our applications would integrate more effectively with the workflow framework. Interaction between the co-authors of this document—domain scientists along with TeraGrid AUSS team members—was essential to achieving a production-level service. The effort invested has been worthwhile, with a significant gain in realized efficiency (Table 2) for our first production run.

In the near term, we have applied for continued XSEDE ECS support aimed at integrating the postprocessing and catalog production stages into the Airavata workflow (see Figure 1). We also want to integrate data movement into the workflow. In addition, we have requested time at two XSEDE facilities (TACC Ranger and UCSD Trestles), and would like to expand the workflow logic to choose execution location in real time, based on queue loads.

We also would like to use Airavata to improve our provenance practices. The workflow system can capture provenance, including information such as when the data set was created, by whom, where, and with what application version and which input parameters. Improved provenance can enable broader forms of sharing, reuse, and long-term preservation of our simulations and the resultant galaxy catalogs. Additional development could include standardizing an API for simulation parameter input and output, so that other codes could be easily implemented in the workflow, such as N-body models that employ modified gravity [17, 15].

In the longer term, we could also expand our scope, generalizing our galaxy catalog construction process into a science gateway that would support broader classes of astrophysical studies. The optical catalogs we create could be augmented by synthetic surveys at other wavelengths, from radio to X-ray, and our focus on galaxies could be expanded to include quasars, galactic stars, and other astrophysical objects.