Real science at the petascale

Radhika S. Saksena, Bruce Boghosian, Luis Fazendeiro, Owain A. Kenway, Steven Manos, Marco D. Mazzeo, S. Kashif Sadiq, James L. Suter, David Wright, Peter V. Coveney
Published 17 May 2009.DOI: 10.1098/rsta.2009.0049

Abstract

We describe computational science research that uses petascale resources to achieve scientific results at unprecedented scales and resolution. The applications span a wide range of domains, from investigation of fundamental problems in turbulence through computational materials science research to biomedical applications at the forefront of HIV/AIDS research and cerebrovascular haemodynamics. This work was mainly performed on the US TeraGrid ‘petascale’ resource, Ranger, at Texas Advanced Computing Center, in the first half of 2008 when it was the largest computing system in the world available for open scientific research. We have sought to use this petascale supercomputer optimally across application domains and scales, exploiting the excellent parallel scaling performance found on up to at least 32 768 cores for certain of our codes in the so-called ‘capability computing’ category as well as high-throughput intermediate-scale jobs for ensemble simulations in the 32–512 core range. Furthermore, this activity provides evidence that conventional parallel programming with MPI should be successful at the petascale in the short to medium term. We also report on the parallel performance of some of our codes on up to 65 636 cores on the IBM Blue Gene/P system at the Argonne Leadership Computing Facility, which has recently been named the fastest supercomputer in the world for open science.

1. Introduction

Ranger (http://www.tacc.utexas.edu/resources/hpcsystems/) comprises 15 744 AMD ‘Barcelona’ quad-core processors housed in 3936 nodes, with a total memory of 123 TB. It has a total shared disk space of 1.73 petabytes (PB) and a local disk of 31.4 TB with global, high-speed Lustre file system, running across 72 I/O servers. This massive I/O capacity is critical for our computational endeavours. Ranger is connected to the US TeraGrid (http://www.teragrid.org/userinfo/hardware/resources.php) backbone via four nodes with 10 Gb Ethernet adapters. It has a theoretical peak performance of 0.58 petaflops (Pflops) and, according to recent benchmarks (TOP500, http://www.top500.org/list/2008/06/100), it delivers a Linpack speed of 0.31 Pflops. More recently, we have begun to use the IBM Blue Gene/P system called Intrepid at the Argonne Leadership Computing Facility (ALCF; http://www.alcf.anl.gov/resources/storage.php), Argonne National Laboratory (ANL), and have found excellent scaling on up to 65 536 cores. Intrepid (http://www.mcs.anl.gov/uploads/cels/files/news/2008/ALCF%20Science%20Brochure.pdf) delivers a theoretical peak of 0.56 Pflops and a Linpack speed of 0.45 Pflops, making it the fastest computer available for open science research. Intrepid consists of 163 840 cores and 80 TB of RAM. It has a total storage capacity of over 6 PB with an aggregate transfer speed of 80 GB s−1. The storage is managed by two parallel file systems—PVFS and GPFS. Networking requirements for our work are met by a dedicated 1 Gb JANET (http://www.ja.net/services/lightpath/index.html) lightpath link (Coveney et al. in press) between the Starlight switch/router facility in the USA and our local storage and visualization resources at University College London. In the following sections, we will describe each of our applications that have been ported to these petascale machines and the scientific research that has been made possible. We first describe the lattice–Boltzmann (LB) code HYPO4D that has been ported to both these petascale resources. Next we describe our work on large-scale simulations of liquid crystals performed using an LB code called LB3D. We then discuss our large-scale molecular dynamics (MD) simulations of clay–polymer nanocomposite systems using LAMMPS (http://www.cs.sandia.gov/∼sjplimp/lammps.html). Finally, we present two biomedical applications that use the MD code NAMD and the LB code HemeLB, respectively. The applications, which all ran concurrently, have been enabled by the petascale resources mentioned here and involve complex workflows and high-throughput computing.

2. Identification of unstable periodic orbits in the Navier–Stokes equations

The unstable periodic orbits (UPOs) of a dynamical system provide us with a fertile representation of the chaotic dynamics on the attractor of that system and allow for the computation of time averages in terms of an exact expansion, based on the dynamical zeta function (Artuso et al. 1990; Cvitanović et al. 2004). A small number of these UPOs should furnish accurate predictions of turbulent flow. The averages obtained in this way are not stochastic in nature; also there is no need to solve the initial-value problem each time we wish to compute a given observable. One of the more relevant applications of this method lies in turbulent hydrodynamics, described by the Navier–Stokes equations (NSE; Temam 1984), and hailed as one of the great unsolved problems of classical physics.

The approach we adopt in this work (Boghosian et al. 2008) is one that efficiently parallelizes time and space, by locating candidate space–time orbits and then applying a fully four-dimensional numerical relaxation algorithm. The algorithm has been successfully tested on low-dimensional dynamical systems such as the Lorenz equations, with favourable results (Boghosian et al. 2008). For the case of the three-dimensional NSE, this approach puts the computational requirements into the petascale as we discuss here.

Within our research, we have developed a scalable LB (Succi 2001) code called HYPO4D, an application to locate and characterize UPOs in the NSE. The code is written in C and parallelized using the Message Passing Interface (MPI), according to a spatial domain decomposition strategy. Asynchronous communications are used to prevent deadlock. We consider periodic boundary conditions so as to exclude wall effects and have an isotropic system. HYPO4D has been tested and deployed on several UK National Grid Service (NGS; http://www.ngs.ac.uk) sites, on HPCx (http://www.hpcx.ac.uk) and on TeraGrid resources with special emphasis on Ranger. Figure 1a shows the linear scaling of HYPO4D on Ranger, with the total system size being simulated kept fixed (‘strong scaling’) as well as with the portion of the system per core kept fixed (‘weak scaling’). Figure 1b shows weak scaling benchmarks, performed on the Blue Gene/P machine at ALCF.

Figure 1
Figure 1

HYPO4D performance data, the number of total site updates per second (SUPS) given as a function of the number of cores used on (a) Ranger (solid curve, fixed total system size: Lx=Ly=2048, Lz=1020; dashed curve, fixed system size per core: 32×32×32) and (b) Blue Gene/P machine at ALCF (solid curve, 32-cube per core; medium-dashed curve, 64-cube per core; short-dashed curve, 128-cube per core).

In the simulations performed so far, the fluid starts from a condition of constant density and zero velocity, and an ABC-type force (Dombre et al. 1986) is then applied at each time step to drive the viscous flow. After having identified some relevant and physically valid sets of parameters on the NGS, we performed simulations of millions of LB time steps on Ranger, and then compared many of these in order to identify minima that could indicate the proximity to a periodic space–time trajectory. One of these minima can be seen in figure 2 for a cubic lattice with dimension L=64 and a Reynolds number of 500, which corresponds to weakly turbulent behaviour. The quantity T gives the distance between the two time steps being compared, while t indicates the value of the time step in LB units. The quantity being plotted isEmbedded Image(2.1)where fi(r,t) are the l distribution functions at each site r of the lattice at time step t and the sum is over the whole lattice. HYPO4D uses the D3Q19 LB model, with l=19. Larger lattices were also simulated on Ranger, including cubic lattices of dimensions L=128 and 512. For the latter, a total of 3.4 million LB time steps were simulated with checkpoints being kept at every 104 steps, for post-processing. This took close to half a million CPU-hours to run, and produced approximately 10 TB of data.

Figure 2
Figure 2

Three-dimensional plot of a minimum centred around t=9.9×105 and T=3×104. The values of Δ(t,T) appear on the z-axis. Darker regions indicate greater similarity between different time steps.

In this novel approach to turbulence, memory is the main bottleneck. For the minimum shown in figure 2, the relaxation procedure will require approximately 1.2 TB for a dataset consisting of 3×104 time steps. On increasing the system size to simulate flows with higher Reynolds numbers, we easily reach the level of petascale resource requirements, memory-wise. The results obtained so far show that the best candidates for even the smaller UPOs will consist of at least several tens of thousands of time steps. We intend to simulate flows as close to the fully developed turbulent regime as possible. Preliminary calculations show that we can aim to identify UPOs in systems with Reynolds numbers as high as 104, with the memory requirement for a single candidate orbit being approximately 0.5 PB. By identifying and cataloguing these UPOs, we will then be able to systematically compute observables, in a non-stochastic manner, which could have important consequences not only in turbulence studies, but also in other areas involving highly nonlinear dynamical systems.

3. Large-scale LB simulations of liquid crystalline rheology

The non-equilibrium behaviour of amphiphilic liquid crystals self-assembled into various periodic phases, such as the lamellar, hexagonal and the cubic diamond, gyroid and primitive phases, is commonly observed in experiments. These amphiphilic liquid crystalline phases are of great technological and scientific interest in materials science, in the synthesis of mesoporous catalytic templates, in biosensing, ultrafiltration and protein crystallization, and as possible intermediates in various important cellular processes (Schwarz & Gompper 2002; Squires et al. 2002). The non-equilibrium dynamics of these liquid crystalline phases can provide insight into the pathways and kinetics of self-assembly in non-cubic↔cubic as well as cubic↔cubic phase transitions that are also being actively studied experimentally (Seddon et al. 2006). Shear flow fields as well as other fields such as roll casting and extensional flow are used in the processing of functional materials for switchable optical, transport and electronic properties as well as ordered nanostructures (Chen et al. 1997). To design functional materials, one must understand the physical basis of the effects of molecular architecture and processing conditions (especially shear frequency, strain and temperature) on the rheological response of the liquid crystalline mesophases (Berghausen et al. 1998; Rey & Denn 2002).

The distribution of defects within a liquid crystal domain, defined as the ‘texture’ of the liquid crystal, plays a significant role in its rheology (Rey & Denn 2002). Kinetic LB simulations on smaller 643 lattice sites lead to self-assembled monodomains of perfect gyroid, diamond and primitive mesophases (Saksena & Coveney 2008a), while real-world liquid crystalline mesophases usually do not exist as monocrystals but rather as multiple oriented grains with grain boundary defects as well as dislocations and defects within the grains. Investigation of the rheology of these defect-containing multi-domain systems requires simulation of larger lattices, ideally comprising 2563 lattice sites or more. Routine production simulation of such large lattices has recently become possible with the arrival of the Texas Advanced Computing Center (TACC) Ranger machine. The simulation of ternary amphiphilic mixtures composed of oil, water and an amphiphilic species is performed using a kinetic LB implementation in the highly scalable, parallel MPI code LB3D (Love et al. 2003). The code has built-in computational steering capabilities, has been hosted in the Application Hosting Environment (AHE; Coveney et al. 2007) and has been deployed on various grid resources in the USA, UK and EU. Here, we report on the close-to-linear parallel scaling behaviour of LB3D on Ranger and on the new UK Cray XT4 machine, HECToR (http://www.hector.ac.uk/), in figure 3.

Figure 3
Figure 3

Lattice SUPS versus number of cores calculated from a 100 time-step simulation on different lattice sizes using the LB3D code on (a) Ranger (circles, 1024×1024×1024 lattice; squares, 2048×2048×1024 lattice; diamonds, 2048×2048×2048 lattice) and (b) HECToR (circles, 1024×1024×1024 lattice; squares, 2048×1024×1024 lattice).

The simulations of the rheological yield curve of a 2563 defect-containing multi-domain gyroid system on Ranger have given us insight into the defect dynamics under shear flow for the first time (Saksena & Coveney 2008b). We found that regions of high stress magnitude arising from the applied strain are localized in the vicinity of defects (figure 4) and we are currently studying the temporal evolution of the defects and the associated stress field. Additionally, we have observed non-Newtonian shear-thinning behaviour (also seen in smaller 643 and 1283 lattice systems) as applied strain rate was increased (Giupponi et al. 2006). The strain velocity profile induced by the Lees–Edwards boundary conditions takes approximately seven times longer to be established, i.e. 70 000 LB time steps from the time shear is first applied, as compared with the smaller 1283 system. This increase in simulation time required to achieve steady state in shear flow simulations with increasing system size further contributed to making such large-scale rheological simulations computationally unfeasible before Ranger became available. Future work in this area will focus on following temporal defect dynamics that involves extremely disk-intensive analysis; scaling of disk I/O will become a major factor in determining code performance as we run simulations on higher and higher core counts. We plan to extend this work to other liquid crystalline phases such as the primitive and diamond phases (Saksena & Coveney 2008a).

Figure 4
Figure 4

A 323 sublattice extracted from a 2563 lattice site simulation domain. Contour plot of the gyroid morphology is superimposed on the volume rendering of the pxz scalar field of the component of the pressure tensor in the shear plane, pxz=−σxz. The maximum stress value σxz, corresponds to the defect close to the centre of the domain.

4. Material properties of clay–polymer nanocomposites

In the materials science domain, we are using large-scale MD simulation to understand the unique chemical and physical aspects associated with polymer-based nanocomposite materials composed of micrometre-sized clay platelets dispersed in a polymer matrix (Pinnavaiam & Beall 2000). These nanocomposites have similar, or improved, mechanical, thermal, barrier and electrical properties to conventional composites, but at considerably lower filler fractions. A fundamental understanding of their properties will lead to widespread usage as with traditional, conventional composite materials (Chen et al. 2008). It is difficult to access the material structures of the clay–polymer nanocomposites experimentally, owing to their small grain size, while at the same time it has proved problematic to predict their materials properties theoretically owing to the great disparity of length scales inherent in these systems. Hence, valuable insights into the structure, behaviour and properties of clay–polymer nanocomposites are provided by atomistic simulations that are large enough to sample all the length scales required to describe the system accurately.

In order to fully exploit petascale resources such as Ranger, we use the open-source, MPI MD code LAMMPS. We find almost linear scaling of LAMMPS on up to 4000 processors with a system approaching 100 million atoms (figure 5). The deviation from linearity depends on the system size, as communication begins to dominate over processing on individual cores.

Figure 5
Figure 5

Scaling properties of LAMMPS on Ranger at TACC, calculated using a clay system containing 85 095 000 atoms. Speed-up=time(N1)/time(N), where N is the number of processors used in the simulation and N1 is the reference simulation, using 256 processors.

We have used LAMMPS to investigate the effects of montmorillonite clay filler on poly(ethylene) glycol in clay–polymer nanocomposites (Suter & Coveney in press). The sizes of these simulations were approximately 3 million atoms, two to three orders of magnitude larger than any previous atomistic clay mineral study. This was facilitated through a global infrastructure of federated grid-based computing resources (the US TeraGrid, the UK NGS and the EU DEISA (http://www.deisa.eu) grids) accessed through our lightweight hosting environment, the AHE (Coveney et al. 2007). We find that the conformation of polymer adsorbed on the surface of the clay is distinct from that of polymer far from the surface, owing to complexation with the counterions found near the charged clay surface. The diffusion of the adsorbed polymer is much smaller, which effectively increases the size of the nanofiller and renders gas permeation more tortuous. We also performed non-equilibrium MD simulations by imposing a strain on the model and analysing the stress response. By partitioning the stress response into clay and different polymer layers, each a certain distance from the clay surface, we find that the adsorbed polymer layer is much stiffer than the other layers of polymer molecules in the system. Our next step is the simulation of physically realistic-sized atomistic models, comprising 100 million atoms, which explicitly include the edges of the clay for the first time, effectively reproducing the experimental system. These simulations will run optimally on 4096 cores as per the scaling performance on Ranger (figure 5). In order to simulate polymer penetration into the clay platelet, which occurs on a time scale outside that possible with a single MD simulation, we will use a parallel tempering technique (replica exchange), which involves simulating multiple ensembles of the system, at different temperatures, on multiple partitions of cores. These replicas are then exchanged according to a Monte Carlo procedure. Since the exchanges occur rarely, parallel tempering can be regarded as an embarrassingly parallel problem. For such a scenario, the entire study will comprise at least 16 replicas of the clay–polymer nanocomposite system each running on 4096 cores for an estimated wall-clock time of two weeks. This unprecedented scale is therefore eminently possible with petascale resources.

5. Drug resistance in HIV-1 proteases and reverse transcriptases

An emerging concept in the biomedical domain is that of the use of high-end computational power to confer clinical decision support for the treatment of a range of illnesses (Sadiq et al. 2008b; Sloot et al. 2008). Applied to HIV therapy, a new approach that uses this paradigm aims to compute the binding affinity of any given inhibitor to any sequence of the corresponding viral protein target with sufficient accuracy to discriminate resistant mutations on a patient-specific basis and with enough rapidity to return on clinically relevant time scales. Therefore, provided the method fulfils such requirements and that suitable computational power is available, the binding affinity of a large array of protein–ligand combinations can be calculated for any given patient and in turn used to optimize treatment for that individual.

We have previously investigated the applicability of binding free-energy methods in fully atomistic molecular simulations as a way of ranking drug-resistant mutations of the inhibitor saquinavir (SQV; Stoica et al. 2008). These used the molecular mechanics Poisson–Boltzmann surface area (MMPBSA; Wang & Kollman 2001; Wang et al. 2001) with the inclusion of configurational entropies from normal mode analysis (Brooks et al. 1995) and gave excellent correlation to experimental values (Maschera et al. 1996; Ermolieff et al. 1997). We have also developed a tool, the binding affinity calculator (BAC), for the automated construction of MD simulations of HIV-1 protease–ligand complexes and calculation of binding free energies (Sadiq et al. 2008a), which further generalizes the protocol to any array of mutations and available ligands.

Here, we used BAC to perform an unprecedented number of fully atomistic simulations of an array of HIV-1 protease variants (wild-type and mutants) with six FDA inhibitors as well as seven cleavage-site natural peptide substrates using both single- and three-trajectory MMPBSA approaches up to a production time of 50 ns each. Peptide substrates were the MA-CA, CA-p2, p2-NC, NC-p1, p1-p6, RT-RH and RH-IN with amino acid sequences SQNY-PIVQ, ARVL-AEAM, ATIM-MQRG, QANF-LGKI, PGNF-LQSR, AETF-YVDG and RKIL-FLDG, respectively (Prabu-Jeyabalan et al. 2002). The G48V, L90M and G48V/L90M mutant proteases bound to the NC-p1 and CA-p2 substrates and the V82A, A431V and V82A/A431V mutations in complex with NC-p1 were also simulated. Inhibitor systems simulated were amprenavir (APV), indinavir (IDV), lopinavir (LPV), ritonavir (RTV) and SQV with the HXB2 wild-type protease sequence (WT) and five MDR mutants, including the hexa-mutant L10I,M46I,I54V,V82A,I84V,L90M (denoted HM) as well as comprising subsets of HM, namely the M46I,I54V,V82A,I84V quatro-mutant (denoted QM) and pairwise mutants L10I,L90M (denoted DM), M46I,I54V (denoted FL) and V82A,I84V (denoted AS). The equilibration period for each run was always 2 ns, followed by a production run for each system up to 10 ns. LPV and RTV sets of systems were extended to 50 ns and the separate three-trajectory MMPBSA method was additionally applied to the LPV system, each up to 50 ns production.

Multiple and/or ensemble MD approaches, employed by projects such as Folding@home (Shirts & Pande 2000; Fujitani et al. 2005; Jayachandran et al. 2006), have been shown to explore conformational space more effectively than single long trajectories. Such approaches, in which multiple replicas of the same system are simulated with slightly varying initial conditions, are particularly suited to petascale resources owing to the overwhelming speed-up in turnaround they can afford over performing single long trajectory simulations of equivalent total simulation time. We simulated all LPV systems using the ensemble MMPBSA method, in which 50 replicas of each LPV-bound system were simulated each to 4 ns production.

The details of the underlying simulation protocol and free-energy calculation method are provided in Sadiq et al. (2008a) and Stoica et al. (2008). Over 3.8 μs of simulation time was achieved within two weeks of elapsed wall-clock time via distribution across a maximum of 150 simultaneously run systems, each system comprising approximately 40 000 atoms and using 32 CPUs. For the first time, we were also able to expand the methodology to the reverse transcriptase system (an approx. 160 000 atom system). The well-established MD code NAMD was used for the production runs, giving a benchmark performance of approximately 8 h ns−1 under the above conditions. Minimum single-trajectory turnaround time (TAT) using 192 cores (at 18 ns d−1) was 416 hours; theoretical minimum ensemble trajectory TAT using 9600 cores (simulating 300 replicas simultaneously at 900 ns d−1) was 48 hours. In practice, we achieved on average 120 ns d−1 with a peak of 300 ns d−1. The theoretical minimum total TAT (simulation+calculation) using this approach was approximately 3 days. In total, over 3.5 TB of trajectory and free-energy analysis data were generated. Importantly, the decoupled nature of each replica means that, using the entire core count of Ranger, approximately 6 μs of MD production simulation data can be generated per day. The protocol used here is thus scalable to any number of cores, being able to turn around approximately 2000 replicas (40 HIV-1 protease–ligand binding affinities) within 3 days on Ranger.

Figure 6 shows the relative ranking of the binding free energy of LPV and RTV studied when bound to wild-type and MDR proteases over a sampling time scale of up to 50 ns. The LPV and RTV systems exhibited excellent peak correlation coefficients (κ) of 0.84 and 0.93, respectively. The results are encouraging and they build on our earlier work (Stoica et al. 2008).

Figure 6
Figure 6

Computed binding free-energy differences compared with experimental results for wild-type and all MDR proteases with inhibitors (a) LPV (black circles, experiment; red squares, 1–10 ns, κ=0.81; green diamonds, 1–20 ns, κ=0.72; blue uptriangles, 1–50 ns, κ=0.84; purple left triangles, 31–50 ns, κ=0.79; orange down triangles, 41–50 ns, κ=0.84) and (b) RTV (black circles, experiment; red squares, 1–10 ns, κ=0.93; green diamonds, 1–20 ns, κ=0.86; blue up triangles, 1–50 ns, κ=0.59; purple left triangles, 31–50 ns, κ=0.40; orange down triangles, 41–50 ns, κ=0.46) using 50 ns single-trajectory runs. For comparison, the default experimental correlation is also shown. Correlation coefficients, κ, between ΔGtheor and ΔGexp for each set are as listed.

6. Whole-brain patient-specific haemodynamics

Understanding cerebral blood flow behaviour is crucial in the diagnosis and treatment of cardiovascular disease. Blood flow simulation in combination with patient-specific medical imaging data can help not only to understand blood flow patterns in the presence of neurovascular pathologies, but also to provide a tool that surgeons can use in a diagnostic and therapeutic capacity for patients according to their unique vascular structure. HemeLB is an application specifically designed to address fluid flow in sparse geometries, which simulates pulsatile blood flow in the brain (Mazzeo & Coveney 2008). Performance benchmarks show that the core LB code of HemeLB scales superlinearly (figure 7a) on some platforms. The HemeLB algorithm uses a data layout that automatically handles lattices with large unstructured regions covered by non-fluid nodes. A single-pass algorithm is used to spatially decompose the vascular tree into partitions of equal size—resulting in balanced communication between MPI ranks. Partial overlapping of communication between subdomains and on-core computation is also exploited to enhance performance, and the computational kernel is highly optimized. Indeed, the parallel performance of this code outstrips any comparable LB codes that we are aware of (Mazzeo & Coveney 2008). HemeLB has been used in combination with MPIg, the latest pre-release of the successor to MPICH-G2 (Karonis et al. 2003) on the TeraGrid, Louisiana Optical Network Initiative (LONI; USA) and NGS2 (UK) networks, achieving near-linear scalability across different computational sites (figure 7b).

Figure 7
Figure 7

(a) Scalability of the LB fluid solver (FS), along with the FS in combination with real-time volume rendering (FS+VR) of HemeLB on the HPCx Power5 cluster, ANL IA64 cluster ‘Abe’ and ANL IBM BlueGene/P cluster. The system tested is a three-dimensional arterial neurovasculature with 7.7×106 fluid sites (filled circles, HPCx and FS; open circles, HPCx and FS+VR; crosses, Abe and FS; filled squares, Abe and FS+VR; asterisks, Blue Gene/P and FS; open squares, BlueGene/P and FS+VR). (b) Scalability on the LONI network using MPIg with 32 cores per site (dashed curve with circles, FS; dashed curve, ideal and FS; dashed curve with filled squares, FS+VR; dashed curve with open squares, FS+IS). The sites refer to three geographically separate IBM P5 clusters. FS+VR does not scale as well due to the global communications that are needed for rendered frame compositions. However, in reality, VR has little impact on the overall speed since it is only conducted at a rate of 25 times s−1.

In order to avoid post-processing of data, we employ in situ visualization, an optimal strategy for large-scale visualization where rendering is done within the same cores that perform the simulation in a unified and tightly coupled approach. Within HemeLB, each subdomain is rendered independently, corresponding to the LB domains. Owing to the requirement of global MPI communications to compose the final image, the rendering aspect of HemeLB at each time step is not as efficient as the LB code alone (figure 7). Computational steering is used to control and monitor physical parameters of the simulation along with visualization parameters (viewpoint, colour maps, etc.), in real time. The code is structured in such a way that the performance of HemeLB is decoupled from the network conditions, such that scenes are rendered only as frames are sent. Thus, the impact of rendering on the overall speed of the simulation is minimal. Depending on the performance of HemeLB on the particular machine as well as the network performance, rendering takes place only every 10–200 LB time steps. As a result, we can volume render in real time, at 10242 pixel resolution and 25 frames s−1 while the LB simulation itself runs at near full speed, between our local desktop visualization client (figure 8) and globally distributed remote resources.

Figure 8
Figure 8

Real-time interactive rendering at 10242 resolution of patient-specific cerebral haemodynamics. (a) Volume rendering of the blood velocity, (b) volume rendering of the internal blood von Mises stress, (c) pressure against the arterial wall and (d) stress against the arterial wall.

The software environment developed for use with HemeLB aims to hide unnecessary grid-computing details from clinicians, while bringing to the forefront the details and processes that they need to be aware of. This is particularly important given the time scales involved in the clinical decision-making process in interventional neuroradiology. During a surgical procedure, three-dimensional rotational angiography is used to rapidly create a volume highlighting the vasculature (6003 voxels, 440 μm resolution). Flow data are needed within 15–20 min of this acquisition. Within this time, data are segmented into the neurovasculature and boundary conditions are applied, the patient-specific data are staged using the AHE and the job is then executed within a reservation, as a cross-site or single-site run. Finally, as the job commences running, the client connects to the running simulation, which is used for steering and visualization purposes. The availability of an urgent computing system, or pre-reserved slots on remote petascale resources, is imperative. The scalability properties of HemeLB also significantly contribute to the rapid turnaround required. For example, using 1024 cores on the NCSA (http://teragrid.ncsa.uiuc.edu) Abe cluster, a single pulsatile cycle of an entire arterial neurovasculature (10243 voxels, 0.245 mm resolution, 106 fluid sites) takes approximately 9 min, well within the clinical requirement.

To the best of our knowledge, this is the first method developed with the ability to simulate whole-brain haemodynamics. In collaboration with neuroradiologists at the National Hospital for Neurology and Neurosurgery, London, UK, we are currently carrying out simulations on real medical cases. The longer-term goal is to determine the direct clinical utility of not only cerebral haemodynamics, but also real-time patient-specific medical-simulation methods in general, where validation, verification and clinical studies will be required.

7. Conclusion

In this paper, we describe five novel scientific applications that have been enabled by the availability of petascale resources. We have observed excellent scalability of our codes on up to tens of thousands of cores without any specific optimizations for individual machine architectures. This computational power has allowed us to perform unprecedented simulations in order to investigate complex phenomena over a range of applications, concurrently. One of the major operational issues associated with running such ‘capability’ jobs on large core-count machines is the low mean time to failure compared with the simulation run-time. Future growth in petascale resource usage within our work will involve integrating fault tolerance (FT) into scientific application codes, for example by using fault-tolerant MPI (FT-MPI; Chen et al. 2005), so that they can survive multiple bursts of process failures with low overhead.

Acknowledgments

We are grateful for funding from EPSRC grants EP/C536452/1, EP/E045111/1 and EP/F00561X/1, the Virolab project (EU IST STREP Project 027446) and the US NSF MRAC grants TG-DMR070013N and TG-DMR070014N. We acknowledge a recently awarded DOE INCITE grant entitled ‘Large scale condensed matter and fluid dynamics simulations’ at ALCF, which is supported by the US Department of Energy under contract DE-AC02-06CH11357. We are grateful to John Boisseau and the TACC support team for facilitating access to Ranger and to Peter Beckmann and Ramesh Balakrishnan for facilitating access to the Surveyor and Intrepid systems at ALCF.

Footnotes

  • One contribution of 16 to a Theme Issue ‘Crossing boundaries: computational science, e-Science and global e-Infrastructure I. Selected papers from the UK e-Science All Hands Meeting 2008’.

References

  1. ↵
    1. Artuso R.,
    2. Aurell E.,
    3. Cvitanović P.
    1990 Recycling of strange sets: I. Cycle expansions. Nonlinearity. 3, 325–359. doi:10.1088/0951-7715/3/2/005.
    OpenUrlCrossRefWeb of Science
  2. ↵
    1. Berghausen J.,
    2. Zipfel J.,
    3. Lindner P.,
    4. Richtering W.
    1998 Shear-induced orientations in a lyotropic defective lamellar phase. Europhys. Lett. 43, 683–689. doi:10.1209/epl/i1998-00417-3.
    OpenUrlCrossRefWeb of Science
  3. ↵
    Boghosian, B. M., Coveney, P. V., Fazendeiro, L., Latt, J., Tam, J. & Tang, H. 2008 New variational principles for locating periodic orbits of differential equations. Preprint.
  4. ↵
    1. Brooks B.R.,
    2. Janezic D.,
    3. Karplus M.
    1995 Harmonic analysis of large systems. I. Methodology. J. Comput. Chem. 16, 1522–1542. doi:10.1002/jcc.540161209.
    OpenUrlCrossRefWeb of Science
  5. ↵
    1. Chen Z.-R.,
    2. Kornfield J.A.,
    3. Smith S.D.,
    4. Grothaus J.T.,
    5. Satkowski M.M.
    1997 Pathways to macroscale order in nanostructured block copolymers. Science. 277, 1248–1253. doi:10.1126/science.277.5330.1248.
    OpenUrlAbstract/FREE Full Text
  6. ↵
    Chen, Z., Fagg, G. E., Gabriel, E., Langou, J., Angskun, T., Bosilca, G. & Dongarra, J. 2005 Fault tolerant high performance computing by a coding approach. In PPoPP '05: Proc. 10th ACM SIGPLAN Symp. on Principles and Practice of Parallel Programming, Chicago, IL, USA, pp. 213–223. New York, NY: ACM. (doi:10.1145/1065944.1065973)
  7. ↵
    1. Chen B.,
    2. Evans J.R.G.,
    3. Greenwell H.C.,
    4. Boulet P.,
    5. Coveney P.V.,
    6. Bowden A.A.,
    7. Whiting A.
    2008 A critical appraisal of polymer–clay nanocomposites. Chem. Soc. Rev. 37, 568–594. doi:10.1039/b702653f.
    OpenUrlCrossRefPubMed
  8. ↵
    1. Coveney P.V.,
    2. Saksena R.S.,
    3. Zasada S.J.,
    4. McKeown M.,
    5. Pickles S.
    2007 The Application Hosting Environment: lightweight middleware for grid-based computational science. Comput. Phys. Commun. 176, 406–418. doi:10.1016/j.cpc.2006.11.011.
    OpenUrlCrossRefWeb of Science
  9. ↵
    Coveney, P. V. et al. In press. Large scale computational science on federated international grids: the role of switched optical networks. Future Gen. Comput. Syst. (doi:10.1016/j.future.2008.09.013)
  10. ↵
    Cvitanović, P., Artuso, R., Mainieri, R., Tanner, G., Vattay, G., Whelan, N. & Wirzba, A. 2004 Chaos: classical and quantum. See http://www.chaosbook.org/.
  11. ↵
    1. Dombre T.,
    2. Frisch U.,
    3. Greene J.M.,
    4. Hénon M.,
    5. Mehr A.,
    6. Soward A.M.
    1986 Chaotic streamlines in the ABC flows. J. Fluid Mech. 167, 353–391. doi:10.1017/S0022112086002859.
    OpenUrlCrossRef
  12. ↵
    1. Ermolieff J.,
    2. Lin X.,
    3. Tang J.
    1997 Kinetic properties of saquinavir-resistant mutants of human immunodeficiency virus type 1 protease and their implications in drug resistance in vivo. Biochemistry. 36, 12 364–12 370. doi:10.1021/bi971072e.
    OpenUrl
  13. ↵
    1. Fujitani H.,
    2. Tanida Y.,
    3. Ito M.,
    4. Jayachandran G.,
    5. Snow C.D.,
    6. Shirts M.R.,
    7. Sorin E.J.,
    8. Pande V.S.
    2005 Direct calculation of the binding free energies of FKBP ligands. J. Chem. Phys. 123, 084 108 doi:10.1063/1.1999637.
    OpenUrl
  14. ↵
    1. Giupponi G.,
    2. Harting J.,
    3. Coveney P.V.
    2006 Emergence of rheological properties in lattice Boltzmann simulations of gyroid mesophases. Europhys. Lett. 73, 533–539. doi:10.1209/epl/i2005-10438-x.
    OpenUrlCrossRefWeb of Science
  15. ↵
    1. Jayachandran G.,
    2. Shirts M.R.,
    3. Park S.,
    4. Pande V.S.
    2006 Parallelized-over-parts computation of absolute binding free energy with docking and molecular dynamics. J. Chem. Phys. 125, 084901 doi:10.1063/1.2221680.
    OpenUrlCrossRefPubMed
  16. ↵
    1. Karonis N.T.,
    2. Toonen B.,
    3. Foster I.
    2003 MPICH-G2: a grid-enabled implementation of the Message Passing Interface. J. Parallel Distrib. Comput. 63, 551–563. doi:10.1016/S0743-7315(03)00002-9.
    OpenUrlCrossRef
  17. ↵
    1. Love P.J.,
    2. Nekovee M.,
    3. Coveney P.V.,
    4. Chin J.,
    5. González-Segredo N.,
    6. Martin J.M.R.
    2003 Simulations of amphiphilic fluids using mesoscale lattice–Boltzmann and lattice–gas methods. Comput. Phys. Commun. 153, 340–358. doi:10.1016/S0010-4655(03)00200-5.
    OpenUrlCrossRefWeb of Science
  18. ↵
    1. Maschera B.,
    2. Darby G.,
    3. Palú G.,
    4. Wright L.L.,
    5. Tisdale M.,
    6. Myers R.,
    7. Blair E.D.,
    8. Furfine E.S.
    1996 Human immunodeficiency virus. Mutations in the viral protease that confer resistance to saquinavir increase the dissociation rate constant of the protease–saquinavir complex. J. Biol. Chem. 271, 33 231–33 235. doi:10.1074/jbc.271.52.33231.
    OpenUrl
  19. ↵
    1. Mazzeo M.D.,
    2. Coveney P.V.
    2008 HemeLB: a high performance parallel lattice–Boltzmann code for large scale fluid flow in complex geometries. Comput. Phys. Commun. 178, 894–914. doi:10.1016/j.cpc.2008.02.013.
    OpenUrlCrossRefWeb of Science
  20. ↵
    1. Pinnavaiam T.J.,
    2. Beall G.W.
    Polymer–clay nanocomposites. In Wiley 2000 Chichester, UK:Wiley
  21. ↵
    1. Prabu-Jeyabalan M.,
    2. Nalivaika E.,
    3. Schiffer C.A.
    2002 Substrate shape determines specificity of recognition for HIV-1 protease: analysis of crystal structures of six substrate complexes. Structure. 10, 369–381. doi:10.1016/S0969-2126(02)00720-7.
    OpenUrlCrossRefPubMed
  22. ↵
    1. Rey A.D.,
    2. Denn M.M.
    2002 Dynamical phenomena in liquid-crystalline materials. Annu. Rev. Fluid Mech. 34, 233–266. doi:10.1146/annurev.fluid.34.082401.191847.
    OpenUrlCrossRefWeb of Science
  23. ↵
    1. Sadiq S.K.,
    2. Wright D.,
    3. Watson S.J.,
    4. Zasada S.J.,
    5. Stoica I.,
    6. Coveney P.V.
    Automated molecular simulation based binding affinity calculator for ligand-bound HIV-1 proteases. J. Chem. Inf. Model. 48, 2008a 1909–1919. doi:10.1021/ci8000937.
    OpenUrlCrossRefPubMedWeb of Science
  24. ↵
    1. Sadiq S.K.,
    2. et al.
    Patient-specific simulation as a basis for clinical decision-making. Phil. Trans. R. Soc. A. 366, 2008b 3199–3219. doi:10.1098/rsta.2008.0100.
    OpenUrlAbstract/FREE Full Text
  25. ↵
    1. Saksena R.S.,
    2. Coveney P.V.
    Self-assembly of ternary cubic, hexagonal, and lamellar mesophases using the lattice–Boltzmann kinetic method. J. Phys. Chem. B. 112, 2008a 2950–2957. doi:10.1021/jp0731506.
    OpenUrlPubMed
  26. ↵
    Saksena, R. S. & Coveney, P. V. 2008b Massively parallel lattice–Boltzmann simulations for the investigation of cubic phase rheology. In Proc. TeraGrid '08 Conf., 9–13 June.
  27. ↵
    Schwarz, U. S. & Gompper, G. 2002 Morphology of condensed matter: physics and geometry of spatially complex systems. Lecture Notes in Physics, vol. 600, pp. 107–151. Berlin, Germany: Springer.
  28. ↵
    1. Seddon J.M.,
    2. Squires A.M.,
    3. Conn C.E.,
    4. Ces O.,
    5. Heron A.J.,
    6. Mulet X.,
    7. Shearman G.C.,
    8. Templer R.H.
    2006 Pressure-jump X-ray studies of liquid crystal transitions in lipids. Phil. Trans. R. Soc. A. 364, 2635–2655. doi:10.1098/rsta.2006.1844.
    OpenUrlAbstract/FREE Full Text
  29. ↵
    1. Shirts M.,
    2. Pande V.S.
    2000 Screen savers of the world unite!. Science. 290, 1903–1904. doi:10.1126/science.290.5498.1903.
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Sloot P.,
    2. Coveney P.,
    3. Bubak M.T.,
    4. Vandamme A.-M.,
    5. Ó Nualláin B.,
    6. van de Vijver D.,
    7. Boucher C.
    Multi-science decision support for HIV drug resistance treatment. In Studies in Health Technology and Informatics Global Healthgrid: e-Science meets biomedical informatics. vol. 138 2008 Amsterdam, The Netherlands:IOS Press
  31. ↵
    1. Squires A.M.,
    2. Templer R.H.,
    3. Seddon J.M.,
    4. Woenckhaus J.,
    5. Winter R.,
    6. Finet S.,
    7. Theyencheri N.
    2002 Kinetics and mechanism of the lamellar to gyroid inverse bicontinuous cubic phase transition. Langmuir. 18, 7384–7392. doi:10.1021/la0259555.
    OpenUrlCrossRefWeb of Science
  32. ↵
    1. Stoica I.,
    2. Sadiq S.K.,
    3. Coveney P.V.
    2008 Rapid and accurate prediction of binding free energies for saquinavir-bound HIV-1 proteases. J. Am. Chem. Soc. 130, 2639–2648. doi:10.1021/ja0779250.
    OpenUrlCrossRefPubMedWeb of Science
  33. ↵
    1. Succi S.
    The lattice Boltzmann equation for fluid dynamics and beyond. In Oxford University Press 2001 Oxford, UK:Oxford University Press
  34. ↵
    Suter, J. L. & Coveney, P. V. In press. Computer simulation study of the materials properties of intercalated and exfoliated poly(ethylene)glycol clay nanocomposites. Soft Matter (doi:10.1039/b822666k)
  35. ↵
    1. Temam R.
    Navier–Stokes equations. In North-Holland 1984 Amsterdam, The Netherlands:North-Holland
  36. ↵
    1. Wang W.,
    2. Kollman P.A.
    2001 Computational study of protein specificity: the molecular basis of HIV-1 protease drug resistance. Proc. Natl Acad. Sci. USA. 98, 14 937–14 942. doi:10.1073/pnas.251265598.
    OpenUrl
  37. ↵
    1. Wang W.,
    2. Donnii O.,
    3. Reyes C.M.,
    4. Kollman P.A.
    2001 Biomolecular simulations: recent developments in force fields, simulations of enzyme catalysis and protein–ligand, protein–protein and protein–nucleic acid non-covalent interactions. Annu. Rev. Biophys. Biomol. Struct. 30, 211–243. doi:10.1146/annurev.biophys.30.1.211.
    OpenUrlCrossRefPubMedWeb of Science
View Abstract
Back to top
Email
Print
Manage alerts
Citation tools
Download

Article reuse

Get Permission

See related subject areas:

Related articles

Cited by

Celebrating 350 years of Philosophical Transactions

Anniversary issue with free commentaries, archive material, videos and blogs.