Keynote Addresses

Putting the hardware together to form a TeraScale simulation environment is only a small part of the challenge. Building the software environment to enable applications to run efficiently on such a system is quite another challenge. This talk will outline computer science research under investigation at the Advanced Computing Laboratory at Los Alamos in support of application development for TeraScale platforms. The technologies which will be discussed include scalable run-time systems, object-oriented frameworks, clustering technology, and scalable visualization. In particular, discussion of how these technologies must work together in an integrated fashion to enable requisite performance on clustered shared memory processor systems will be examined. Mixed-model approaches to scalable simulation outside the normal message-passing approach and techniques for the encapsulation of concurrency will be a primary focus of this presentation. These technologies will be described in the context of how they enabled actual applications which scale across thousands of processors in multiple programs at the Advanced Computing Laboratory and how they are aligned with the current procurement strategy of choice for high-end computing in the United States Department of Energy.

We present an overview of some major issues in computational simulation of coupled, multicomponent systems, as exemplified by the simulation of solid propellant rockets. High-fidelity simulation of such a complex system requires simultaneous simulation of fluid, solid, and combustion components and their interactions with each other. Topics addressed include algorithms for data transfer between components with diverse discretizations, consistent time stepping, software integration, and scalability on parallel computers.

A Computational Enzymology, which combines quantum and molecular mechanics (QM/MM) methods, is used to simulate the complete minimum energy surface and reaction pathway for the interconversion of malate and oxaloacetate catalyzed by the enzyme malate dehydrogenase (MDH). A reaction mechanism for the proton and hydride transfers associated with MDH is deduced from the topology of the calculated energy surface. The putative mechanism consists of (1) a sequential reaction with the proton transfer preceding the hydride transfer (malate to oxaloacetate direction), (2) the existence of two transition states with energy barriers of about 10 and 15 kcal/mol for the proton and hydride transfers, respectively, and (3) reactant (malate) and product (oxaloacetate) states that are nearly isoenergetic. A simulation analysis of the calculated energy profile shows that solvent effects due to the protein matrix dramatically alter the intrinsic reactivity of the functional groups involved in the MDH reaction. The enzyme effectively changes the reaction

imidazole + malate + NAD imidazolium + oxaloacetate + NADH

from an exothermic reaction () in the gas phase to a nearly isoenergetic one () in the protein_solvent environment of MDH. An energy decomposition analysis indicates that specific MDH residues (Asp-150, Arg-81, Arg-87, and Arg-153) in the vicinity of the substrate make significant energetic contributions to the stabilization of the proton transfer and destabilization of the hydride transfer. This suggests that these amino acids play an important role in the catalytic properties of MDH, which is consistent with site-directed mutagenesis experiments.

Fellows' Talks

For the last four years, I have been researching problems in computational ecology. Specifically, I have considered the application of parallel discrete event simulation to problems in demography. The consideration of discrete event simulation stems from two recent trends in ecological modeling, namely individual oriented modeling and spatially realistic modeling. In an individual oriented model, we model a population by modeling the physiological and social interactions of its members. Recent work has shown that spatial distribution of organisms tends to explain some demographic stability. This has led some ecological modelers including those interested in individual oriented models--to focus on the effect of spatial realism.

With few exceptions, modelers implement spatially explicit, individual oriented models as discrete time models. Discrete time models assume we can represent interactions by updating a model at regular intervals, but ecosystems are not generally synchronous when considering system behavior at the individual level. As a consequence, discrete event models may be more realistic and more efficient to implement.

In the initial phase of my research, I reviewed parallel simulation techniques and implemented a spatially explicit model using the Time Warp Operating System (Time Warp is a major parallel simulation algorithm). This experience led me to two conclusions: the performance of parallel simulation algorithms are not characterized well enough to comment on their appropriateness for a given problem and ecological models must be developed as discrete event simulations before we attempt to fit ecological models to a particular parallel algorithm.

The goal of my research is two fold. First I am working, in conjunction with ecologists, to develop discrete event ecological models. The primary product of this research is the Ecological Simulation System (ESS) a framework for population modeling based on Zeigler's Theory of Modeling and Simulation. The second goal of my research is a performance analysis of realistic ecological simulation, implemented with parallel simulation algorithms. The end of this research is to develop a general simulation tool for ecological modeling.

The ocean is a stratified medium. As a result, waves can propagate horizontally at density interfaces or at an angle to the horizontal within continuously stratified regions. A significant portion of internal wave energy along coastal regions is generated by the interaction of the surface tides with the bottom topography. The surface, or barotropic, tide perturbs the density profile as different density fluids are forced over topography at the edges of the continental shelves. These perturbations generate internal waves at tidal frequency, or internal tides, that propagate both out to the open ocean as well as on to the continental shelf. Under the right conditions, internal waves can focus enough energy at a point such that either a shear or convective instability can form and cause them to break. This breaking process generates turbulence and mixing that can cause significant sediment transport and nutrient redistribution and can modify acoustic wave propagation and optical clarity within the littoral ocean. Furthermore, breaking internal waves are thought to account for a significant portion of the dissipation within the ocean, concentrating most of it at coastal boundaries. This project focuses on a large eddy simulation of the internal wave breaking process at an interface between two fluids of different densities. In order to create internal breaking waves, the instability can be forced by either the use of characteristic focusing or by wave steepening due to topographic effects. Characteristic focusing relies on the dispersive nature of internal waves. By sending out short, slow waves followed by long, fast waves, an internal wave chirp signal can be generated in which the fast waves overtake the slow waves and focus to a point at which they generate a breaking instability. On the other hand, sending an internal wave train over sloped topography causes unstable growth in wave amplitude that leads to a breaking instability as well.

The present numerical study is performed in conjunction with a laboratory scale study in order to verify the results. With these results, the goal is to parameterize the dissipative mechanisms of internal wave breaking in order to allow a full-scale simulation of internal wave breaking in Monterey Bay on a Beowulf cluster of Pentium III workstations. The model we use incorporates a free surface as well as a nonhydrostatic pressure solver. Typically, numerical simulations of coastal processes do not solve for the nonhydrostatic pressure as the hydrostatic approximation holds for most cases. In Monterey Bay, however, measurements have shown that the internal-wave energy cascade transfers enough energy to the nonhydrostatic regime to make a nonhydrostatic model a necessity. This increases the required resolution to roughly 8 million grid points for a reasonable simulation of Monterey Bay, making parallel computation a necessity as well.

Nowadays, computers are used in many areas traditionally dependent on human expertise, such as medicine. Computer-assisted orthopeadic surgery (CAOS) is the new area of computer application in medicine. It consists of three components: modeling, planning and execution. One of the most common sources of preoperative information in CAOS is medical imaging. Advances in medical imaging technology (computed tomography (CT), magnetic resonance imaging (MRI), etc.) and in computer-based image processing and modeling have made it possible to geometrically model and

visualize anatomical structures in patients. These developments have given rise to the computer-based preoperative planners that are now appearing. The surgeon can interactively determine the operative procedure by manipulating computer-generated models of anatomical structures, surgical tools and/or implants. The computer can perform the simulations of the operative plan, giving the surgeon feedback on the proposed surgical procedure.

Existing preoperative planners do not relate the operation plan with the biomechanics of anatomical structures. Understanding the mechanics of orthopaedic surgery is particularly important, especially for cementless press-fitted implants. This work is the first attempt to incorporate patient-specific biomechanical simulation into a geometrical preoperative surgical planner. This can provide the surgeon with qualitatively new tools for computer-assisted surgery.

The complexity of anatomical structures requires the use of numerical methods, such as the FE method, to perform biomechanical simulations. Existing mesh generation techniques cannot automatically produce FE meshes of three-dimensional models of arbitrary geometry. Moreover, geometric modeling of patient-specific anatomical structure, which underlies existing preoperative planners, is still under development. A new method to automatically obtain three-dimensional FE models from CT scan data will be developed. This method utilizes mapping techniques, commonly applied for registration of MRI of the human brain. The method provides segmentation as a byproduct, which can be compared to the geometric planner's segmentation as a validation step. In the future, surface models obtained from FE models can also be used within geometrical preoperative planners. FE simulation of the acetabular component of Total Hip Replacement (THR) will be performed using patient-specific three-dimensional pelvic model of anatomical geometry.

Although this research is motivated by a specific problem of biomechanical analysis of the acetabular component of THR, the methods to be developed within this framework are general and can be applied to other biomechanical and medical imaging problems.

Classical Geostatistics methods have been designed to use mainly statistical knowledge about natural variables and they lack the ability to incorporate important forms of knowledge like physical laws and scientific theories into the mapping process. On the other hand, the powerful and versatile Bayesian maximum entropy (BME) method of Modern Geostatistics can accomplish such a task, rigorously and efficiently. In this work, BME is used to incorporate the Darcy law of subsurface hydrology in the spatial mapping of a hydraulic head field. The hydraulic map thus obtained is physically meaningful as well as numerically more accurate than that obtained using classical methods (e.g., kriging). Moreover, taking advantage of the Darcy law the BME hydraulic head mapping can involve other related soil properties, like hydraulic conductivity. The approach leads to very accurate hydraulic head solutions, and may be also applied to study the inverse problem, in which one seeks to estimate hydraulic conductivity from hydraulic head measurements.

The Element-Free Galerkin method has already been demonstrated to accurately model dynamic fracture problems. However, it did not have an obvious means of modeling contact across crack faces. The means of modeling contact will be briefly outlined, and results shown.

The Partition of Unity Finite Element Method (PU-FEM) has recently come to light, and been applied to static crack problems. It is a consistent means of extrinsically enriching a sub-set of finite elements in the model with linearly independent functions. When these functions are chosen to be a set which can construct the near-tip asymptotic crack fields, then an arbitrary crack can be accurately modeled, by even a relatively coarse mesh.

Since PU-FEM seems to be a robust, accurate method in statics, I have attempted to apply it to dynamic problems. The formulation of the problem has to be dealt with carefully, since a propagating crack causes the shape functions themselves to be time-dependent. The formulation will be explained, and results shown.

Modeling cracks and crack growth with the finite element method is difficult due to the need for both the construction and regeneration of a fine mesh near the crack tip. We present a technique to enrich a standard displacement-based finite element approximation near a crack by incorporating a discontinuous field and the near tip asymptotic fields with a partition of unity method. In comparison to modeling the crack geometrically with the mesh, the crack is modeled physically with the special enrichment functions. This technique provides for both coarse mesh accuracy and arbitrary crack alignment within the mesh. Crack growth is modeled in a straightforward fashion by redefining the physical extent of the discontinuity: no remeshing is necessary. The advantages of the present technique over meshless and boundary element methods are discussed. Several example problems of fracture in two dimensional elasticity and Mindlin-Reissner plates are provided to illustrate the accuracy and utility of the new formulation.

Recently Wavelet Transform based compression techniques have demonstrated, superior performance over Fourier (and other traditional) Transform based compression techniques in terms of high quality, low bit-rate compression, for many types of signals. Even more recently it has be shown that Adaptive Wavelet Packets (variant of the Wavelet Transform) based techniques have the promise of greatly outperforming the standard Wavelet Transform techniques. A potential hindrance to the use of the Adaptive Wavelet Packet techniques, in very low bit-rate applications, is the amount of overhead information needed to encode what is the "best" wavelet packet tree structure for processing each signal segment (block) in a sequence of many signal segments (blocks) to be processed. For example, in the compression of 2-D image blocks, if 4-level Wavelet Packet Tree Structures are utilized, then, 66 bits of overhead information are needed. For 5-level tree structures, 261 overhead bits are needed. The number of overhead bits (along with the versatility of the compression algorithm) increases greatly with depth of the tree structures considered. We have developed an Adaptive Wavelet Packet method that avoids sending this overhead information all together. We call this the Predictive Adaptive Wavelet Packet method. Initial results from this research indicate that the new Predictive Adaptive Wavelet Packet (PAWP) method can achieve a rate/distortion performance that is better than that of regular Adaptive Wavelet Packet (AWP) methods, in which overhead information is needed. The trade-off is that about twice as much computation is required for PAWP as for AWP. Results have been acquired for speech signals and image sequences.

Rarely in science is one interested in quantities that can be measured directly. More often, the interesting parts of a system are only accessible by means of a model relating them to (experimentally measurable) external features. Furthermore, since the number of measurements that can be taken is necessarily finite, the problem of deciding which ones to take must inevitably be addressed. This is the optimal sampling problem.

X-ray tomography provides a good example of this. Suppose an image of a certain area inside the body is to be reconstructed from a fixed number of views. From which angles should the views be collected so as to permit the best reconstruction of the region of interest?

We will formulate the optimal sampling problem in a very general setting, consider what is meant by the "best" sampling scheme, and look at some examples.

Droplet splashing occurs in a wide variety of applications, including spray painting processes, injection of liquid fuel in internal combustion engines, fire suppression systems, and many others. The subject is definitely not new and has received attention by numerous researchers for over a century. Due to its complexity, however, characterized by large deformations over very small time scales, propagations of instabilities in the emerging liquid crown and sheet, formation of jets along the rim, and subsequent jet break, the study of the splashing phenomena is far from complete.

Our interests in this area are the determination of the velocity and size characteristics of secondary droplets created during splashing, as well as the calculation of quantities that described the overall splashing event (e.g. the amount of energy lost due to viscous dissipation). With an emphasis towards automotive applications, the single droplet splashing study presented here provides a foundation for the investigation of spray impingement on containment walls, where thousands of droplets collide against a solid surface.

Our approach is based on describing the products of splashing in terms of a joint probability density function (jpdf) of an assumed form, and solving this jpdf by applying the conservation laws to the splashing process. Drop size calculations from this model are compared to experimental drop size data yielding acceptable agreement. The distribution of incident droplet energy among the kinetic and surface energy of the secondary droplets, the residual film surface energy, and the time integrated viscous dissipation are calculated as a function of the incident Weber number and the impingement angle. Information is also obtained on the incident angle dependence of the impact forces.

A fundamental problem for bootstrap current-driven magnetic islands is to understand the dynamics and profile properties of density and electron temperature. Although electron temperature flattening inside the island separatrix has been observed, it is uncertain the degree to which density profiles show a similar flattening. To study the evolution of density and electron temperature in a helical island geometry, we solve the coupled parallel momentum and continuity equations in the presence of a slower-growing magnetic perturbation that simulates the evolution of a neoclassical tearing mode. It is assumed that sound wave propagation along field lines is primarily responsible for equilibration, of density over perturbed flux surfaces. Electron temperature equilibration, on the other hand, occurs on perturbed flux surfaces as a result of a rapid parallel heat flux. Properly accounting for this heat flux requires a solution of an appropriate kinetic equation that takes into account free-streaming and collisional effects in a helical magnetic geometry. An analytic closure for the heat flux is constructed based on a multiple time and spatial scale scale, Chapman-Enskog-like analysis. This heat flux is then inserted into a temperature evolution equation which is solved in the presence of an evolving helical magnetic island.

Many nonequilibrium problems in physics and engineering have natural mathematical formulations as hyperbolic systems with relaxation source terms. Such systems pose interesting challenges for numerical approximation. The source terms cause the system to be dispersive; generally, both the eigenvalues and eigenvectors of the system are modified by the relaxation processes which drive the system towards equilibrium. Another difficult numerical issue arises when the relaxation source terms operate on much smaller time scales than the advection terms; the problem is then said to be stiff. Unfortunately, many problems simultaneously exhibit behavior in both the stiff and non-stiff limits, as well as in between. It is therefore desirable to develop high-resolution numerical algorithms which are uniformly accurate irrespective of the scales of the relaxation source terms. For example, if the data are such that the advection terms permit a time step large enough to bring the flow into local equilibrium, one would wish to see an accurate solution of the equilibrium problem emerge without resolving the details of the transition.

In recent years, there has been much development on such uniformly accurate methods. A variety of splitting schemes, characteristic schemes, and staggered, centered schemes have been proposed. The most promising of these rely on staggered grids to eliminate the difficulty of resolving the Riemann problem, which no longer has a self-similar solution. However, one drawback of centered schemes is that one may lose definition of contact discontinuities, which are likely to be very important in reactive flows. The question still remains as to whether there is a simple, natural extension of the upwind approach which includes the effects of the relaxation source terms and resolves all discontinuities equally well. The close coupling between the advection and relaxation has lead us to explore an upwind scheme based upon an approximate solution to the Riemann initial value problem. Currently, experiments are being performed which apply this method to a model linear system as well as to equations of extended hydrodynamics.

The equations of resistive magnetohydrodynamics pose many challenges to computational science. In particular, the presence of widely separated fast and slow scales imposes severe Courant-Friedrichs- Lewy condition on the time step and consequently makes computing long term evolution of low-beta plasmas difficult. It is the long term dynamics, however, that we are interested in the study of solar plasmas and laboratory simulations of prominence eruptions. In this talk we will discuss how new exponential propagation methods can be used to overcome these difficulties in numerical modeling of spheromak plasma experiments and present example simulations which demonstrate this.

Molecular dynamics programs simulate the behavior of biomolecular systems, leading to insights and understanding of their functions. However, the computational complexity of such simulations is enormous. Parallel machines provide the potential to meet this computational challenge. To harness this potential, it is necessary to develop a scalable program. It is also necessary that the program be easily modified by application-domain programmers. The program NAMD2 seeks to provide these desirable features. It uses spatial decomposition combined with force decomposition to enhance scalability. It uses intelligent periodic load balancing, so as to maximally utilize the available compute power. It is modularly organized, and implemented using a parallel C++ dialect, so as to enhance its modifiability. It uses a combination of numerical techniques and algorithms to ensure that energy drifts are minimized, ensuring accuracy in long running calculations. NAMD2 uses a portable run-time framework that also supports interoperability among multiple parallel paradigms. As a result, different components of applications can be written in the most appropriate parallel paradigms. NAMD2 runs on most parallel machines including workstation clusters.

Lithium transition metal oxides (Li_x MO₂ with M = transition metal) are ideally suited to serve as cathode materials in the next generation of rechargeable lithium batteries. A lithium battery consists of an anode (which is metallic lithium) and a cathode (the lithium transition metal oxide) separated by an electrolyte. In use, lithium ions are released at the anode, travel through the electrolyte and are intercalated into the cathode material. Lithium transition metal oxides exhibit a high voltage, can undergo large variations in lithium concentration and are stable under repeated removal and insertion of lithium ions. Currently the most important cathode material is Li_x CoO₂ . In this work we have investigated the thermodynamic phase stability of Li_x CoO₂ as a function of the lithium concentration. We have also calculated the mobility of the lithium ions within the CoO₂ crystal structure. Both quantities are essential to enable an intelligent optimization and design of existing and new cathode materials.

The organometallic vapor phase epitaxy (OMVPE) process is widely used to produce compound semiconductors for applications ranging from solar cells to cellular phones. The OMVPE process is a chemical vapor deposition (CVD) process where thin films are synthesized from organometallic precursors. As device dimensions decrease, processing issues that affect the material properties (i.e. stringent demands on film thickness, morphology and impurity incorporation) become increasingly important. Models of OMVPE can be used in fundamental studies of the processing of thin films, as well as in predictive simulations to aid in process design and control. The difficulties in modeling the process arise from the complex nature of OMVPE; processes occur on widely different length and time scales. In this work, models on various length and time scales for the OMVPE process will be developed and methodologies for linking these models will be formulated. The "linked" models will be shown to make accurate predictions on various length scales. The complex coupling of gas-phase transport phenomena with gas-phase and surface chemical kinetics means that gaining an understanding of the entire OMVPE process requires more than one type of modeling approach. Predictions on the macroscopic length scale (growth rate, film uniformity, film composition) and on the "mesoscopic" length scale (surface morphology) are both needed to determine optimal operating conditions for OMVPE reactors. Models have been developed to make predictions on each of the relevant length scales:

• Macroscopic reactor scale models have been developed to understand the fluid flow, heat and mass transfer combined with detailed gas-phase chemical kinetic mechanisms to predict the type and concentration of precursors arriving at the growth front. Finite Element simulations have been developed in our group to solve these problems.
• "Mesoscopic" surface models have been developed to study the evolution of surface morphology during growth. Kinetic Monte Carlo methods have been used to investigate the importance of individual steps in surface mechanisms, as well as determine the growth mode on the surface. Parallel algorithms are also in the process of being developed in order to speed the computations of the surface problem.

A methodology is developed to "link" reactor scale and surface models in order to gain an overall understanding of the OMVPE process. Case studies of GaAs growth using a variety of precursors is used to validate the models separately, as well as to validate the "linked" model. Using the reactor scale model alone cannot give predictions of surface morphology. The surface model used alone is missing a needed input of gas-phase flux to the surface. The "linked" model is able to make extremely accurate predictions on both the macroscopic length scale (growth rate) and mesoscopic length scale (surface morphology).

The steady advance of both algorithmic and computational power in the last decade is leading us to an age when Molecular Simulation will be able to routinely, and quantitatively predict the fluid phase diagrams of industrially important molecules. This talk will briefly discuss the current state of the art of Configurational-bias Gibbs Ensemble Monte Carlo simulation, and present recent applications of this technique to the single and multi-component vapor-liquid phase coexistence of linear and branched alkanes at sub and super-critical temperatures and pressures. The partition constants of solutes between a squalene alkane liquid phase, and a Helium vapor phase are presented as an example of how we can compute the relative retention times of solutes in a chromatographic column.

Phosphine has significant potential as a precursor in chemical vapor deposition to form in situ phosphorus-doped epitaxial and polycrystalline silicon films. Due to this potential, we study phosphine's interaction with the Si(100)-2x1 surface. We characterize the potential energy surfaces of various mechanisms for adsorption and decomposition of phosphine using ab initio quantum chemical calculations. We find molecular adsorption occurs on the lower silicon atom in the tilted Si-Si dimer structure, with a binding energy of 16 kcal/mol. The molecular adsorption of phosphine is followed by dissociation to PH₂(a) + H(a), which occurs by three competing pathways. The most thermodynamically favored pathway involves hydrogen loss to the same dimer, which is exothermic by 48 kcal/mol relative to PH₃(g). The lowest energy barrier for hydrogen loss, however, involves hydrogen loss across the trench with an energy barrier 5 kcal/mol below the PH₃(g) state. The final pathway involves hydrogen loss to a neighboring dimer in the same row, whose energy barrier and thermodynamic stability lie in between those of the other two pathways.

Fellows' Posters

We address a very simple question here: can a sparse and recurrently-connected network of (very simplified) neuron-like elements learn, in an unsupervised fashion, to generate a reliable temporal sequence of activity patterns in response to a weak externally-originated input? The structure of this problem arises in several contexts related to the general theory of learning in neural networks, the biology of neural computation, and the cognitive phenomenon of trace classical conditioning. An approximate parameterization relating the reliability, temporal length, and spatial organization of the generated patterns to the architectural parameters of model networks is developed (via numerical simulations and qualitative analysis), as are the implications of these relationships for the surrounding scientific contexts.

Rarely in science is one interested in quantities that can be measured directly. More often, the interesting parts of a system are only accessible by means of a model relating them to (experimentally measurable) external features. Furthermore, since the number of measurements that can be taken is necessarily finite, the problem of deciding which ones to take must inevitably be addressed. This is the optimal sampling problem.

X-ray tomography provides a good example of this. Suppose an image of a certain area inside the body is to be reconstructed from a fixed number of views. From which angles should the views be collected so as to permit the best reconstruction of the region of interest?

We will formulate the optimal sampling problem in a very general setting, consider what is meant by the "best" sampling scheme, and look at some examples.

Modeling cracks and crack growth with the finite element method is difficult due to the need for both the construction and regeneration of a fine mesh near the crack tip. We present a technique to enrich a standard displacement-based finite element approximation near a crack by incorporating a discontinuous field and the near tip asymptotic fields with a partition of unity method. In comparison to modeling the crack geometrically with the mesh, the crack is modeled physically with the special enrichment functions. This technique provides for both coarse mesh accuracy and arbitrary crack alignment within the mesh. Crack growth is modeled in a straightforward fashion by redefining the physical extent of the discontinuity: no remeshing is necessary. The advantages of the present technique over meshless and boundary element methods are discussed. Several example problems of fracture in two dimensional elasticity and Mindlin-Reissner plates are provided to illustrate the accuracy and utility of the new formulation.

Initial work is presented on the application of adaptive mesh refinement to the simulation of multiphase flow and transport in the subsurface. In this work, adaptive mesh refinement is combined with higher-order Godunov methods and a hybrid mixed finite element discretization.

Olefin metathesis is an important reaction for many chemical processes including polymerization. It has been shown that the Grubbs ruthenium-based catalyst is very efficient at olefin metathesis with high functional group tolerance. It has been so successful that we would quickly deplete the world of its supply of ruthenium if we were to use it for all its applications. Looking at related elements that have a higher natural abundance leads to iron. Iron should have similar electronic behavior which could lead to a successful iron based olefin metathesis catalyst.

The mechanism of the ruthenium-based Grubbs catalyst of the formula L₂X₂Ru=CHR has been examined. Experimentally it has been shown that L=PCy₃,PⁱPr₃, X=Cl, and R=CHPh₂ produce an effective catalyst. When X=Cl and R=H in our theoretical study, L was varied from a simple phosphine to a carbene and an sp² nitrogen donating ligand. Understanding gained from this mechanism has allowed for possible performance enhancement of the ruthenium-based catalyst. In addition, the iron-based analog is showing promise as the ligand effects are being understood and tuned to the needs of the iron-based reaction pathway.

The ocean is a stratified medium. As a result, waves can propagate horizontally at density interfaces or at an angle to the horizontal within continuously stratified regions. A significant portion of internal wave energy along coastal regions is generated by the interaction of the surface tides with the bottom topography. The surface, or barotropic, tide perturbs the density profile as different density fluids are forced over topography at the edges of the continental shelves. These perturbations generate internal waves at tidal frequency, or internal tides, that propagate both out to the open ocean as well as on to the continental shelf. Under the right conditions, internal waves can focus enough energy at a point such that either a shear or convective instability can form and cause them to break. This breaking process generates turbulence and mixing that can cause significant sediment transport and nutrient redistribution and can modify acoustic wave propagation and optical clarity within the littoral ocean. Furthermore, breaking internal waves are thought to account for a significant portion of the dissipation within the ocean, concentrating most of it at coastal boundaries. This project focuses on a large eddy simulation of the internal wave breaking process at an interface between two fluids of different densities. In order to create internal breaking waves, the instability can be forced by either the use of characteristic focusing or by wave steepening due to topographic effects. Characteristic focusing relies on the dispersive nature of internal waves. By sending out short, slow waves followed by long, fast waves, an internal wave chirp signal can be generated in which the fast waves overtake the slow waves and focus to a point at which they generate a breaking instability. On the other hand, sending an internal wave train over sloped topography causes unstable growth in wave amplitude that leads to a breaking instability as well.

The present numerical study is performed in conjunction with a laboratory scale study in order to verify the results. With these results, the goal is to parameterize the dissipative mechanisms of internal wave breaking in order to allow a full-scale simulation of internal wave breaking in Monterey Bay on a Beowulf cluster of Pentium III workstations. The model we use incorporates a free surface as well as a nonhydrostatic pressure solver. Typically, numerical simulations of coastal processes do not solve for the nonhydrostatic pressure as the hydrostatic approximation holds for most cases. In Monterey Bay, however, measurements have shown that the internal-wave energy cascade transfers enough energy to the nonhydrostatic regime to make a nonhydrostatic model a necessity. This increases the required resolution to roughly 8 million grid points for a reasonable simulation of Monterey Bay, making parallel computation a necessity as well.

Long-term patient support with oxygenators results in hematologic complications as blood elements are exposed to unnatural flow regimes and artificial surfaces. Platelet consumption and activation is significantly increased, potentially resulting in hemorrhagic and thrombotic complications. Thrombotic deposition (blood clotting) within the device may lead to stroke or result in oxygenator failure as the pressure drop across the device increases. These complications may be minimized by eliminating areas of slow and recirculating flow, and by reduction of the overall oxygenator size. Oxygenator design improvement is currently an iterative process requiring fabrication of prototypes and extensive testing. Computational fluid dynamics (CFD) offers an investigator the potential to interpret the effect of oxygenator geometry upon blood flow characteristics, and indirectly, oxygenator hematologic biocompatibility.

Blood flow and oxygen transfer in a commercial hollow fiber membrane oxygenator has been modeled using a widely available CFD package. A 3D geometric model and surface mesh representing the Medtronic Maxima Plus was generated using I-DEAS (SDRC, Milford, OH), a commercial computer-aided design and engineering (CAD/CAE) package. The volumetric mesh was created in TGRID (Fluent, Inc., Lebanon, NH) and the governing equations were solved in FLUENT/UNS (Fluent, Inc., Lebanon, NH). Blood was modeled as a continuum with constant density and viscosity applicable at the flow rates under investigation. Since modeling blood flow around thousands of individual fibers in the oxygenator is currently an intractable problem, the fiber bundle was described as a porous region using a modified Ergun relation to approximate the flow resistance.

The oxygen flux to blood within the oxygenator was modeled using a mass transfer coefficient and the oxygen gradient between the fibers and the surrounding fluid. Oxygen transfer rates were computed on an element by element basis. The effect of local velocities and blood saturation was thus incorporated into the determination of local oxygen transfer. The variable mass transfer coefficient was determined by utilization of a correlation between the local Sherwood number and the local Reynolds and Schmidt numbers. Oxygen uptake by hemoglobin affected the oxygen diffusivity and flux and was accounted for by use of the Hill equation. The blood flow and oxygen transfer characteristics of the Medtronic Maxima Plus at clinical flowrates were computed using the analysis described above. Regions of low velocity flow were qualitatively correlated with clinically observed areas of thrombus formation within the oxygenator. In conclusion, CFD offers a methodology to investigate alternative oxygenator geometries which may lead to more efficient and biocompatible devices in the future.

The Element-Free Galerkin method has already been demonstrated to accurately model dynamic fracture problems. However, it did not have an obvious means of modeling contact across crack faces. The means of modeling contact will be briefly outlined, and results shown.

The Partition of Unity Finite Element Method (PU-FEM) has recently come to light, and been applied to static crack problems. It is a consistent means of extrinsically enriching a sub-set of finite elements in the model with linearly independent functions. When these functions are chosen to be a set which can construct the near-tip asymptotic crack fields, then an arbitrary crack can be accurately modeled, by even a relatively coarse mesh.

Since PU-FEM seems to be a robust, accurate method in statics, I have attempted to apply it to dynamic problems. The formulation of the problem has to be dealt with carefully, since a propagating crack causes the shape functions themselves to be time-dependent. The formulation will be explained, and results shown.

Systems that emulate human decision-making in an interactive environment include computer games and training programs. These systems typically employ such techniques as rule-based algorithms and utility function maximization. In certain applications, these traditional approaches have met with great success, with last year's chess victory by Deep Blue being perhaps the most famous example.

However, such systems have limited utility, because an experienced human operator will eventually discern the machine's rules and learn to out-maneuver them. Skilled human opponents typically adapt to one another's "rules," as well as to unexpected opportunities and set-backs, and are able continually to maneuver around one another in novel ways. In addition, efforts to create sub-expert computer opponents, such as for training purposes, usually consist of a crude stupefaction of the expert program. This often results in a player that makes mistakes no sane human novice ever would, and which is thus sub-optimal for training.

This project seeks to incorporate adaptive neural network models of human cognitive processing into a commercial game of strategy, in order to help build computer systems that more realistically mimic human players. This approach includes studies on human novices, and also draws from the literature on decision-making, reinforcement learning, and classical game theory. The ultimate goal is to develop a computer program that gradually adapts its playing style in a manner more like that of an improving human novice.

Researchers on this project, in collaboration with a commercial computer game manufacturer, have chosen a popular n-player board game as a testbed. Throughout the course of development, new computer game-players will be designed alongside and tested against the manufacturer's state-of-the-art rule-based system. The manufacturer has agreed to share proprietary software in exchange for early knowledge of novel cognitive information processing models that may be applied to the company's products.

For the last four years, I have been researching problems in computational ecology. Specifically, I have considered the application of parallel discrete event simulation to problems in demography. The consideration of discrete event simulation stems from two recent trends in ecological modeling, namely individual oriented modeling and spatially realistic modeling. In an individual oriented model, we model a population by modeling the physiological and social interactions of its members. Recent work has shown that spatial distribution of organisms tends to explain some demographic stability. This has led some ecological modelers including those interested in individual oriented models--to focus on the effect of spatial realism.

With few exceptions, modelers implement spatially explicit, individual oriented models as discrete time models. Discrete time models assume we can represent interactions by updating a model at regular intervals, but ecosystems are not generally synchronous when considering system behavior at the individual level. As a consequence, discrete event models may be more realistic and more efficient to implement.

In the initial phase of my research, I reviewed parallel simulation techniques and implemented a spatially explicit model using the Time Warp Operating System (Time Warp is a major parallel simulation algorithm). This experience led me to two conclusions: the performance of parallel simulation algorithms are not characterized well enough to comment on their appropriateness for a given problem and ecological models must be developed as discrete event simulations before we attempt to fit ecological models to a particular parallel algorithm.

The goal of my research is two fold. First I am working, in conjunction with ecologists, to develop discrete event ecological models. The primary product of this research is the Ecological Simulation System (ESS) a framework for population modeling based on Zeigler's Theory of Modeling and Simulation. The second goal of my research is a performance analysis of realistic ecological simulation, implemented with parallel simulation algorithms. The end of this research is to develop a general simulation tool for ecological modeling.

Nowadays, computers are used in many areas traditionally dependent on human expertise, such as medicine. Computer-assisted orthopeadic surgery (CAOS) is the new area of computer application in medicine. It consists of three components: modeling, planning and execution. One of the most common sources of preoperative information in CAOS is medical imaging. Advances in medical imaging technology (computed tomography (CT), magnetic resonance imaging (MRI), etc.) and in computer-based image processing and modeling have made it possible to geometrically model and

visualize anatomical structures in patients. These developments have given rise to the computer-based preoperative planners that are now appearing. The surgeon can interactively determine the operative procedure by manipulating computer-generated models of anatomical structures, surgical tools and/or implants. The computer can perform the simulations of the operative plan, giving the surgeon feedback on the proposed surgical procedure.

Existing preoperative planners do not relate the operation plan with the biomechanics of anatomical structures. Understanding the mechanics of orthopaedic surgery is particularly important, especially for cementless press-fitted implants. This work is the first attempt to incorporate patient-specific biomechanical simulation into a geometrical preoperative surgical planner. This can provide the surgeon with qualitatively new tools for computer-assisted surgery.

The complexity of anatomical structures requires the use of numerical methods, such as the FE method, to perform biomechanical simulations. Existing mesh generation techniques cannot automatically produce FE meshes of three-dimensional models of arbitrary geometry. Moreover, geometric modeling of patient-specific anatomical structure, which underlies existing preoperative planners, is still under development. A new method to automatically obtain three-dimensional FE models from CT scan data will be developed. This method utilizes mapping techniques, commonly applied for registration of MRI of the human brain. The method provides segmentation as a byproduct, which can be compared to the geometric planner's segmentation as a validation step. In the future, surface models obtained from FE models can also be used within geometrical preoperative planners. FE simulation of the acetabular component of Total Hip Replacement (THR) will be performed using patient-specific three-dimensional pelvic model of anatomical geometry.

Although this research is motivated by a specific problem of biomechanical analysis of the acetabular component of THR, the methods to be developed within this framework are general and can be applied to other biomechanical and medical imaging problems.

Recently Wavelet Transform based compression techniques have demonstrated, superior performance over Fourier (and other traditional) Transform based compression techniques in terms of high quality, low bit-rate compression, for many types of signals. Even more recently it has be shown that Adaptive Wavelet Packets (variant of the Wavelet Transform) based techniques have the promise of greatly outperforming the standard Wavelet Transform techniques. A potential hindrance to the use of the Adaptive Wavelet Packet techniques, in very low bit-rate applications, is the amount of overhead information needed to encode what is the "best" wavelet packet tree structure for processing each signal segment (block) in a sequence of many signal segments (blocks) to be processed. For example, in the compression of 2-D image blocks, if 4-level Wavelet Packet Tree Structures are utilized, then, 66 bits of overhead information are needed. For 5-level tree structures, 261 overhead bits are needed. The number of overhead bits (along with the versatility of the compression algorithm) increases greatly with depth of the tree structures considered. We have developed an Adaptive Wavelet Packet method that avoids sending this overhead information all together. We call this the Predictive Adaptive Wavelet Packet method. Initial results from this research indicate that the new Predictive Adaptive Wavelet Packet (PAWP) method can achieve a rate/distortion performance that is better than that of regular Adaptive Wavelet Packet (AWP) methods, in which overhead information is needed. The trade-off is that about twice as much computation is required for PAWP as for AWP. Results have been acquired for speech signals and image sequences.

A novel technique is proposed for analyzing transient electromagnetic scattering from three-dimensional inhomogeneous dielectric targets. The electromagnetic volume equivalence principle is invoked to construct an integral equation in terms of the electric flux density throughout the scatterer. The integral equation expresses a consistency relation for the fields in the scatterer's interior. To solve this integral equation, the electric flux is expanded in volumetric Rao-Wilton-Glisson basis functions on a tetrahedral mesh describing the scatterer. Using a Galerkin procedure, the integral equation is converted into a linear system of equations in terms of the basis function expansion coefficients. This linear system is solved using a marching-on-in-time scheme, which recursively evaluates the volume electric flux distribution at a given time step based on the external field that impinges on the scatterer and the flux distributions at previous time steps. For validation purposes, sphere scattering data obtained using this algorithm are compared to analytical Mie series solutions. In addition, results for scattering from more complicated objects are compared to those obtained with a frequency domain method of moments solution via Fourier transformation. The acceleration of this algorithm using the plane wave time domain algorithm is currently being studied.

Gamma-ray bursts (GRBs) are the most energetic cosmological events in the universe since the Big Bang. During their few seconds of existence, GRBs may produce as much as energy as millions of galaxies.

In the late 1960s, the U.S. Department of Defense and the U.S. Atomic Energy Commission launched a network of Vela satellites to monitor the nuclear test ban treaty by looking for gamma radiation from atomic explosions. Brief, randomly distributed explosions of gamma radiation from space were detected. These "gamma ray bursts" were announced in 1973, and since then, scientists have been trying to explain the physical mechanisms that could create such cataclysmic events.

The Antarctic Muon and Neutrino Detector Array (AMANDA) is a neutrino detector located at the South Pole. Neutrinos, which are chargeless and nearly massless, are notoriously difficult to detect. On occasion, one will interact with rock or ice in the earth and produce a muon. The relativistic motion of the muon creates an electromagnetic shock wave called Cherenkov radiation.

AMANDA consists of hundreds, and eventually thousands, of photomultiplier tubes (PMTs) deployed between one and two kilometers deep in the ice. The AMANDA PMTs detect the Cherenkov radiation from the relativistic muons. By recording the timing and topology of PMTs triggered in an event, we can reconstruct the paths and directions of the muons and neutrinos.

The detection of neutrinos associated with GRBs would provide important information on the physics of GRBs and whether they may be responsible for another astrophysical enigma, ultra high energy cosmic rays.

AMANDA uses data from the Compton Gamma Ray Observatory (CGRO), an orbiting gamma radiation satellite detector, to assist our search for neutrinos associated with GRBs.

The AMANDA detector recorded several hundred gigabytes of raw data in 1997. This poster will describe how that data was filtered for various research goals at the National Energy Research Scientific Computing Center (NERSC) at Lawrence Berkeley National Laboratory. I will also describe how the data filtered for our GRB studies was reconstructed and analyzed at the University of Wisconsin-Madison, and how data from the CGRO satellite was incorporated in our work. The current status of GRB neutrino research will be summarized.

A fundamental problem for bootstrap current-driven magnetic islands is to understand the dynamics and profile properties of density and electron temperature. Although electron temperature flattening inside the island separatrix has been observed, it is uncertain the degree to which density profiles show a similar flattening. To study the evolution of density and electron temperature in a helical island geometry, we solve the coupled parallel momentum and continuity equations in the presence of a slower-growing magnetic perturbation that simulates the evolution of a neoclassical tearing mode. It is assumed that sound wave propagation along field lines is primarily responsible for equilibration, of density over perturbed flux surfaces. Electron temperature equilibration, on the other hand, occurs on perturbed flux surfaces as a result of a rapid parallel heat flux. Properly accounting for this heat flux requires a solution of an appropriate kinetic equation that takes into account free-streaming and collisional effects in a helical magnetic geometry. An analytic closure for the heat flux is constructed based on a multiple time and spatial scale scale, Chapman-Enskog-like analysis. This heat flux is then inserted into a temperature evolution equation which is solved in the presence of an evolving helical magnetic island.

Active flutter suppression has been a topic of interest for some time and computational methods have been used to simulate the time response of Aeroservoelastic systems. Though these simulations are valuable for validations and other non-linear problems, time-integration methods do not easily provide all of the information required for controller design for an aeroservoelastic system. For such purposes it is necessary to produce a reduced order model of the system to be controlled. In the current context, one of the requirements for a reduction technique is that it be derived from a computational model analogous to the ones used for time integration. Several approaches have been recently developed to tackle this problem, notably methods designed to either reduce the fluid part of the coupled system or the complete aeroelastic system.

The purpose of an active control system for flutter suppression is to stabilize the aeroelastic coupled modes that are unstable. Therefore it seems natural to use the second model reduction technique as the basis for the representation of the complete system to be controlled. However, the aerodynamic of the actuation by way of the motion of a control surface has a characteristic time that is rather small when compared to the period of one oscillation of the coupled system. Therefore we suggest that the two aforementioned model reduction techniques can be successfully combined for the analysis and design of an aeroelastic system. In this poster, we present an approach that allows the designer of an active control system, or active flutter suppression system, to obtain a reduced order system that combines the representation of the aeroelastic modes to be controlled and of the aerodynamic of an actuation by control surface motion. We also present a simple example of the technique, as well as our current progress regarding theoretic stability analysis and a more complex aeroservoelastic system example.

Many nonequilibrium problems in physics and engineering have natural mathematical formulations as hyperbolic systems with relaxation source terms. Such systems pose interesting challenges for numerical approximation. The source terms cause the system to be dispersive; generally, both the eigenvalues and eigenvectors of the system are modified by the relaxation processes which drive the system towards equilibrium. Another difficult numerical issue arises when the relaxation source terms operate on much smaller time scales than the advection terms; the problem is then said to be stiff. Unfortunately, many problems simultaneously exhibit behavior in both the stiff and non-stiff limits, as well as in between. It is therefore desirable to develop high-resolution numerical algorithms which are uniformly accurate irrespective of the scales of the relaxation source terms. For example, if the data are such that the advection terms permit a time step large enough to bring the flow into local equilibrium, one would wish to see an accurate solution of the equilibrium problem emerge without resolving the details of the transition.

In recent years, there has been much development on such uniformly accurate methods. A variety of splitting schemes, characteristic schemes, and staggered, centered schemes have been proposed. The most promising of these rely on staggered grids to eliminate the difficulty of resolving the Riemann problem, which no longer has a self-similar solution. However, one drawback of centered schemes is that one may lose definition of contact discontinuities, which are likely to be very important in reactive flows. The question still remains as to whether there is a simple, natural extension of the upwind approach which includes the effects of the relaxation source terms and resolves all discontinuities equally well. The close coupling between the advection and relaxation has lead us to explore an upwind scheme based upon an approximate solution to the Riemann initial value problem. Currently, experiments are being performed which apply this method to a model linear system as well as to equations of extended hydrodynamics.

A few problems related to groundwater flow are considered. One particular problem concerns groundwater mound formation and extension following a flood. A flood, which can follow a break-through of a dam, can have a very significant environmental impact if contaminated floodwater enters and then extends far into the soil.

The mathematical model for groundwater flow, which is described by the porous medium equation, has many interesting properties. For the case when the pores do not retain water, it can be shown analytically that there exist special self-similar solutions. The graphs of such solutions at different times are related by a similarity transformation. Such solutions are not simply special analytic solution, but they also provide asymptotics for a class of initial and boundary value problems. However, for a porous medium with capillary retention, the existence of self-similar solutions cannot be established analytically. My co-author and I performed numerical computations that showed that special self-similar solution exists for this problem too.

We also considered the problem of controlling the groundwater mound extension after the flood by using forced drainage. It has been shown analytically that the porous medium equation with forced drainage also has self-similar solutions for some special forced drainage rates. However, there are no self-similar solutions that model the constant forced drainage rate. Numerically, we computed non self-similar solutions with constant drainage and showed that, in this case, it is possible to extinguish completely the water mound that appeared after the flood.

We also analyzed the problem of flooding and forced drainage for the case of porous medium with cracks. In many practical cases, this is a more realistic model than a homogeneous porous medium. The groundwater moves much faster in cracks than in the pores. However, also the volume of cracks is very small compared to the volume of pores. Thus, the inclusion of the cracks in the model significantly changes the results obtained above.

Neutrophils extend pseudopodia in the presence of N-formyl-methyl-leucyl-phenylalanine (fMLP), which is a by-product of bacterial metabolism. One theory that attempts to explain the mechanisms by which the outward pushing of the pseudopodia evolves, is the swelling model. In this model new actin filament networks, which form underneath the site of fMLP binding, swell due to interfilament repulsion and filament-solvent attraction, thus forming a pseudopodia. Using finite element analysis, we have implemented the swelling model to simulate pseudopodia extension and retraction and neutrophil locomotion in a micropipet. Our results indicate that the swelling model can explain a number of observations seen in the experimental setting. These include data on the kinematics of pseudopodia extension and retraction, the distribution of the cytoskeleton, the geometry of extending and retracting pseudopodia, and the physical force generated by a neutrophil during active movement.

The steady advance of both algorithmic and computational power in the last decade is leading us to an age when Molecular Simulation will be able to routinely, and quantitatively predict the fluid phase diagrams of industrially important molecules. This talk will briefly discuss the current state of the art of Configurational-bias Gibbs Ensemble Monte Carlo simulation, and present recent applications of this technique to the single and multi-component vapor-liquid phase coexistence of linear and branched alkanes at sub and super-critical temperatures and pressures. The partition constants of solutes between a squalene alkane liquid phase, and a Helium vapor phase are presented as an example of how we can compute the relative retention times of solutes in a chromatographic column.

The adsorption of phosphine on the Si(100)-2x1 surface is a difficult system to model due to the non-local character of the electron transfer arising from the formation of the Si-P dative bond. We first model this non-local effect with a couple of clusters, Si₉H₁₂ and Si₂₁H₂₀, which model the surface along the dimer row; and with a third cluster, Si₂₃H₂₄, which models dimers separated by the surface trench. We characterize some of errors involved in the cluster calculations, associated with cluster constraints, the basis set size, and the level of theory. We compare the fully relaxed cluster with an over-constrained cluster to estimate upper and lower bounds of the potential energy surface. We compare the effects of a very large basis set with the split basis set we used for the majority of our computations. Finally, we calculate the errors in the level of theory by comparing the Becke3LYP DFT method we used with the QCISD(T) method. We find that with this non-local chemistry, cluster size is important while constraint effects are not usually important. The errors due to the basis set and the level of theory are small for large clusters.

In recent years, complex, spatial robotic manipulators and mechanisms have attracted the interest of machine designers because they can combine high dexterity with complex three dimensional trajectories, relatively large workspaces and high accuracy. Examples of such devices include prosthetic devices, high speed assembly machines and walking robots. The kinematics analysis and synthesis are two important steps in the design of these mechanical devices by studying the geometrical properties of them. The analysis problem studies the mechanism capabilities and performance when its type and geometric parameters are given. The synthesis problem is the process of designing a mechanism to accomplish a specific task. The mathematical formulation of both problems lead to complex systems of multivariable polynomial equations. Due to their nonlinearity, multiple solutions exist. It is not uncommon that a kinematics problem possesses hundred or even thousand of solutions. These multivariate polynomial systems are, therefore, not easy to solve and many important design problems remain unsolved to this date.

The methods used in solving polynomial systems can be classified as either analytical or purely numerical methods. Analytical methods are those that, by carefully creating and subtracting new equations from old ones, eliminate all but one variable in the multivariate polynomial system and reduce it to a single univariate polynomial, which is then solved by a numerical method. These methods involve the manipulation of large amount of intermediate symbolic information and large integer coefficients. Examples of these methods are the classical elimination method and the Grobner Basis Method which are based on modern theory of commutative algebra. Purely numerical methods rely on iterative algorithms which compute the numerical solutions directly from the polynomial systems. One of the most popular numerical methods is the polynomial continuation method which computes all the possible solutions by tracking all the possible solution paths. Thus, both classes of method involve large amount of computation and storage of information. Solving such problems is not possible without the careful use of modern computing technology. The main objective of this research is to apply existing computational methods as well as develop generic methodologies for the kinematic analysis and design of robotics systems and mechanisms. Using these computational methods, analysis and synthesis algorithms for different types of mechanisms will be developed and used in design automation and control of them.

I have extended an existing simple Lagrangian model developed by Yamazaki and Haury (1993) which uses a nonparametric density estimation technique to study animal aggregation. The essence of the original model is its two parameters: a, which describes the perception distance of the organisms, and b, which relates to their motivation to aggregate. In the original model an attractive component caused individuals to move toward higher densities of organisms with strength b. This revised model also includes a tendency to be repulsed, instead of attracted, by densities of organisms above a certain threshold level, represented by a third parameter, g. Numerical simulations were performed in two-dimensional space using 100 organisms and a range of values for the three parameters. While including repulsion leads, in general, to more dispersal of organisms with time, the interplay between perception distance and motivation to aggregate can be quite complex with some surprising results, including initial super-diffusion.

Molecular dynamics programs simulate the behavior of biomolecular systems, leading to insights and understanding of their functions. However, the computational complexity of such simulations is enormous. Parallel machines provide the potential to meet this computational challenge. To harness this potential, it is necessary to develop a scalable program. It is also necessary that the program be easily modified by application-domain programmers. The program NAMD2 seeks to provide these desirable features. It uses spatial decomposition combined with force decomposition to enhance scalability. It uses intelligent periodic load balancing, so as to maximally utilize the available compute power. It is modularly organized, and implemented using a parallel C++ dialect, so as to enhance its modifiability. It uses a combination of numerical techniques and algorithms to ensure that energy drifts are minimized, ensuring accuracy in long running calculations. NAMD2 uses a portable run-time framework that also supports interoperability among multiple parallel paradigms. As a result, different components of applications can be written in the most appropriate parallel paradigms. NAMD2 runs on most parallel machines including workstation clusters.

Monitoring atmospheric air pollutants is an issue of primary importance in our society. Toxic air pollutants are poisonous substances in the air that come from natural or manmade sources, and can harm the environment or our health. Two key issues need to be addresses when mapping air pollutants: a) the distribution of air pollutants exhibits a high variability in both space and time, and b) the measurements of air pollutants at monitoring stations may be uncertain. The Bayesian Maximum Entropy (BME) method of modern Geostatistics provides an adequate and rigorous framework for this type of analysis. BME uses the theory of Space/Time Random Field S/TRF that accounts for the space/time variability of air pollutants. Additionally BME incorporates uncertain information (also called soft data) in a rigorous manner, providing more accurate mapping estimates than the classical kriging methods when a combination of hard (exact) and soft information is used. This approach is demonstrated on a case study concerned with mapping ambient PM10 in the state of North Carolina.

This poster will discuss the use of Krylov methods in magnetohydrodynamic simulations. These methods use nonlinear exponential propagation to perform time evolution and approximate the action of the right-hand-side differential operator and its Jacobian by projections to the Krylov subspace via the Arnoldi process. The primary advantage of Krylov methods is that they allow one to compute the evolution of equations using a time step which exceeds the CFL limit imposed by the fast modes in the system.

Another advantage of the method is the use of Krylov subspace projections to propagate the right-hand-side differential operator, which avoids the inversion of large complex matrices and makes the method highly parallelizable. We will describe this numerical technique and discuss its performance for the applications in plasma physics.

Droplet splashing occurs in a wide variety of applications, including spray painting processes, injection of liquid fuel in internal combustion engines, fire suppression systems, and many others. The subject is definitely not new and has received attention by numerous researchers for over a century. Due to its complexity, however, characterized by large deformations over very small time scales, propagations of instabilities in the emerging liquid crown and sheet, formation of jets along the rim, and subsequent jet break, the study of the splashing phenomena is far from complete.

Our interests in this area are the determination of the velocity and size characteristics of secondary droplets created during splashing, as well as the calculation of quantities that described the overall splashing event (e.g. the amount of energy lost due to viscous dissipation). With an emphasis towards automotive applications, the single droplet splashing study presented here provides a foundation for the investigation of spray impingement on containment walls, where thousands of droplets collide against a solid surface.

Our approach is based on describing the products of splashing in terms of a joint probability density function (jpdf) of an assumed form, and solving this jpdf by applying the conservation laws to the splashing process. Drop size calculations from this model are compared to experimental drop size data yielding acceptable agreement. The distribution of incident droplet energy among the kinetic and surface energy of the secondary droplets, the residual film surface energy, and the time integrated viscous dissipation are calculated as a function of the incident Weber number and the impingement angle. Information is also obtained on the incident angle dependence of the impact forces.

Lithium transition metal oxides (Li_x MO₂ with M = transition metal) are ideally suited to serve as cathode materials in the next generation of rechargeable lithium batteries. A lithium battery consists of an anode (which is metallic lithium) and a cathode (the lithium transition metal oxide) separated by an electrolyte. In use, lithium ions are released at the anode, travel through the electrolyte and are intercalated into the cathode material. Lithium transition metal oxides exhibit a high voltage, can undergo large variations in lithium concentration and are stable under repeated removal and insertion of lithium ions. Currently the most important cathode material is Li_x CoO₂ . In this work we have investigated the thermodynamic phase stability of Li_x CoO₂ as a function of the lithium concentration. We have also calculated the mobility of the lithium ions within the CoO₂ crystal structure. Both quantities are essential to enable an intelligent optimization and design of existing and new cathode materials.

It has been more than two decades since the introduction of Smoothed Particle Hydrodynamics (SPH) by Lucy [1] and Gingold & Monaghan [2]. Initially used to simulate astrophysical problems including galaxial collisions and the dynamics of star formation, SPH has recently been adapted to many relevant chemical engineering problems, including heat and mass transfer [3], molecular dynamics [4], and fluid & solid mechanics [5].

Because SPH is a Lagrangian method, novel SPH algorithms will have a definite computational advantage over conventional methodologies: the elimination of a domain grid structure, which for complex flows, must be continuously rebuilt and realigned. In conventional methods, a significant amount of computation time is spent in creating and redistributing mesh points, and typical computational fluid dynamics calculations are inefficient and require mass data storage. The absence of a grid also means that 3-D calculations are as easy to implement as 1-D. Additionally, the extension of the SPH methodology to multi-phase systems could be inherently important to technologies such as particulate control and collection, combustion/chemical reaction processes, and the design of micro/nano-electromechanical systems.

To simulate any physical system, the best methodology could be an N-particle dynamics simulation. Even though an N-particle simulation is logically straightforward, it is extremely difficult to implement in engineering systems because of the enormous number of degrees of freedom, and strong, long-range interactions between particles make N-particle simulations virtually impossible. By reducing the degrees of freedom and averaging the N-particle dynamics, one often adopts continuum theory to tackle complex engineering problems. The continuum-mechanical approach is best-suited for the numerical simulation of single-phase systems. However, in handling hybrid systems composed of continua/dispersed particles or chemically-reacting flows, continuum theory may not be the optimal scheme for numerical studies. An alternative scenario for the description of physical systems is the SPH technique, which casts either the N-particle system or the continuum into a finite number "n" (n << N) of "pseudo-particle" aggregates. SPH not only reduces the degrees of freedom to a numerically-manageable number (from the N-particle picture), but it also eases the implementation of parallel algorithms (from the continuum picture) because the partial differential equation representation (field view) is converted into ordinary differential equations (particle view). Recently, there have been two frontier achievements in the application of SPH to fluid mechanics. The original work was initiated by Takeda, et. al. [6], who converted the Navier-Stokes equations to the SPH framework and focused on relatively high-Reynolds number flows. Morris, et. al. [7] extended Takeda's analysis to low-Reynolds number situations and was the first to use an (anti)-symmetrized form for the viscous flow term, including artificial viscosity. In this paper, we will generalize and unify these two methodologies by focusing on and clarifying the following three unresolved issues found in the previous work. The first issue addresses boundary conditions. Both Takeda, et. al. and Morris, et. al. introduced fluid "ghost" particles outside of the actual fluid domain to simulate the fluid/solid boundary. The use of ghost particles may be reasonable and efficient in simple cases (e.g., rectangular coordinates) but is extremely difficult for more complex problems (curved coordinates, etc.). Therefore, there is a dire need to invent a new, physical, and numerically-elegant boundary condition methodology. This will be especially helpful when one expands the SPH fluid mechanics technique to include multi-phase flows.

1. L.B. Lucy, Astron. J. 82, 1013 (1997).
2. R. A. Gingold and J. J. Monaghan, Monthly Notices Royal Astron. Soc. 181, 375 (1977).
3. W. Moran, W. G. Hoover, and S. Bestiale, J. Stat. Phys. 109, 709 (1993).
4. W. G. Hoover and C. G. Hoover, Adv. Comp. Meth. Materials Modeling 180, 231 (1993).
5. L. Libersky, A. G. Petschek, T. C. Carney, J. R. Hipp, and F. A. Allahidadi, J. Comp. Phys. 109, 67 (1993).
6. H. Takeda, S. M. Miyama, and M. Sekiya, Progress Theoretical Phys. 92, 939 (1994).
7. J. P. Morris, P. J. Fox, and Y. Zhu, J. Comp. Phys. 136, 214 (1997).

The organometallic vapor phase epitaxy (OMVPE) process is widely used to produce compound semiconductors for applications ranging from solar cells to cellular phones. The OMVPE process is a chemical vapor deposition (CVD) process where thin films are synthesized from organometallic precursors. As device dimensions decrease, processing issues that affect the material properties (i.e. stringent demands on film thickness, morphology and impurity incorporation) become increasingly important. Models of OMVPE can be used in fundamental studies of the processing of thin films, as well as in predictive simulations to aid in process design and control. The difficulties in modeling the process arise from the complex nature of OMVPE; processes occur on widely different length and time scales. In this work, models on various length and time scales for the OMVPE process will be developed and methodologies for linking these models will be formulated. The "linked" models will be shown to make accurate predictions on various length scales. The complex coupling of gas-phase transport phenomena with gas-phase and surface chemical kinetics means that gaining an understanding of the entire OMVPE process requires more than one type of modeling approach. Predictions on the macroscopic length scale (growth rate, film uniformity, film composition) and on the "mesoscopic" length scale (surface morphology) are both needed to determine optimal operating conditions for OMVPE reactors. Models have been developed to make predictions on each of the relevant length scales:

• Macroscopic reactor scale models have been developed to understand the fluid flow, heat and mass transfer combined with detailed gas-phase chemical kinetic mechanisms to predict the type and concentration of precursors arriving at the growth front. Finite Element simulations have been developed in our group to solve these problems.
• "Mesoscopic" surface models have been developed to study the evolution of surface morphology during growth. Kinetic Monte Carlo methods have been used to investigate the importance of individual steps in surface mechanisms, as well as determine the growth mode on the surface. Parallel algorithms are also in the process of being developed in order to speed the computations of the surface problem.

A methodology is developed to "link" reactor scale and surface models in order to gain an overall understanding of the OMVPE process. Case studies of GaAs growth using a variety of precursors is used to validate the models separately, as well as to validate the "linked" model. Using the reactor scale model alone cannot give predictions of surface morphology. The surface model used alone is missing a needed input of gas-phase flux to the surface. The "linked" model is able to make extremely accurate predictions on both the macroscopic length scale (growth rate) and mesoscopic length scale (surface morphology).

Void growth is one of the primary causes of failure in small circuit elements sustaining high current densities, and is driven by vacancy diffusion, electromigration, surface tension, stress, and grain geometry. I will present a 2-D model (appropriate for thin metallic conducting lines) which neglects stress and grain geometry; in particular, I will demonstrate the use of level set methods to track interfaces, and will describe our method of solving elliptic PDE on an irregular domain using a uniform grid which does not line up with the domain boundary.

Models of evolutionary ecology often assume fixed, externally imposed selection pressures, while global climate models often ignore the effects of evolution. In reality, plant and animal species work together to modify the planet's atmosphere and terrain, while these modifications affect the selective pressures that determine what future species will do. The feedback between biotic modification of the environment and environmental selective presures on the biosphere is an underexplored aspect of climate and ecological modeling. We describe a model which allows us to explore these feedbacks, and present some model results.