Multiphysics simulations

2.1.2 Fission reactor fuel performance

Simulating the performance of light water nuclear reactor fuel involves complex thermomechanical processes between fuel pellets, which are made of fissile material, and the protective cladding that surrounds the pellets (Hansen et al., 2009b). Figure 3 shows the relationship between the fuel pellet and cladding and the helium-filled gap separating them. Advancing fuel design and analysis capability requires progress both in the modeling of the fuel material behavior and in the development of a multidimensional computational framework that couples these models in a consistent, accurate, and efficient manner to describe how the fuel geometry and behavior change over the lifetime of the fuel.

OPEN IN VIEWER

Figure 3.

Illustration of a multiscale multiphysics calculation to predict the behavior of reactor fuel cladding containing hydrogen inclusions (hydrides). In this process, an initial guess of the stress state in the fuel cladding and the temperature history are calculated by using a macroscale calculation of the fuel rod. These states, defined at the integration points of the multiscale calculation, are passed to concurrent phase field mesoscale calculations in order to determine the morphology of the hydride inclusions. This morphology is then passed back to the mesoscale calculation and used to calculate the stress profile within the cladding.

The base physics begins with heat diffusion that partially drives the mechanical response of the fuel pellet and cladding system (Newman et al., 2009a). Thermomechanical response problems are often solved in an operator-splitting fashion, but this approach can result in convergence and stability issues as the pellet and cladding come into contact as the fuel ages. Newman et al. (2009b) solve the fuel performance problem in a fully coupled manner on a single mesh. Further, an oxygen diffusion equation is added since, under hyperstoichiometric conditions, both the thermal conductivity and mechanical properties of the fuel pellet are impacted by the diffusion of oxygen within the matrix. Additional material models are used to describe other details of the thermal and mechanical behavior of the fuel and cladding, such as swelling of the fuel due to the accumulation of fission products within it (Williamson and Knoll, 2009). The expansion of the pellet and, indirectly, the width of the gap separating the pellet and the cladding are functions of the thermal profile within the pellet, the age of the fuel (the degree of fission-product-driven swelling), and the creepdown of the cladding onto the outside surface of the pellet.

Given the significant radial and axial displacement of the pellet relative to the cladding, it is infeasible to mesh the helium-filled gap between the two, especially in a multiple pellet simulation. Since the gap is long and narrow, Rayleigh numbers are unlikely to be high enough within the gap to support significant convection. However, it is important to accurately model the conductive heat transfer across the gap, since the gap is one of the principal factors that influence performance of the fuel. Thus, a fully coupled thermomechanical contact model was developed to model the evolution of the gap within the fully coupled multiphysics code (Hansen, 2011).

Irradiation of the fuel and the cladding affects the thermal and mechanical properties of both. It is difficult to describe these effects on the coarse scale of the finite-element code described above. Modeling the irradiation of the material is more easily done at lower-length scales, such as the mesoscale. One might employ a phase-field simulation approach to calculate the evolution of the thermal conductivity of the fuel under irradiation. Use of this model can be challenging, however, since the characteristics of the irradiation damage change with fuel temperature. Tonks et al. (2009, 2010) couple a mesoscale solution of the Cahn–Hilliard equation to the macroscale approach to form a multiscale simulation. The mesoscale simulation requires the coupling of a temperature field from the macroscale as input to the simulation. Further, the thermal conductivity calculated at the mesoscale impacts the temperature field calculated on the macroscale. A similar approach, shown in Figure 3, may be used to assess the safety of handling used fuel after drying and storage. Here, hydrogen inclusions (hydrides), which may weaken the cladding surrounding the fuel pellets, can form during the drying of fuel. To predict cladding performance, a macroscale structural and thermal response analysis is linked to a mesoscale calculation that computes the structure of the hydride inclusions (the hydride morphology). The stress and temperature at the integration points serve as boundary conditions for the mesoscale phase field calculation hosted there; the temperature and temperature gradient in time drive the precipitation of hydrogen from solid solution; and the stress state during precipitation affects the orientation of the hydride inclusions. The resultant hydride morphology calculated at the mesoscale is then upscaled for input to the mechanical response macroscale calculation that predicts the cladding behavior under loading scenarios of potential concern.

The equations that describe reactor fuel performance are solved by using a fully coupled, JFNK approach (Gaston et al., 2009a,c; Newman et al., 2009b) (see Algorithm 3 and Section 3.1.2). A fully implicit time-integration strategy is used to support accurate time integration needed over short transient events, such as the initial warmup of the reactor and power transients during operation, and to allow long timesteps to be taken under the long quasi-static periods of normal operation. All the phenomena within the reactor fuel are tightly interrelated, particularly the mechanical aspects of the fuel. As the fuel warms up, the pellet expands, thereby narrowing the gap and increasing the pressure of the gases between the pellet and cladding. The narrowing gap increases the heat transfer efficiency across it, which cools the pellet, checking its expansion. Attempting to decouple these processes in order to simplify the solution process can result in convergence issues and even “numerical chatter”, where the pellet and cladding chatter on incipient contact instead of arriving at an equilibrium solution (Hansen, 2011). In the gap-integration technique presented here, fluxes and tractions are exchanged between the pellet and cladding using the mortar finite-element method. The integration space and Lagrange multiplier implementation employ the Moertel package (Gee and Hansen, 2012) in Trilinos.

2.1.3 Conjugate heat transfer and neutron transport coupling in reactor cores

It is well understood that nuclear fission energy is a key component of our current and future energy needs. It is therefore urgent to develop nuclear reactors that are safe and environment-friendly. Accurately predicting neutron densities in a reactor core is a critical aspect of design and safe operation. Reactor cores are designed to balance power generation with efficient heat transfer and fuel burnup objectives, optimizing energy density relative to design limits for key global properties such as peak temperature and temperature gradient, neutron fluence, and pumping power (pressure drop).

With the availability of leadership-class supercomputers, high-fidelity simulations of different physical phenomena (such as neutron transport, coolant flow, structural mechanics, and materials science) inside complex reactor geometries appear a feasible and attractive option. However, the modeling and simulation (M&S) tools need to accurately represent the coupling among these physical processes while preserving the accuracy in each physics. The existing M&S tools are based on various approximations (to reduce the computational burden) and rely on various correlations derived from experiments that can be applied for a narrow range of design space. Therefore, what is needed is to develop a modern, integrated design tool that allows automated data transfer among different design activities, and facilitates enhanced global design optimizations. Work in the Simulation-based High-efficiency Advanced Reactor Prototyping (SHARP) project (Siegel et al., 2007; Pointer et al., 2012a) and the Center for Exascale Simulation of Advanced Reactors (CESAR) (Rosner et al., 2012) focuses on developing such tools for reactor design.

On relatively fast timescales (compared with the physics of neutronics/fuel coupling) and within certain parameter regimes, methods for computing the accurate coupling between the neutron transport and conjugate heat transfer, for fixed fuel and structural configuration, become the key problem of interest for reactor analysts. Historically, reactor design and analysis tools have been based on weak coupling between bulk fluid characterizations and few-group diffusion or nodal neutronics methods.

The SHARP and CESAR projects are investigating high-fidelity fluid/neutronic coupling methods for a range of reactor phenomena. Many of the calculations require large numbers of energy groups and are crucial for simulating various accident scenarios, such as ejected control rod and loss of coolant flow without shutdown. These accidents, which involve rapid thermal–mechanical changes in the reactor core geometry, require several minutes to hours of time data. The existing tools do model these types of problems, but confidence in the solutions is difficult to assess. Higher-fidelity solution capabilities will not only allow better understanding of the behavior of existing reactor systems but will also predict the behavior of many newly proposed systems with as yet untested design characteristics (i.e. those without direct experimental data to back up the approximations made in the simplified modeling). High-fidelity solutions can also be used to enhance the existing simplified models by providing improved methodologies for obtaining reactivity coefficients and heat transfer coefficients.

The neutron transport equation is volumetrically coupled nonlinearly to thermal hydraulics via fuel and coolant temperature and density. Both neutronics and thermal hydraulics are coupled to the structural behaviors of fuel rods (and other components as shown in Figure 4). At present, the SHARP and CESAR projects are studying the coupling between neutronics and thermal hydraulics. The details of the discrete formulations of the governing equations are reported in companion papers (Fischer et al., 2007, 2008; Kaushik et al., 2009) and are not essential to understanding the intermodule coupling problem.

OPEN IN VIEWER

Figure 4.

Illustration of multiphysics coupling in SHARP. Left: Schematic representation of interrelationships among the different modules of SHARP. Here Ω n denotes the neutronics (UNIC) mesh, and Ω t h denotes the thermal-hydraulics (Nek) mesh (Siegel et al., 2007). Right: Predicted coolant velocity distribution in a light-water reactor (Pointer et al., 2012b). Evaluating reactor design performance accurately depends on resolving multiphysics coupling between neutron transport and coolant flow through the complex reactor geometries.

When the physics is split into separate components, these couplings take the form of specific data values that are interpolated from source to target meshes and sent between the modules through well-defined interfaces. To solve the governing equations, SHARP includes two physics modules: (1) Nek (Fischer et al., 2002, 2012), a spectral-element code that solves the three-dimensional, incompressible (Boussinesq), time-dependent Navier–Stokes equation with conjugate heat transfer on unstructured, higher-order quadrilateral meshes; and (2) UNIC (Kaushik et al., 2009; Smith et al., 2011), an unstructured, deterministic neutron transport code that incorporates parallel even parity, sweeping, and ray-tracing algorithms.

Separating physics modules into distinct components and implementing coupling as interactions between these components impose a number of requirements on the overall design of the SHARP framework, discussed in Siegel et al. (2007). SHARP couples physics components and other services through the spatial discretization or mesh abstraction provided by MOAB (Tautges et al., 2004, 2011b,c), which is introduced in Section 4.1.3.

New research in the CESAR project is building on this foundation to develop a coupled, next-generation nuclear-reactor-core simulation tool capable of efficient execution on exascale computing platforms. The approach incorporates tight coupling of extreme-fidelity models of neutron transport and conjugate heat transfer.

2.1.4 Multiscale methods in crack propagation

Most physical phenomena of interest to humankind involve a range of temporal and spatial scales. For example, the brittle fracture of solids under external stresses can lead to the failure of large structures such as bridges or ships, but originates at the atomic scale with the propagation of microcracks. Likewise, tidal currents in bays extend over many miles, but their behavior is dictated by the water viscosity, determined from the molecular-scale interactions of water molecules. A computational approach that captures the complexity of such physical systems must be able to treat all the relevant scales at the level of accuracy required for a realistic description. In this section and the next, we illustrate the issues involved in such multiscale simulations through two examples that require a broad range of scales, from the continuum to the atomistic and even to the electronic level. (Lu and Kaxiras, 2005, review similar problems and methods for their simulation.) The two examples are chosen to demonstrate both the progress that has been made and important challenges that remain in the areas of hard and soft materials.

In this section, our example involves hard materials, specifically the modeling of crack propagation and fracture of solids. Crack propagation is often the result of chemical impurities (corrosion), which can change the mechanical behavior of a solid from ductile (tough, difficult-to-break solids) to brittle (weak, easily breakable ones). The length scale at which this phenomenon is observed in everyday situations is roughly millimeters and above. At this scale cracks are detectable by the naked eye and, if not arrested, can lead to the failure (occasionally catastrophic) of large structures such as ships, airplanes, and bridges. The presence of chemical impurities, often in minute concentrations (a fraction of a percent or lower), changes the nature of the bonds between the atoms, which is crucial in the region near the tip of the crack where bonds are being stressed to the point of breaking. These atomic-scale changes are ultimately responsible for the macroscopic-scale change in the mechanical behavior of the solid.

A schematic breakdown of the different length and scale regimes relevant to this phenomenon is illustrated in Figure 5. The macroscopic scale (10 − 3 m or above) on which the crack propagation is of practical interest contains a huge number of atoms, of order 10 24 , and can be handled only by continuum approaches such as elasticity theory. The response of the solid to external stresses is determined by the motion of defects such as dislocations, which provide the mechanism for deformation. Their distribution and motion at the mesoscopic scale (10 − 5 m) can be described by theories that explicitly take into account these defects and the pattern of stress they create within the solid (Peierls, 1940; Nabarro, 1947; Hirth and Lothe, 1992; Devincre and Kubin, 1997). The origin of dislocations and the nature of the interactions between them can be described only by considering their structure at the microscopic scale (10 − 7 m) in terms of atomic degrees of freedom interacting by springlike forces (Bulatov et al., 1998). Furthermore, the presence of chemical impurities and their effect on atomic bonding, and by extension of the forces that hold atoms together near the tip of a crack, must include a quantum-mechanical description of the ion–electron interactions that can capture the chemical properties of different types of atoms. This scale can be treated in the context of density functional theory using, for example, methodologies such as the Car–Parrinello approach (Car and Parrinello, 1985). A first attempt to capture all these scales for the case of fracture in silicon, a prototypical brittle solid, was the multiscale atomistic ab initio dynamics approach, which was able within the same method seamlessly to couple the quantum-mechanical scale to the classical atomistic and to the continuum (Abraham et al., 1998). A direct link between the atomistic scale and the scale of defects (dislocations, domain boundaries, etc.) dynamics was accomplished by the quasi-continuum method (Tadmor et al., 1996), which enabled realistic simulations of the large deformation of materials (Smith et al., 2001; Tadmor et al., 2002). These examples addressed, in an effective manner, many of the computational issues in the coupling of length scales for multiscale simulations of materials.

OPEN IN VIEWER

Figure 5.

Illustration of the length and timescales in crack propagation. Left: Schematic representation of the length and timescales involved. Right: Multiscale simulation of crack propagation in silicon, a prototypical brittle solid, which includes all the relevant length scales; the color maps show the distribution of stress (highest at the tip of the crack), with stress waves passing seamlessly from one region to the next (Abraham et al., 1998).

In both this hard-material example and the soft-material example, which is described in the next section, there are a number of common issues that affect the accuracy and tractability of simulations. These include boundary conditions, model stability, error propagation, and timestep limitations. These issues are discussed at the close of the next section.

2.1.5 Multiscale methods in ultrafast DNA sequencing

In this section, we examine models for soft materials by applying multiscale methods to the problem of ultrafast sequencing of DNA through electronic means, which is a concept being vigorously pursued by several experimental groups (Kasianowicz et al., 1996; Meller et al., 2000; Li et al., 2003; Storm et al., 2005; Dekker, 2007; Fyta et al., 2011). The idea is to form a small pore, roughly of the same diameter as the DNA double helix, and detect the DNA sequence of bases by measuring the tunneling current across two electrodes at the edges of the pore during the translocation of the biopolymer from one side of the pore to the other. Estimates of the rate at which this sequencing can be done are on the order of 10 kb per second, which translates to sequencing the entire human genome in a couple of days. While this goal has not been achieved yet, the successful demonstration of this concept has the potential to produce long-lasting changes in the way medicine is practiced.

This system is also complex, involving several scales that are shown schematically in Figure 6. At the coarsest scale, the biopolymer in solution needs to be directed toward the nanopore. The scale of the biopolymer in solution is on the order of micrometers, and the timescale involved in finding the pore is on the order of fractions of a second. At these scales, the system can be reasonably modeled as consisting of a continuous polymeric chain in a uniform solvent under the influence of an external non-uniform but continuous field that drives it toward the pore. Understanding the translocation process involves modeling the polymer at scales set by the persistence length ( ~ 50 nm for DNA), while including the effect of the solvent. The coupling of the molecular motion to the fluid motion has been successfully demonstrated, with the solvent described by the lattice-Boltzmann approach to modeling fluid dynamics (Benzi et al., 1992; Succi, 2001) and the polymer motion described by Newtonian dynamics under the influence of local friction from the fluid motion (Fyta et al., 2006). This lattice-Boltzmann plus molecular-dynamics (LBMD) description has successfully reproduced details of the observed translocation process in a realistic manner, as observed in experiments for both single-file and multi-file translocation. An extension of the LBMD method to blood flow in the main arteries of the heart, including the explicit dynamics of red blood cells, was able to capture both the effects of non-Newtonian blood flow on the arterial walls (the “endothelial shear stress” linked to heart disease) and the formation of correlated cell motion (Peters et al., 2010). The coupling to the atomistic and electronic scales, required to capture the whole process of electronic sequencing during translocation, remains a significant challenge. The electron current can be calculated from quantum-transport theories, once the relevant electronic states are known. To this end, detailed calculations of the electronic properties of DNA bases need to be carried out, at the level of quantum chemistry, including all the details of the atomistic structure dictated by the morphology during the translocation process. Assembling the entire picture in a coherent model is a serious challenge to current computational capabilities and is likely to require significant conceptual and algorithmic innovations in the way we handle the phase-space of complex systems.

OPEN IN VIEWER

Figure 6.

Illustration of the length and timescales in ultrafast DNA sequencing by electronic signals during translocation through nanopores. Left: Schematic representation of the process and the various scales involved. Right: multiscale simulation coupling the molecular scale of DNA, represented by the blue spheres at the molecular persistence-length scale, to the motion of the fluid solvent whose local velocity is given by the colored contours.

The examples outlined in this and in the previous section illustrate several of the issues that arise in the context of multiscale simulations for hard and soft materials. The physics coupled in both cases is the motion of atoms and electrons, at the finest scale, as dictated by boundary conditions on the system at macroscopic scales, where continuum theories can be applied. The two ends of the spatial scale are intimately coupled, with the atomic motion constrained by the imposed boundary conditions and, in turn, determining the response of the system to the externally applied stress. The information passed between the scales is essentially the forces on the key components at each scale, which are derivatives of the appropriate energy functionals. These functionals have explicit dependence on degrees of freedom that are coupled across the scales.

Two key sources of errors and instabilities in these multiscale simulations are the differences in the energy functionals across two scales and the loss of information in passing from the finer to the coarser scale. The differences in energy functionals lead to the appearance of “ghost forces” since the physical quantities that underlie each functional may have a range much larger than the boundary region between the two scales. A case in point is the range of electronic wavefunctions, which can extend over much larger distances than what is practical to capture with a multiscale method involving a quantum-mechanical description at the finest length scale (for metallic solids, these wavefunctions may extend over the entire range of the solid). The second limitation, the loss of information due to coarsening of degrees of freedom, is intrinsic to the multiscale problem and can be addressed to some extent by performing finite-temperature simulations at the atomistic scale. These, however, are both extremely computationally demanding and restricted to a single temperature value. A more general framework is required to overcome this limitation.

The most severe limitation of multiscale methods that include atomistic or electronic scales is that the timestep is dictated by the dynamics at the smallest scale, which for atomic motion is on the order of 10 − 15 s. Interesting methodologies for extending the timescale have been proposed (see e.g., Voter, 1997a,b), but they have their own limitations in the size of system that can be treated or the range by which the timescale can be extended. This issue remains one of the most challenging theoretical problems in performing multiscale simulations for materials problems that can be of direct relevance to real applications.

2.1.6 Magnetic confinement fusion

Magnetic fusion is a long-term solution for producing electrical power for the world, and the large international device ITER (2012) being constructed will produce net energy and a path to fusion energy, provided the computer modeling is accurate. Fusion computation is enormously difficult because of the wide range of temporal and spatial scales in fusion devices (Post et al., 2004; Tang et al., 2009). Consequently, the fusion community has, over decades, developed a suite of models for tokamaks (and other confinement devices).

A model that quantitatively accounts for all the interacting physical processes at work in a tokamak discharge has come to be known as a whole device model (WDM). Such a model must account for nonlinearly interacting phenomena, including sources of energy and particles into the plasma; transport of energy, particles, and momentum through the plasma; the interface between the hot central plasma and the cool edge plasma and the material walls; and the interaction of the plasma with external magnetic coils, circuits, and control systems. While work of this kind has been going on in fusion research for a long time, limitations of theory and computer power traditionally have restricted such modeling to much-simplified models for the physical processes. Including this physics in a single, integrated, multiphysics, multiscale simulation presents many challenges of physics, applied mathematics, and computer science. Beyond the mathematics and computer science required to realize the physics components in an integrated simulation, there are important mathematical issues related to model coupling and stable, accurate timestepping algorithms, as well as computer science issues of data communication and efficient management of very large-scale computer resources.

As an example of one significant WDM project, FACETS, the Framework Application for Core-Edge Transport Simulations (Cary et al., 2009, 2011; Hakim et al., 2012), has the goal of providing modeling of a fusion device from the core to the wall. Figure 7 shows the cross-section of a tokamak plasma, which has a toroidal (donut) shape. The separatrix (so labeled) separates hot plasma in the center from the cooler, surrounding plasma. The central plasma is hot because it can move primarily only along the magnetic field lines, and those always remain in the core of the system. Thus, the plasma is never in direct contact with material walls, which would cool the plasma. In the edge region of the plasma, just outside the separatrix, the magnetic field lines (projections on a poloidal plane shown in green) wrap around the core but then terminate on the walls. This contact with the walls leads to much cooler plasma.

OPEN IN VIEWER

Figure 7.

Schematic diagram of a tokamak cross-section. The thick black line shows the material wall. The region inside the separatrix (closed brown contour) is the core, while outside is the edge. The physics in the core, edge, and material walls occur at different timescales and include different physical phenomena. Coupled simulations are needed to self-consistently understand the physics of the tokamak in order to optimize the fusion gain.

In the two regions, the plasma is analyzed with different approximations. In the core, the velocity distributions of the particles are nearly Maxwellian, and bulk parameters, such as temperature and density, are nearly constant along field lines, and hence, on flux surfaces, the surfaces on which the field lines remain. This constancy in flux surfaces implies that the plasma only has one-dimensional variation, on the dimension that crosses flux surfaces. In the edge, the plasma varies strongly along a field line, from the region in the midplane near the core to the ends, which are in contact with the walls. Hence, in this region, the plasma has two-dimensional variation. Being cooler, the edge region can have a high density of atomic species, whereas in the core, the plasma is for the most part fully ionized.

Just inside the separatrix is a transition region, where the cool, two-dimensional plasma of the edge morphs into the hot, one-dimensional plasma of the core. The atomic physics processes become less important, and the character of the turbulence, which controls the plasma confinement, changes. In this region a pedestal may form. At the pedestal there can almost be a transport barrier, such that the temperature of the plasma can be very different across it. Good pedestals are needed for the beneficial H-mode of operation.

Important to the dynamics of this system is how the edge plasma interacts with the wall. The fuel for fusion is hydrogen (actually, its isotopes, deuterium and tritium). The nuclei of these atoms can diffuse into the wall and be released at a later time. These processes are important to fusion confinement devices because they act as boundary conditions on the edge plasma, determining, in part, the influx of particles.

FACETS uses a component-based approach to couple the various physics. Currently, the core, edge, and wall equations are evolved separately and coupled with an explicit coupling scheme; preliminary work on implicit coupling is promising. The coupling between the core and the edge involves a transformation of the two-dimensional edge data to a boundary condition for the one-dimensional core equations, while a reverse transformation is required to transfer fluxes from the core to the edge. The edge–wall coupling is a pointwise coupling in which each cell on the edge boundary is coupled to a wall component.

Coupled whole-device simulations of tokamaks are required in order to understand and optimize fusion gain. The current approach to coupling is essentially code-to-code coupling between existing physics components that were not originally designed for this purpose. Hence, the challenge is not only to retrofit these components to work with one another but also to design robust and accurate algorithms that evolve the coupled physics. Further, developing new WDM components that take into account emerging high-performance architectures is a significant challenge that needs to be addressed in order to satisfy the modeling requirements as ITER and other tokamaks become operational.

2.1.7 Subsurface science

Problems in hydrogeology encompass a broad range of multiphysics and multiscale issues. In this section, we discuss several solution and coupling approaches in the area of subsurface fluid flow, solute transport, and reactions.

First, we consider continuum-scale subsurface reactive transport. For several decades researchers have studied the coupling of fluid flow, multicomponent solute transport, and geochemical reactions. In the 1980s, researchers (Rubin, 1983; Cederberg et al., 1985; Lichtner, 1985) developed formulations for coupling multispecies chemistry with multicomponent transport. Yeh and Tripathi (1989) and Steefel and MacQuarrie (1996) provide overviews of these techniques, which can be divided into three general categories: (1) the global implicit approach (GIA), where solute transport and geochemical reaction are tightly coupled and solved simultaneously in a single system of nonlinear equations; (2) the sequential iteration approach (SIA), where solutes are transported independently, then coupled through geochemical reaction, with the process iterating until the global residual equations converge to a prescribed tolerance; and (3) the sequential non-iterative approach (SNIA), where the solutes are transported independently and coupled through geochemical reaction, but with the process not iterated to convergence. The SNIA approach is more commonly referred to as “operator splitting.” Steefel and MacQuarrie (1996) present a concise overview of these three approaches and provide examples illustrating the trade-offs of each.

Through the years, researchers have argued in favor of SNIA or SIA over the fully coupled GIA because of GIA’s excessive computing requirements (Yeh and Tripathi, 1989). However, both SIA and SNIA, with less coupling, present their own challenges regarding accuracy and performance (Valocchi and Malmstead, 1992; Barry et al., 1996; Kanney et al., 2003a,b; Carrayrou et al., 2004). SNIA is known to introduce operator-splitting error in the solution, and SIA may result in prohibitive iteration counts when geochemical species are tightly coupled to the immobile phase (e.g., strong sorption, mineral precipitation–dissolution with fast reaction kinetics). An interesting modification can be made to SIA by selectively coupling a subset of strongly coupled reacting species and solving the resulting subsystem of equations with GIA within the SIA framework (Robinson et al., 2000). This selective coupling can reduce SIA iteration counts by two orders of magnitude and model runtimes by more than one order of magnitude.

Recent advances in high-performance computing (HPC) have enabled the use of GIA to solve three-dimensional problems comprising large numbers of chemical species. The feasibility of GIA has been demonstrated on a groundwater-flow and reactive-transport problem involving tens of millions of degrees of freedom simulated on thousands of processor cores (Hammond and Lichtner, 2010). With the use of increasing numbers of processor cores (e.g., 10,000 or more), simulating these larger nonlinear systems of equations requires the development of new solver algorithms to preserve nonlinear and linear iterative solver iteration counts and to minimize communication requirements. Physics-based preconditioners have been investigated for reactive transport (Hammond et al., 2005), where the GIA Jacobian operator is split into two separate operators for global transport and local reaction. These two transport and reaction operators are then applied sequentially in the preconditioning step, much in the same way that operator splitting is applied as a solution technique. This physics-based preconditioner performs well for subsurface problems composed of primarily aqueous reactions, but more poorly for scenarios where one would expect larger operator-splitting error (e.g., fast reaction kinetics and solute retardation through sorption and mineral precipitation).

Next, we go beyond the continuum scale and consider multiscale approaches that couple pore-scale to continuum-scale models. At the pore scale, the geometry of solid grains and pore spaces is explicitly defined, and the physics represented are Navier–Stokes equations for single-phase fluid flow, advection, and diffusion as transport processes, and various aqueous and heterogeneous reactions (both equilibrium and kinetic). Pore-scale models take various forms, including traditional mesh-based computational fluid dynamics (CFD) methods, lattice-Boltzmann methods, pore network models, and fully Lagrangian smoothed-particle-hydrodynamics (SPH) particle methods (Gingold and Monagham, 1977).

In this multiscale approach, at the continuum scale, the system is treated as a “porous medium” without consideration of explicit grains and pores. Continuum-scale physics includes conservation of mass combined with Darcy’s law for single-phase fluid flow, mean advection, and dispersion (including diffusion and effects of unresolved velocity variations) as transport processes, and a similar set of chemical reactions. Most continuum-scale models use traditional numerical methods (finite difference, finite element, finite volume) to discretize and solve the system of PDEs, as discussed above.

Since efforts in this treatment of multiscale–multiphysics coupling in subsurface reactive-transport problems have been initiated only recently, the literature contains only a small number of examples. These use both concurrent and hierarchical approaches. In concurrent methods, different physics are solved on separate subdomains with exchange of information and boundary condition matching at the boundaries between the subdomains (or in some cases on an overlapping region). Balhoff et al. (2007) used a concurrent approach to couple a pore network model of flow within a sand-filled fracture to a continuum model of flow in the surrounding porous matrix. Tartakovsky et al. (2008) have developed a coupled pore-continuum model of a diffusion/reaction system in which SPH is used at both scales; this allows straightforward model coupling at the domain boundary but is difficult to adapt to situations in which advective flow exists. Mortar methods (Balhoff et al., 2008) also provide a means of coupling subdomains with different physics by solving a boundary-matching problem on the mortar grid. Battiato et al. (2011) used an iterative boundary-matching scheme to couple pore- and continuum-scale models of reactive transport including advective flow. Alternatively, hierarchical methods solve both models simultaneously over common domains and pass information between scales using various approaches. Tartakovsky and Scheibe (2011) have successfully applied a hierarchical method called dimension reduction with numerical closure to the simulation of mixing-controlled mineral precipitation. In this approach, conceptually similar to the heterogeneous multiscale method of E et al. (2003), macroscopic properties and states are defined in terms of pore-scale variables, and short bursts of pore-scale simulation are used to provide numerical closure for a given set of local conditions. Then, a larger timestep is taken with the macroscopic model (eliminating a significant amount of pore-scale simulation) before it is necessary to update the numerical closure to reflect new macroscale conditions.

2.1.8 Surface and subsurface hydrology

A blueprint for modeling fully integrated surface, subsurface, and land-surface processes, originally proposed 40 years ago (Freeze and Harlan, 1969), is now becoming a reality. Although truly coupled models have appeared only recently in the literature (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004; Jones et al., 2006; Kollet and Maxwell, 2006; Qu and Duffy, 2007; Kollet and Maxwell, 2008), a growing library of models and community of modelers now contribute considerably to our understanding of the coupled terrestrial hydrologic and energy cycles. These coupled and integrated hydrologic models are being used on a range of important issues. The models have common features in that they solve some form of the nonlinear diffusion equation (e.g., the Richards equation) for subsurface water flow and some form of the shallow water equations (e.g., kinematic or diffusive wave) for surface flow. Additionally, these hydrologic models are being coupled to land-surface energy, biogeochemistry/ecology models and to atmospheric models (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004; Maxwell and Miller, 2005; Kollet and Maxwell, 2006; Maxwell et al., 2007; Kollet and Maxwell, 2008; Maxwell et al., 2011). Updated numerical and computational technologies have, in part, enabled new approaches for modeling these coupled interactions. While these models all take different numerical, discretization, and coupling approaches, they share the common goal of modeling the terrestrial hydrologic and energy cycle as an integrated system, rigorously and mathematically. Research addressing these issues encompasses a range of scales and includes a variety of processes (Ebel et al., 2009). Hydrologic simulation models break into two broad solution categories: a globally implicit formulation, and operator-splitting or Gauss–Seidel methods. Among the globally implicit hydrologic models, there are two distinct formulations for the solution of surface and subsurface flow (Ebel et al., 2009): first-order exchange (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004) and continuity of pressure (Kollet and Maxwell, 2006). Integrated hydrologic models solve equations for surface flow (generally the shallow wave equations) and subsurface flow (generally the Richards equation) simultaneously (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004; Jones et al., 2006; Kollet and Maxwell, 2006; Qu and Duffy, 2007; Kollet and Maxwell, 2008; Camporese et al., 2009; Kumar et al., 2009; Park et al., 2009b) and represent both coupling approaches and a range of numerical implementations of the solution (finite-element methods and finite-difference/finite-volume methods). While not widespread, integrated hydrologic models have been designed to run on large-scale computing infrastructures and, in one case, demonstrate very good weak scaling out to large numbers of processors (Kollet et al., 2010).

ParFlow solves the variably saturated Richards equation in three dimensions. Overland flow and groundwater surface-water interactions are represented through a free-surface overland-flow boundary condition, which routes ponded water as via the kinematic wave equation (Kollet and Maxwell, 2006). The coupled subsurface–overland equations are solved simultaneously using a highly efficient and robust Newton–Krylov method with multigrid preconditioning, which exhibits excellent parallel scaling and allows the simulation of large-scale, highly heterogeneous problems (Ashby and Falgout, 1996; Jones and Woodward, 2001). Simultaneous solution ensures full interaction between surface and subsurface flows; this is particularly important in regions with spring-fed streams, large baseflows, and rivers that alternate between gaining and losing.

As one example of multiscale hydrological modeling, ParFlow solves the coupled water and energy balance at the land surface by incorporating a modified version of the Common Land Model (CLM3.0) (Dai et al., 2003), which calculates evaporation from the vegetation canopy and ground surface; transpiration and plant water-uptake over the soil column; snow accumulation and melting; and latent, sensible, and ground heat fluxes. Vegetation and land-energy processes are calculated at each timestep by using an operator-splitting approach that ensures full mass and energy conservation (Maxwell and Miller, 2005; Kollet and Maxwell, 2008). A fully coupled hydrologic and climate model, PF.WRF was developed by combining ParFlow with WRF (Maxwell et al., 2011). The general solution procedure begins with explicit advancement of the atmospheric code. Using the operator-splitting approach, surface fluxes that are relevant to ParFlow, such as infiltration and evaporation rates, are integrated within the atmospheric model over the ParFlow timestep. These accumulated surface fluxes are applied to ParFlow at the start of its implicit time advancement. The subsurface moisture field calculated by ParFlow is passed directly to the Noah land-surface model (Mitchell, 2005) within WRF for use in the next timestep. Noah is advanced for each WRF timestep to provide the necessary land-surface fluxes, but soil moisture values are now specified by ParFlow. ParFlow is incorporated as a subroutine in Noah, with communication over all Noah soil layers; subsurface hydrology in Noah is entirely replaced by ParFlow. The combined model is fully mass and energy conservative.

2.1.9 Climate modeling

Earth system modeling, specifically coupled global climate models (GCMs), is one of the truly definitive multiphysics applications in several respects. GCMs encompass physical behavior ranging in spatial scales from molecular dynamics to global atmospheric wave patterns. Moreover, GCMs capture climate variations covering months to 10,000 years over regional to global areas of the Earth’s surface. A number of GCMs are under active development around the world, and the Community Earth System Model (CESM), which has been developed by hundreds of scientists and contributes to the Intergovernmental Panel of Climate Change Assessment Report (Solomon et al., 2007), is a good example to illustrate the multiphysics nature of GCM. In most cases, the multiphysics characteristics of the CESM apply to other GCMs.

The CESM comprises five Earth system components (atmosphere, land surface, ocean, sea ice, and ice sheets) and has provision for future components, as shown in Figure 8. These components periodically exchange data as fluxes at component boundaries through a flux coupler, which coordinates the transfer of coupled physics across model components, except for ice sheets, which currently receive one-way surface–mass balance values from the land model. Typically, coupling occurs every simulated day, although recent high-resolution coupled model development has indicated the necessity of exchanging this information more frequently (Dennis et al., 2012). The source-model flux data is interpolated to the target component by using second-order or conserved mapping weights within the coupler. Coastal “masks” which account for the boundary between land surface and ocean are created as a data file to determine whether atmosphere model grid cells couple to the land surface or ocean below and vice versa. These masks need to be created for each combination of component grids used in a coupled-model run.

OPEN IN VIEWER

Figure 8.

A schematic diagram illustrating the multi-fidelity structure of the Earth-system model. Except for the ice sheet model, all components in the model mutually couple through the flux coupler, which is shown as the hub of the diagram. The ice sheet model does not couple in the same way. It is coupled only to the land model, receiving surface-mass balance values. As is shown, this coupling is one way: the land model receives no feedback from the ice sheet model.

Within each model component, resolved physics are represented by discretized equations referred to as the dynamical cores, or dycores. Currently, all the physics and dynamics variables exist on the scale of the three-dimensional discretized grid, although the land model performs computations using statistical representations of subgrid-scale land-use characteristics. The atmospheric component in CESM, the Community Atmosphere Model (CAM), has several dycore options that scientists can choose based on desired attributes such as computational efficiency and local conservation properties. Each model component uses a different strategy for solving the model equations. For details, the reader is referred to the model documentation (Gent et al., 2011). Unlike CAM, which has several choices for spatial discretization, the ice sheets, sea-ice, and ocean components in CESM use a finite-difference scheme on regularly spaced grids. Their time-integration methods are diverse, arising from the different equations and physical scales being treated in each component. These methods include subcycled fully explicit (sea-ice), semi-implicit to “step-over” barotropic dynamics (ocean), and fully implicit velocity with subcycled advection (ice-sheets) methods. Several of these components have developmental dycores using an unstructured grid to allow local refinement (Ringler et al., 2011) and a testbed of fully implicit algorithms to overcome stability limitations (Evans et al., 2009).

As mentioned above, CAM has several dycore options for scientists to select when configuring the coupled climate model for experiments. These options include a spectral-element-discretized cubed-sphere configuration, which uses a subcycled fully explicit numerical scheme. Other dycore options include an Eulerian spectral and spherical finite-volume discretization, which both use an operator-splitting, semi-implicit time-integration method to separate waves of different scales. In this case, a common splitting of advection and gravity waves adapted from the geophysical fluid dynamics community is employed because it takes advantage of a natural separation at coarse resolution, and favorable numerical properties of the resulting equations. The atmospheric dynamical core resolves horizontal waves ranging from vertically propagating gravity waves ( ~ 200 m/s) to low-frequency global modes within the background flow ( ~ 2–20 m/s). Vertical motions are smaller in scale and are represented by the hydrostatic approximation, which assumes that the vertical pressure gradient is largely balanced by gravity and is valid to scales down to individual cloud systems of about 100 km in extent.

Climate is uniquely challenging because the most important drivers of global climate change (chemistry, moisture, and radiation) vary over microscopic spatial scales but impact the global-scale net energy balance of the Earth system. The global-scale atmospheric flow is coupled to local chemical, cloud, and radiative forcings; and the time scale for this coupling ranges from the diurnal cycle to decades. Subgrid-scale effects are treated either as prescribed values from datasets or as separate, vertical-column, horizontally independent models that are updated as forcings to the dynamical equations. These independent models are themselves multiphysical. For example, the chemistry model in CESM (Horowitz et al., 2003) uses a split explicit–implicit solution method for the slow and fast species reactions, respectively.

2.1.10 Radiation hydrodynamics

Physical processes that involve radiation or hydrodynamics individually are part of many phenomena of everyday life and everyday scientific research. However, phenomena where radiation and hydrodynamics couple and stand together on an equal footing, where strong coupling between the two components affects the evolution of each, rarely occur terrestrially but are common astrophysically. The energetics involved explain this apparent dichotomy. In astrophysics, RHD plays a part in highly energetic phenomena, such as stellar core collapse, supernova shock propagation, supernova remnants, stellar pulsations, and accretion flow. Where radiation-hydrodynamic phenomena occur terrestrially, they are restricted to high-energy density environments: inertial-confinement fusion, laboratory astrophysics, and high-temperature shock-wave and explosive phenomena.

In addition to RHD being a multiphysics field, in that it strongly couples the evolution of distinct physical components (matter and radiation), many problems in RHD introduce other multiphysics aspects. For example, a single problem can involve any number of the four fundamental interactions across multiphysics domains and can mix classical- and modern-physics paradigms (e.g., classical or relativistic space–time and gravity, quantum statistics, electroweak, and strong nuclear interactions).

Numerical RHD models are typified by a number of common features (Pomraning, 1973; Mihalas and Mihalas, 1984; Castor, 2004): conservation PDEs play a fundamental role (e.g., mass, species, energy, momentum). As a result of the energies involved, shock waves are frequently present in such systems, leading to discontinuities in physical quantities. Additionally, material interfaces, which are part of many models of terrestrial setups, introduce another form of discontinuity into the problem. Some RHD problems can be classified as reactive-flow problems, as evidenced by complex and strong coupling terms between components. Figure 9 shows an example of the complexity of these systems: interacting shock waves are propagating in a simulation of a terrestrial shock-tube experiment as part of the CRASH project (Drake et al., 2011). A beryllium foil is driven into a xenon-filled shock tube by a drive laser to speeds in excess of 100 km/s. This process generates a radiative primary shock wave, which in turn generates a wall shock as the plastic walls of the shock tube are heated and ablated.

OPEN IN VIEWER

Figure 9.

Interacting shock waves propagating in a simulation of a high-energy-density physics shock-tube experiment (Drake et al., 2011). The figure shows a slice representing a half cross-section of the cylindrical shock tube with the central axis at the bottom of the figure. In addition to the shock fronts, a number of material interfaces are shown, which represent boundaries between the multiple materials used in the experiment (beryllium and xenon within a polyimide-plastic shock tube). The horizontal extent of this schematic represents roughly 0.6 mm of a shock tube that is about 4.0 mm in total length.

Other complexities also are common in RHD models. Frequently, the equations for the radiation are nonlinear, requiring solution by nested iterative procedures, such as preconditioned Newton–Krylov methods. (These nonlinearities can be caused by strong coupling between matter and radiation as well as by quantum-mechanical degeneracy effects.) Coupling between physical components can be a sensitive issue because, in complex phenomena, coupling varies in strength as a function of time, position, and/or energy and thus can require different solution approaches in different situations. High dimensionality is also an issue. Proper treatment of the radiation requires solution of the Boltzmann transport equation over six-dimensional phase space as a function of time. Almost always, finite computer resources force a compromise, leading to reduced dimensionality in some fashion (Graziani, 2005). Such compromises can affect self-consistency of coupling (e.g., radiation–matter momentum exchange in flux-limited diffusion approximations).

Since numerical solution of RHD problems predates the computer age, solution techniques are many and varied. Today, however, solutions typically involve the solution of Euler’s equations for the matter component. (Reynolds numbers are typically very high in such systems.) Standard methods include finite-difference and finite-volume methods. Lagrangian, Eulerian, and hybrid approaches (e.g., ALE methods) are used. Such techniques also translate to solution of the radiation (Bowers and Wilson, 1982; Bruenn, 1985; Turner and Stone, 2001; Swesty and Myra, 2009; Van der Holst et al., 2011). Here, however, finite-element methods also find use (Lewis and Miller, 1993).

When a similar discretization is used for both matter and radiation, coupling between components reduces to a problem of obtaining consistency and accuracy of paired coupling terms. However, when different discretizations are used for matter and radiation (e.g., finite volume for matter and finite element for radiation), then accuracy of the transfer operator can become problematic. Unlike frequently encountered multiphysics couplings discussed in this report, matter–radiation coupling occurs, in some strength, over the entire physical domain of interest; it is a volumetric coupling. Methods for improving the accuracy of this coupling are an area of active research.

RHD models commonly use a mixture of explicit and implicit approaches within the context of operator splitting. Choices stem from the timescales involved. Typically, hydrodynamic timescales are much greater than radiative timescales. This situation leads naturally to use of explicit methods for the matter and implicit methods for the radiation (usually sparse iterative approaches for the latter). However, for problems with long evolutionary timescales (e.g., neutrino emission for cooling neutron stars (Reynolds et al., 2008)), implicit approaches increasingly are being used for both components.

The difficulties inherent in adequate treatment of multiphysics coupling in RHD models are illustrated by a sample list of the myriad of relevant timescales involved: acoustic, radiation-transport and radiation-source, electron–ion equilibration, electron diffusion, and viscous (if present). Arrayed in relative size, these may span a range of 10 15 . A challenge for the application modeler is to determine which of the associated processes are to be followed in detail and which can be assumed to occur so quickly that an equilibrium approximation is adequate.

2.1.11 Geodynamics

Geodynamics is the study of dynamic processes in the Earth from crust to core. Processes range from earthquakes to plate motions, large-scale convection to microscale melt migration. Physical coupling occurs over a wide range of simulations in geodynamics. Short-term crustal behavior is generally modeled with some combination of a porous, viscous, elastic, or plastic material model (Watts, 2001). Beneath the crust, incompressible mantle flow calculations couple the advection–diffusion of energy/temperature to the Stokes equations (Schubert et al., 2001). The complexity may be increased by introducing the advection of chemical components with low or zero diffusivity, which is often achieved by using particle-based methods (Van Keken et al., 1997). Compressible flow, which normally uses an elastic approximation to filter out seismic waves and timescales, may also introduce reference equations of state derived from models of seismic wave propagation (King et al., 2010). These simulations are often multiscale, including narrow plume conduits in large domains (Tarduno et al., 2009; Davies et al., 2011), along with physical discontinuities in material properties at a variety of depths (Christensen and Yuen, 1985). At greater depths, in the outer core, models couple turbulent flow from the Navier–Stokes equations to thermal and magnetic models of rapidly convecting liquid iron (Glatzmaier and Roberts, 1995; Kuang and Bloxham, 1999).

The change of behavior and physical timescales with depth generally leads to different models being used for different shells of the Earth. Boundary conditions are then used to parameterize the fluxes between models. As an example, owing to high convective vigor in the outer core, the lower thermal boundary condition on the mantle may be treated as a Dirichlet condition (Blankenbach et al., 1989), while the motions of crustal plates may be treated as a rigid shell, with prescribed shear at the upper boundary of the mantle (Bunge et al., 1998). Depending on the application, this outer boundary may also be treated as a free surface, supporting dynamic topography, which involves coupling the solution of a kinematic boundary condition to the flow solution in the global domain. This approach potentially introduces instabilities in incompressible flows (Kaus et al., 2010), as shown in Figure 10. Smaller-scale, higher-resolution studies may occur on overlapping domains. As an example, one could focus on the production of melt in subcrustal regions, such as mid-ocean ridges or subduction zones. The differing temporal and spatial scales of pore-filling fluid and melt lead to calculations that are often multiphase, coupling Darcy’s flow with Stokes flow in the mantle or poro-visco-elastic flow in the crust (Spiegelman, 1993a,b; Katz et al., 2007). Experimental, theoretical, and numerical work attempts to link these high-resolution studies to larger domains through parameterizations of rheology and melting (Katz et al., 2003).

OPEN IN VIEWER

Figure 10.

The “drunken sailor” instability. In a simple case of a compositionally (grey shading) driven Rayleigh–Taylor instability, an explicit coupling between the incompressible Stokes equations and the kinematic free-surface on the top boundary leads to a sloshing instability in the velocity field (vectors) (Kaus et al., 2010).

2.1.12 Particle accelerator design

A key component technology in particle accelerators is the accelerator cavity, which is used to generate electromagnetic fields that control acceleration of the subatomic particles in the accelerator beam. To achieve high efficiency, a modern accelerator cavity uses superconducting devices to generate the fields that drive the cavity at its resonant frequency.

In radio-frequency (RF) cavity design, finding a sufficiently accurate resonant frequency for real geometry has become possible with the help of a higher-order, high-fidelity finite-element method (Morgan, 1996; Gupta et al., 1997; Jin, 2002; Sun et al., 2002). One of the main sources of operational failures of RF cavities is strongly related to the excessive thermal heating arising from high power or high currents in operation (Reece et al., 2007; Cheng et al., 2011). As the power and current increase for future experiments, optimizing the performance and cost-effectiveness of RF cavities becomes an increasingly important technological challenge for designing next-generation accelerators. TEM3P, a multiphysics analysis tool using finite-element methods, has been developed to capture thermal, structural, and electromagnetic effects such as cavity wall heating (Akcelik et al., 2009), structural deformations (Akcelik et al., 2008a), and Lorentz force detuning (Liepe et al., 2001). Parallel implementation of TEM3P utilizes supercomputers for large-scale computations (Lee et al., 2004; Liao et al., 2006; Li et al., 2007) with high geometric fidelity and high numerical accuracy. TEM3P solves the thermoelastic equation, in which the heat flux and radiation pressure are interpolated from the electromagnetic solution file from Omega3P (Lee et al., 2009), a parallel finite-element code for frequency domain analysis of accelerator cavities.

The design process starts with accurate electromagnetic solutions obtained from Omega3P analysis. Omega3P solves Maxwell’s equations with higher-order curved elements for highly accurate solutions. TEM3P reads the electromagnetic solution file (in NCDF format), which is interpolated to produce accurate physical quantities of interest for the given multiphysics application. Heat flux is calculated for the thermal solver in TEM3P (Akcelik et al., 2009), while radiation pressure is calculated for the structural solver in TEM3P (Akcelik et al., 2008a). Interpolation of the electromagnetic fields accounts for the general case, in which the nodal points in the TEM3P mesh do not coincide with the nodal points in the Omega3P mesh. Once the desired heat flux or radiation pressure is obtained, it is applied as a boundary condition to the existing thermal model for an RF heating problem or to the existing structural model for a Lorentz Force detuning problem.

Internally, TEM3P couples its thermal solution to its structural solver. When a temperature distribution is known, thermal stresses due to the high-temperature gradient can be readily calculated from the structural solver. As shown in Figure 11b, the high-temperature gradient in the RF gun is strong enough to cause geometric distortions so that the resulting frequency shifts can be larger than the operating frequency bandwidth (Akcelik et al., 2008a).

OPEN IN VIEWER

Figure 11.

Multiphysics coupling analysis of a Linac Coherent Light Source (LCLS) gun.

Commonly, in multiphysics analysis, the primary source of heat flux is RF heating. In this case, an electro-thermo-structural coupling, as shown in Figure 12 is used to provide viable solutions that address frequency shift problems. TEM3P has also been applied to thermal instability problems in superconducting cavities (Akcelik et al., 2009). Thermal instability can lead to serious problems such as thermal quenching, which terminates the superconducting state of a cavity (Kim, 2003). In the face of this instability, operation is usually not possible, or only at a significantly reduced efficiency if it is possible (Lesrel et al., 1997). Superconductor materials exhibit highly nonlinear thermal characteristics, and nonlinear boundary conditions need to be applied in order to obtain an accurate solution. A robust nonlinear solver is implemented in TEM3P to address these problems arising from the nonlinear material properties.

OPEN IN VIEWER

Figure 12.

Computational domain (left) and pressure and velocity contours (right) for coupling of Navier–Stokes flow on the surface and Darcy flow in subsurface porous media.

The current version of TEM3P provides steady-state solutions for both thermal and structural problems. Multiphysics problems with multiple spatial and temporal scales are not yet addressed in TEM3P. As designs become more concerned with fully time-dependent multiphysics coupling, further research and development become increasingly important.

Shape optimization studies have been performed with adjoint gradient information (inverse eigenvalue problems) (Akcelik et al., 2005, 2008b). Because adjoint variables enable gradient calculation at a given design point regardless of the number of design variables, achieving a near-optimal shape in a reasonable turnaround time has become feasible (Choi et al., 2008; Lee, 2010; Xiao et al., 2011). With the successful implementations of multiphysics coupling, applying its solution to optimization problems could open up new research areas. The solutions from the multiphysics coupling could provide insights in the complex physical interactions not resolvable from single-physics analysis. Thus, combining an adjoint-based gradient optimization with multiphysics coupling could potentially play an important role in achieving near-optimal solutions in multiphysics environments.

2.2.1 Choices and challenges in coupling algorithms

Spatial interfaces remain a primary challenge, especially between subdomains using different computational approaches and discretizations (e.g., an interface between two subdomains, one discretized by finite elements and the other by finite volumes). Many commonly used interpolation approaches are non-commutative and thus can lead to non-conservation of supposedly conserved quantities. A number of applications we have described use mortar methods to bridge this difficulty (see Sections 2.1.1, 2.1.2, and 2.1.7). Methods for coupling multiphysics components in space are discussed in Section 3.2.1.

The collection of applications described in Section 2.1 suggests that explicit methods are still prevalent across many fields, especially for “production” applications. This situation is understandable: compared with implicit approaches, explicit methods are easier to implement, maintain, and enhance from a software-engineering perspective, as discussed in Section 4.1.1. Maintaining the stability of an explicit method mandates that the programmer understand at least some the important timescales involved in the problem. Implicit methods are free of splitting errors that are introduced by operator splitting, and may allow the use of much larger timesteps than an explicit method. However, when using an implicit method, timesteps must be carefully selected to control at least local time-integration error in order to avoid “silently” producing inaccurate answers by potentially taking large timesteps over important processes having fast timescales. This is of special concern when, as occurs frequently, the equations being solved are nonlinear. (Nonlinearities can occur in nearly all the application examples above, but notably those described in Sections 2.1.1, 2.1.3, 2.1.6, 2.1.7, 2.1.8, 2.1.9, 2.1.10, and 2.1.12.) In summary, the explicit–implicit debate becomes one of accuracy versus efficiency. Given the increasing use of implicit methods, it is necessary to ask how we should approach judging the accuracy of a method. For example, it is easy to integrate y ′ ( t ) = y ( t ) 2 from y ( 0 ) = 1 in a stable way that completely misses the singularity at t = 1 and continues smoothly with unbounded error into t > 1 . How do we know we are not skipping over physics? What is the collection of metrics we can use to provide confidence in our answers? Which integrated metrics are of greatest value? Which pointwise metrics are reliable?

A number of the applications discussed in this report have timescales for important physical processes that may vary by many orders of magnitude. (See Sections 2.1.4, 2.1.5, 2.1.9, 2.1.10, and 2.1.11.) The challenge for the model builder is to have a clear picture of what quantities need to be tracked in the problem being solved. Such insight helps to determine which of the associated processes need to be followed in detail and which can be assumed to occur too quickly for their evolutionary details to be important. An additional challenge is to write that intelligence into the code itself. Section 3.2.2 discusses methods for coupling multiphysics components in time, including implicit–explicit (IMEX) and adaptive techniques, which can help address such multiphysics challenges in temporal coupling. Additional opportunities are considered in Section 5.

Since diverse timescales are common in multiphysics applications, it is worth asking whether an approach using UQ methods may yield insight into the relation between timestep selection and results obtained; UQ issues are discussed in Section 3.5.

A related issue involves quantifying certain cultural knowledge that allows one to “cheat” a little on constraints when quantities of interest are not affected. An example is the careful violation of the Courant–Friedrichs–Levy (CFL) condition in explicit codes.

Another interesting challenge is spatial coupling of the discrete to the continuum. This issue arises in problems in which the macroscopic evolution of the system relies on detailed knowledge of phenomena occurring on a microscopic level. Examples include the modeling of crack propagation and the ultrafast sequencing of DNA discussed in Sections 2.1.4 and 2.1.5. An additional class of problems in this vein is non-equilibrium fluid flow: adequate description of state variables in the fluid may rely upon, for example, the result of a simulation at the molecular, atomic, or even electronic or nuclear scale. Problems in RHD (see Section 2.1.10) can fall into this class.

Discrete-to-continuum coupling may also arise at interfaces between diverse materials where, in one material, a microscopic model is required for a satisfactory physical description, while in the other material, a classical continuum description is adequate. Whenever these sorts of couplings arise, one or more additional mesoscale descriptions that link the two extremes has found use. (Section 2.1.4 describes an application that benefits from this approach.) We discuss discrete–continuum coupling more fully in Section 3.3.

Any need for adaptivity in algorithms throughout a given multiphysics simulation also presents challenges. As an example, consider a multiphysics application in which a key piece of science changes in the middle of the model’s evolution (e.g., a rate may become much larger or smaller as a result of some change in a linked quantity). This situation often leads to the component equations having a changed character. For example, what might have previously been a parabolic problem might now take on hyperbolic behavior. This circumstance may require a change of solution methods (introduced in Section 3.1), anywhere from the preconditioner to the whole solver approach. Equally extreme are phase transitions: solution methods are likely to be different for a subdomain now containing a solid phase instead of a gaseous phase. These issues, as well as related challenges in supporting changes in modeling and algorithms over the lifetimes of multiphysics applications, are considered from a software-design perspective in Section 4.

2.2.2 Software engineering

The current status of multiphysics software is discussed in Section 4.1, including current practices for multiphysics coupling in production codes, common needs, and examples of libraries, frameworks, and tools that support multiphysics simulations in HPC environments. The development of large-scale multiphysics applications raises many of the same issues as does development of any large-scale application. One of the most fundamental issues is the “buy” versus “build” choice. In the multiphysics world, however, this question takes on increased importance because of the added complexity. The “build” approach frequently requires large collaborations of diverse expertise. However, the sources where one may “buy” are often similarly scattered among disciplines and institutions. Neither approach alone guarantees robust, accurate code. Risk mitigation also plays a role in this choice. On the one hand, third-party code is not designed with the end-user’s specific purpose in mind. Learning the code (and perhaps learning that the code is not suitable) can be time-consuming. Even if such a third-party code is suitable, additional investment is required to maintain multiple, interacting software packages. A decision to “build” can be equally risky. In addition to the obvious issues of funding, assembling, and directing a collaborative team, and then having the team deliver on time, there is an ever-present danger that the drive for science leads to compromises in code accuracy, especially where the subtle coupling of multiphysics models occurs. Challenges, near-term needs, and risks in multiphysics software are discussed in Sections 4.2 and 4.3, while longer-term opportunities and insertion paths are considered in Sections 5 and 6.

2.2.3 Analysis and verification

Analysis and verification of coupled models involve issues that span both algorithms and software and have been only partially addressed by the community overall. Sections 3.4 and 3.5 discuss error estimation and UQ for multiphysics problems from an applied mathematics perspective, while Section 4.2 considers issues in software testing. Verification of coupled codes is challenging. Unit tests do not test the synergistic effects of coupling, which is the very aspect of simulation that probably requires the greatest testing. In addition, large multiphysics codes have too many coupling combinations to allow practical coverage of the code space by any conceivable test suite. Furthermore, standard reference problems are lacking. Given this situation, verification tests tend to devolve into regression tests, whose value further devolves when existing physics is changed or new physics is introduced.

3.1.1 Methods for systems of linear equations

Consider the system of linear equations

A u = f A ∈ R N d × N d u , f ∈ R N d

4

where N d is the number of degrees of freedom. Equation (4) can represent an entire multiphysics problem or some subproblem encountered in a multiphysics solution algorithm. Methods for solving (4) fall into two categories: direct and iterative methods. Direct methods generally rely on some canonical decomposition of A and scale poorly (in both storage and operations) when applied to large problems, even when sparsity in A is effectively exploited. For many years, researchers relied on stationary iterative methods, such as classical Jacobi, Gauss–Seidel, or successive overrelaxation methods. While effectively alleviating the large storage requirements of direct methods and being easy to implement, stationary methods lack robustness and generally do not converge fast enough to be considered practical for realistic multiphysics problems. Current best practices for the efficient solving of (4) are based almost exclusively on multigrid methods and preconditioned Krylov subspace methods.

Krylov subspace methods

Let u 0 ∈ R N d be an initial approximation to the solution of (4), let r 0 = f − A u 0 be the initial residual, and for each integer k ≥ 1 let

K k ( r 0 , A ) = s p a n { r 0 , A r 0 , … , A k − 1 r 0 }

be the Krylov subspace generated by r 0 and A . Consider first the case where A in (4) is symmetric and positive definite. In this case, A induces an inner product ⟨ u , v ⟩ A = u T A v and corresponding proper norm ∥ u ∥ A = ⟨ u , u ⟩ A . For symmetric positive definite A , the conjugate-gradient (CG) method (Hestenes and Stiefel, 1952) generates a sequence u k ∈ u 0 + K k ( r 0 , A ) of approximations to the solution u ⋆ of (4) such that ∥ u k − u ⋆ ∥ A is minimized over u 0 + K k ( r 0 , A ) and the residuals r k = f − A u k are orthogonal to K k ( r 0 , A ) . Because of these optimality properties, CG is guaranteed to converge in at most N d iterations in exact arithmetic. Remarkably, these optimality properties can be achieved by using simple, three-term recurrence relations and a small number of auxiliary vectors. Hence, CG was initially seen as a direct method with low storage that could compete with factorization-based direct solution methods. It was soon realized, however, that the key orthogonality relationships could not be maintained with finite precision, and interest in CG waned until the early 1970s. At that time, it was pointed out (Reid, 1971) that CG could still be used effectively as an iterative method. Ideas for accelerating convergence using preconditioning appeared soon afterward (Axelsson, 1972; Concus et al, 1976). See Golub and O’Leary (1989) for an early history, and Freund et al. (1992) and Saad and Van der Vorst (2000) for more complete recent overviews of Krylov subspace methods.

When paired with a suitable preconditioner, CG is a robust and efficient method for solving (4) when A is symmetric positive definite. In multiphysics problems, however, few circumstances are encountered where A is symmetric positive definite. Numerous generalizations of CG to non-symmetric problems that seek to exploit the optimality properties of CG have been proposed. Two principal approaches have been taken. One is to maintain r k ⊥ K k ( r 0 , A ) through an explicit Gram–Schmidt orthogonalization process. This approach constructs an orthonormal basis of K k ( r 0 , A ) and leads to the Arnoldi family of iterative methods, of which the generalized minimum residual (GMRES) method (Saad and Schultz, 1986) is the best-known variant. The main disadvantages of this approach are the storage requirements for the basis of K k ( r 0 , A ) and growth in the operation count for incrementally constructing that basis. Practical approaches to dealing with this difficulty amount to setting a maximum dimension of the Krylov subspace. Alternatively, a three-term recurrence relationship can be used. This approach is unable to maintain orthogonality, but it is inexpensive and can often be effective when paired with a suitable preconditioner. This strategy leads to the Lanczos-based family of methods. BiCGStab (Van der Vorst, 1992) is perhaps the most widely used method in this class. For general non-symmetric A , maintaining orthogonality while using short recurrence relationships is not possible (Faber and Manteuffe, 1984). Flexible Krylov variants (see e.g., Saad, 1993; Golub and Ye, 1999; Szyld and Voge, 2001; Vogel, 2007; Chen et al., 2012), where the preconditioner may change across iterations, offer promise for multiphysics scenarios. High-performance implementations of Krylov subspace methods are widely available (Balay et al., 2012; Falgout et al., 2011; Heroux et al., 2005).

Preconditioning

The key observations for Krylov subspace methods are that { K k ( r 0 , A ) } k = 1 N d is a nested sequence of subspaces that approximates the range of A and that the efficiency of a Krylov subspace method is determined by how well (4) is represented by its restriction to a low-dimensional Krylov subspace. Preconditioning is a technique for transforming (4) to an equivalent system where a lower-dimensional Krylov subspace provides a better approximation. This amounts to the introduction of a preconditioner P ≈ A − 1 . Typical formulations include left, right, and split preconditioning:

P A u = P f , ( l e f t ) A P P − 1 u = f , ( r i g h t ) P 1 A P 2 P 2 − 1 u = P 1 f , ( s p l i t )

Observe, in particular, that for left preconditioning, the sequence of approximations { u k } now resides in the Krylov subspaces { K k ( P r 0 , P A ) } k = 1 N d . To illustrate with an extreme case, if P = A − 1 , then P r 0 = u ⋆ − u 0 (i.e. the initial error), and the Krylov subspace method converges in one iteration. Thus, the art in selecting a preconditioner is managing the trade-off between the cost of applying the preconditioner and the reduction in the number of required iterations.

Many strategies are available for constructing a preconditioner. Early efforts were based on the notion of incomplete factorizations (Meijerink and Van der Vorst, 1977), and much research into black-box, algebraic preconditioning techniques is built around this idea. Using a stationary iterative method as a preconditioner was first explored by Axelsson (1972), and preconditioning with a multigrid method was first explored by Kettler (1982).

Parallel computing introduces new requirements for preconditioners that have a high ratio of computation to communication. A popular preconditioning approach is to use sparse direct solvers or incomplete factorizations in serial on each parallel process, in which each process solves only the portion of the overall problem that resides on its subdomain. Since these methods follow from the natural parallel decomposition of a problem into tasks, they are referred to as domain decomposition methods (Smith et al., 1996; Toselli and Widlund, 2005). Such solvers make good reuse of existing serial solver infrastructure, can attain excellent flop performance, work well for problems posed on unstructured meshes, and naturally couple multiphysics processes tightly within each solve.

Consider the full computational domain Ω as being given by a set of p non-overlapping subdomains Ω 1 , … , Ω p ; that is, Ω = ∪ i = 1 p Ω i . Then we can partition the unknowns of (4) into two subsets: the first consists of unknowns belonging to a single subdomain, u I , and the second consists of unknowns that lie on inter-subdomain boundaries and couple subdomains together, u Γ . Applying this partition to the matrix A as well results in a block linear system of the form

A I , I A I , Γ A Γ , I A Γ , Γ u I u Γ = f I f Γ

5

where the matrix A I , I is block diagonal. We further define restriction operators R i that take a full set of unknowns, u , and return only the set of unknowns on the interior of Ω i for each i = 1 , … , p .

The simplest domain decomposition method is block Jacobi, where all the interprocessor couplings are ignored (that is, Γ is empty) and (5) is solved by using only the decoupled subdomain solves

u k + 1 = u k + ∑ i = 1 p R i T A i − 1 R i f − A u k

6

where A i is a block of A I , I belonging to Ω i . A block Gauss–Seidel strategy, which sacrifices a lower degree of parallelism for faster convergence, can also be formulated, provided some coloring of the subdomains { Ω i } i = 1 p is available. For problems with indefinite or singular A i , artificial Dirichlet, Neumann, or Robin boundary conditions are typically employed at interprocessor interfaces. Because of the resulting decoupling between each subdomain, these methods do not perform well at large scales, especially for problems whose Green’s functions have global support, for example elliptic and parabolic PDE problems. As a result, various methods have been proposed to increase the coupling between subdomains and attain improved scaling performance.

Additive Schwarz methods (Dryja and Widlund, 1987; Cai, 1989; Chan and Zou, 1994; Smith et al., 1996) increase coupling by extending each subdomain to overlap with neighboring subdomains, increasing the size of each A i matrix and redefining each R i appropriately, then using (6). Extensions of this method include modification of the prolongation matrix R i T (Cai and Sarkis, 1999) and inclusion of coarse-grid solves (Bramble et al., 1990). A multiplicative variant, using the same components of the block Jacobi strategy above, is also possible.

Other classes of domain decomposition methods (e.g., FETI (Farhat and Rouz, 1991), BDD (Mandel, 1993) and BDDC (Dohrmann, 2003)) instead increase subdomain coupling by using a non-empty interface set Γ and following a Schur complement-like formulation to solve the full problem

u k + 1 = u k + I − A I , I − 1 A I , Γ 0 I A I , I − 1 0 0 S − 1 I 0 − A Γ , I A I , I − 1 I f − A u k

where S = A Γ , Γ − A Γ , I A I , I − 1 A I , Γ is the Schur complement matrix. Primary differences in methods following this approach lie in their choice of artificial boundary conditions at interior subdomain interfaces and how much of Γ they include, typically choosing only a small subset of the interface problem within S − 1 at each iteration. Extensive support for constructing domain-decomposed preconditioners is available in high-performance libraries (Balay et al., 2012; Heroux et al., 2005).

3.1.2 Methods for systems of nonlinear equations

Perhaps the simplest and most robust approach for solving the systems of nonlinear equations (3) is the fixed-point iteration, also known as the Picard or nonlinear Richardson iteration. These methods are performed through reformulation of the root-finding problem (3) into a fixed-point problem u = G ( u ) , through definition of a fixed-point iteration function, for example

G ( u ) := u − α F ( u )

7

where α > 0 is a fixed-point damping parameter that is typically chosen to be less than one. Fixed-point methods then proceed through the iteration

u k + 1 = G ( u k )

8

with the goal that ∥ u k + 1 − u k ∥< ε . A particular attraction of the fixed-point iteration is that if the iteration function is a contraction, that is, if there exists some γ ∈ ( 0 , 1 ) such that

∥ G ( u ) − G ( v ) ∥≤ γ ∥ u − v ∥

for all vectors u and v in a closed set containing the fixed-point solution u ⋆ , the fixed-point iteration is guaranteed to converge based on the Banach fixed-point theorem (also called the contraction mapping theorem (Banach, 1922; Granas and Dugundji, 2003)). This convergence is typically linear, however, and can result in slow convergence even when starting with a good initial guess at the solution.

Fixed-point iterations can converge slowly, in particular when the contraction constant γ is close to one. Newton’s method (Algorithm 3) offers faster (up to quadratic) convergence. Direct computation of δ u in Algorithm 3 may be prohibitively expensive for large-scale problems. Inexact Newton methods (Dembo et al., 1982) generalize Algorithm 3 by allowing computation of δ u with an iterative method, requiring only that ∥ J ( u k − 1 ) δ u + F ( u k − 1 ) ∥< ε k for some carefully selected set of tolerances { ε k } ; see Eisenstat and Walker (1994, 1996) for theoretical details. Newton–Krylov methods are important variants of inexact Newton methods in which δ u is computed with a Krylov subspace method. This choice is advantageous because the only information required by the Krylov subspace method is a method for computing Jacobian–vector products J ( u k ) v . A key consequence of this reliance on only matrix–vector products is that if these directional derivatives may be approximated, the Jacobian matrix J ( u k ) is never itself needed. The JFNK method exploits this approach, using a finite-difference approximation to these products

J ( u k ) v ≈ F ( u k + σ v ) − F ( u k ) σ

9

where σ is a carefully chosen differencing parameter (Brown, 1987; Kelley, 1995) and F is sufficiently smooth. This important specialization facilitates use of an inexact Newton method by eliminating the need to identify and implement the Jacobian. Of course, efficiency of JFNK depends critically on preconditioning the inner Krylov subspace method, and the need to track changes in the Jacobian places a premium on preconditioners with low setup cost; see Anderson et al. (1999), Gropp et al. (2000), Pernice and Tocci (2001), Chacón and Knoll (2003), Knoll and Keyes (2004), Mousseau (2004), Knoll et al. (2005), Reynolds et al. (2006) and Yeckel et al. (2009). Implementations of inexact Newton methods are available in several high-performance software libraries (Balay et al., 2012; Heroux et al, 2005; Hindmarsh et al., 2005).

Fixed-point and inexact Newton methods can be used to solve multiphysics problems in a fully coupled manner, but they can also be used to implement multiphysics coupling strategies such as Algorithms 1 and 2. With an implicit method for one of the components, one can directly eliminate it in a linear or nonlinear Schur complement formulation. The direct elimination process for u 1 in the first equation of the equilibrium system (1a), given u 2 , can be symbolically denoted

u 1 = G ( u 2 )

with which the second equation (1b) is well defined as

F 2 ( G ( u 2 ) , u 2 ) = 0

and each iteration thereof requires subiterations to solve (in principle to a high tolerance) the first equation. Unless the first system is much smaller or easier than the second (a special case of merit), this is not likely to be an efficient algorithm, but it may have robustness advantages.

If the problem is linear

F 1 ( u 1 , u 2 ) = f 1 − A 11 u 1 − A 12 u 2 F 2 ( u 1 , u 2 ) = f 2 − A 21 u 1 − A 22 u 2

10

then F 2 ( G ( u 2 ) , u 2 ) = 0 involves the traditional Schur complement

S = A 22 − A 21 A 11 − 1 A 12

If the problem is nonlinear and if Newton’s method is used in the outer iteration, the Jacobian

d F 2 d u 2 = ∂ F 2 ∂ u 1 ∂ G ∂ u 2 + ∂ F 2 ∂ u 2

is, to within a sign, the same Schur complement.

Similar procedures can be defined for the evolution problem, which, when each phase is implicit, becomes a modified rootfinding problem on each timestep with the Jacobian augmented with an identity or mass matrix.

If the types of nonlinearities in the two components are different, a better method may be nonlinear Schwarz or ASPIN (Cai and Keyes, 2002). In nonlinear Schwarz, one solves component subproblems (by Newton or any other means) for componentwise corrections

F 1 ( u 1 k − 1 + δ u 1 , u 2 k − 1 ) = 0 F 2 ( u 1 k − 1 , u 2 k − 1 + δ u 2 ) = 0

and uses these legacy componentwise procedures implicitly to define modified residual functions of the two fields

G 1 ( u 1 , u 2 ) ≡ δ u 1 G 2 ( u 1 , u 2 ) ≡ δ u 2

One then solves the modified rootfinding problem

G 1 ( u 1 , u 2 ) = 0 G 2 ( u 1 , u 2 ) = 0

by any means. If one uses Algorithm 3 to solve the modified problem, the Jacobian is

∂ G 1 ∂ u 1 ∂ G 1 ∂ u 2 ∂ G 2 ∂ u 1 ∂ G 2 ∂ u 2 ≈ I ∂ F 1 ∂ u 1 − 1 ∂ F 1 ∂ u 2 ∂ F 2 ∂ u 2 − 1 ∂ F 2 ∂ u 1 I

11

which clearly shows the impact of the non-dimensionalized cross-coupling in the off-diagonals, since the diagonal blocks constitute the identity. (This makes sense since the componentwise rootfinding problems are already solved at this outer level.) In practice, the outer Newton method must converge in a few steps for ASPIN to be worthwhile, since the inner iterations can be expensive. In nonlinear Schwarz, the partitioning of the global unknowns into u 1 and u 2 need not be along purely physical lines, as they would be if they came from a pair of legacy codes. The global variables can be partitioned into overlapping subsets, and the decomposition can be applied recursively to obtain a large number of small Newton problems. The key to the convenience of nonlinear Schwarz is that the outer iteration is Jacobian-free and the inner iterations may involve legacy solvers. (This is also true for JFNK, where the inner preconditioner may involve legacy solvers.)

Unless linear preconditioning of the global system is required to complement what is in effect nonlinear preconditioning through the inner solves, only distributed vector operations are needed to migrate legacy codes requiring coupling to this framework.

Error bounds on linear and nonlinear block Gauss–Seidel solutions of coupled multiphysics problems are derived in Whiteley et al. (2011), wherein it is shown how the norm of the product of the off-diagonal blocks in (11) enters the analysis. Consider the linear case (10) for which the Jacobian of the physics-blocked preconditioned system is

I A 11 − 1 A 12 A 22 − 1 A 21 I

Let the product of the coupling blocks be defined by the square matrix C ≡ A 11 − 1 A 12 A 22 − 1 A 21 , the norm of which may be bounded as | | C | | ≤ | | A 11 − 1 | | | | A 12 | | | | A 22 − 1 | | | | A 21 | | . Whiteley et al. (2011) show that, provided | | C | | ≤ 1 , any linear functional of the solution ( u 1 , u 2 ) of (10) solved by physics-blocked Gauss–Seidel (Algorithm 1) satisfies conveniently computable bounds in terms of residuals of the individual physics blocks of (10). The cost of evaluating the bounds involves the action of the inverse uniphysics operators A 11 − 1 and A 22 − 1 on vectors coming from the uniphysics residuals and the dual vectors defining the linear functionals of interest. The extra work required to evaluate the bounds provides confidence that block Gauss–Seidel has been iterated enough to produce sufficiently accurate outputs of interest, be they point values, averages, or fluxes, etc. The required actions of the inverse uniphysics operators are a by-product of an implicit method for each phase. The nonlinear case is similar, except that the Jacobian matrices from each physics phase that comprise C may change on each block Gauss–Seidel iteration.

3.2.1 Methods for coupling multiphysics components in space

Multiphysics applications require some means to exchange information between physics components. Preserving important properties, such as conservation or order of accuracy, is a principal concern. In the simplest case, all the physics components are defined on the same physical domain Ω and use the same mesh Ω h . In this situation, the coupling is defined by the physical model, and no error is incurred in transferring data between meshes. Numerically, some additional interpolations may be needed to align variables with different centerings or to achieve some consistency between different discretization schemes. In addition, the resolution must be fine enough to meet the requirements of all the physics components that share the mesh. More often, each physics component uses a mesh that is tailored to best represent the physics. These meshes can either fully or partially overlap and can also introduce different discrete representations of the geometry of Ω . (A similar need is encountered for arbitrary Lagrangian-Eulerian (ALE) methods, where periodic mesh remapping is needed.) Moreover, different physics components can be defined on distinct physical domains Ω i (where i indexes the physics components) and are coupled across lower-dimensional interfaces Γ i j = ∂ Ω i ∩ ∂ Ω j , i ≠ j , between the domains, which may also evolve in time. In this case, the meshes Ω i h i and Ω j h j may or may not match along the interface. When the restrictions Γ i j h i and Γ i j h j of the meshes to the common interface Γ i j do not match, information transfer across the interface must account for these differences in the discrete representation of Γ i j . In some cases the picture can be further complicated when Γ i j is determined as part of the calculation. The rest of this section provides a brief overview of methods currently used by multiphysics applications, and Section 4.2.3 considers related software issues.

At its core, transferring discrete fields between different meshes is a form of interpolation that is determined by how the field is represented on the source and target meshes. Search algorithms that determine the support of the interpolation on the source mesh are critical for constructing efficient transfer operations. Tree structures such as octrees, alternating digital trees, and k-d trees are used when the points to be evaluated are randomly distributed. Mesh topological adjacencies are used when the points are ordered such that the next point is from a nearby element. When the geometries of the discrete domains do not match, a destination point may lie outside the discrete source domain. In this case extrapolation can be used; ideally, information about the physical domain and boundary conditions will inform the extrapolation process.

In more complex cases, the transfer process must satisfy various kinds of local and global constraints. Examples include matching the integral of the solution and guaranteeing a divergence-free property (Bochev and Shashkov, 2005). In many circumstances, local conservation constraints are imposed. For a mesh transfer T to be conservative in a region R ⊂ Ω , it must satisfy

∫ R T f = ∫ R f

where f is the source grid data. Enforcement in a discrete sense can be represented as

∫ R ϕ i T f = ∫ R ϕ i f

where { ϕ i } is a set of test functions. For interpolation, they are simply delta functions at points on the destination mesh. If f is represented in a Galerkin formulation, { ϕ i } is the basis that is used on the destination mesh. Because f and ϕ i are only piecewise smooth on source or target meshes, the product ϕ i f is in general smooth only within the intersection of source and target elements. Therefore, using quadrature rules over source or target elements to compute ∫ ϕ i f d x may violate the regularity assumption of quadrature rules, and large errors and even instabilities can result from computing these integrals using quadrature rules based only on the source or destination mesh (Jiao and Heath, 2004a; Jaiman et al., 2011b). Knowledge of the complete intersection of the source and target meshes can then be used together with appropriate quadrature rules to construct T (Grandy, 1999; Jiao and Heath, 2004a,b,c; Solin et al., 2010; Jaiman et al., 2011a).

Transferring information across a shared interface between two domains depends on the smoothness of solutions on either side of the interface, discretization schemes, and alignment of the two meshes on either side of the interface as well. When the meshes align at the interface, direct transfer can be used effectively. For an illustration, consider coupling Navier–Stokes flow on the surface and Darcy flow in subsurface porous media (Figure 12) (Chidyagwai and Rivière, 2010). The free-flow region is denoted Ω 1 , and the porous medium is represented as several layers (regions A, B, C) with different permeability and faults (region D). The interface conditions enforce continuity of the normal component of the velocity across the interface, a force balance at the interface, and the Beavers–Joseph–Saffman model relating the tangential free-flow velocity to the shear stress at the interface (Girault and Riviére, 2009). Let a N S represent the bilinear form for the Navier–Stokes flow discretized with a continuous finite-element method, and a D represent the bilinear form for the Darcy equations discretized with a discontinuous Galerkin method to accurately capture irregularities in the subsurface region. A grid is used that covers the entire domain but captures the interface between the two models. This means that the bilinear form a Γ for the interface coupling conditions across the interface can easily be constructed by simply using values of the solution obtained at the interface. The fully coupled problem is then given by

a N S ( u , v ) + a D ( u , v ) + a Γ ( u , v ) = L ( v )

12

where v is drawn from suitable test spaces and L is a forcing term.

A fully coupled solution method is the most stable and accurate approach to solving (12). However, this leads to a large system of nonlinear equations that may be difficult to precondition and expensive to solve. One way to lower the cost of solving (12) is to use a two-grid technique (Chidyagwai and Rivière, 2011). In the first step, (12) is solved on a coarse mesh with resolution H ; in the second step, decoupled “single physics” problems are solved on a finer mesh, with resolution h = H 2 . Boundary data along the interface for the decoupled problems is obtained from the coupled coarse mesh solution. This approach is nearly an order of magnitude faster and produces solutions with accuracy comparable to the fully coupled approach, but it is less stable and accurate. The solution obtained with this algorithm (Figure 12) shows smooth velocity contours, even into the fault regions. For details, see Chidyagwai and Rivière (2011).

The mortar finite-element method is another approach for treating interface conditions (Bernardi et al., 1993, 1994). It is especially effective when the meshes on the two sides of the interface do not match and when some quantity (such as flux or stress) is needed at the interface. The mortar finite-element method yields a variationally consistent formulation and provides an abstract framework for analyzing the weakly coupled discrete model. The algebraic problem can then be stated as a saddle point or, equivalently, as a constrained system where static condensation of the Lagrange multiplier is realized and the coupling conditions are imposed on the spaces. Mortar formulations are widely used in the context of nonlinear contact problems in structural mechanics and replace simple node-to-node or node-to-segment penalty approaches; see, for example, Wohlmuth (2011). The surface traction plays the role of the Lagrange multiplier, which has to satisfy a sign condition in the normal direction, which is due to the non-penetration condition, and the friction law, which combines normal and tangential surface forces. The inequality conditions can be reformulated as a weakly consistent saddle point system in terms of a nonlinear complementarity function. Then semi-smooth Newton methods provide a superlinear solver, and an efficient algorithmic implementation can be based on the choice of biorthogonal basis functions, resulting in a diagonal mass matrix.

Figure 13 illustrates the acoustic-structure interaction for a three-dimensional loudspeaker geometry; see Flemisch et al. (2010). The problem has been formulated in a weakly consistent mortar setting, and highly non-matching meshes are applied for the numerical simulation.

OPEN IN VIEWER

Figure 13.

Schematic view of the geometric and acoustic pressure field in acoustic-structure interaction for a three-dimensional loudspeaker geometry.

Of crucial importance for the stability of the method is a uniform inf–sup condition for the pairing between dual and primal variables. In order to obtain optimal a priori error results, the discrete Lagrange multiplier space must additionally satisfy a best-approximation property. With these two ingredients, one can easily construct locally defined higher-order methods. Moreover, the influence of curvilinear interfaces can be analyzed within the framework of variational crimes, and the weak coupling concept can be easily combined with blending or isogeometric finite elements.

Although mortar techniques were originally designed for non-overlapping domain decomposition methods, the idea of weakly coupled subdomain formulations can be used to design two-scale overlapping subdomain approaches. Globally, a simple, coarse-scale discrete model is used, whereas a more sophisticated model is applied locally on a small patch. Different equations, mesh sizes, and geometries can be considered. Here, one has to distinguish between one- and bidirectional couplings. One-directional coupling can be regarded as a postprocess and allows locally for a higher resolution of the numerical solution. The possibly nonlinear system has the structure of a lower triangular block system and thus can easily be solved in terms of a Gauss–Seidel-type solver, possibly in combination with massive parallel algebraic multigrid or domain decomposition techniques. For bidirectional coupling, iteration between the small-patch and coarse-scale models is required to obtain a physically correct solution.

3.2.2 Methods for coupling multiphysics components in time

Coupling time-dependent continuum models depends on the temporal resolution requirements of the physics components. Typically, each component of the solution is allowed to evolve governed by its own physics. However, there can be wide disparity in the timescales of the constituent physics components, which may in turn be quite different from the timescale of the problem being solved. Replacing instantaneous interactions between different physics with discrete interactions can introduce subtle stability and accuracy issues. Likewise, the multiphysics problem may conserve some global quantity (such as energy), but this property may not be preserved by the simulation.

In evolutionary problems (2a)–(2b), the spatial discretization leads to a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs) in general; that is, ODEs possibly with algebraic constraints. A timestepping approach is then used to advance the solution in time. This strategy is called the method of lines. For exposition brevity, we initially focus on ODEs, with initial conditions u ( t 0 ) and computing values of the solution { u ( t n ) } n = 1 N t at discrete times t 1 < … < t N t .

In an explicit method, future states are given by a formula involving known computed quantities that can be evaluated directly. The simplest example is forward Euler:

u ( t n + 1 ) = u ( t n ) + Δ t f ( u ( t n ) )

where Δ t = t n + 1 − t n . Explicit methods are inexpensive per iteration because they require a fixed number of function evaluations to complete the step; however, they typically have severe stability restrictions and, therefore, are suited only for non-stiff problems.

A fully implicit method advances the solution in time by using current information and an inversion process with no explicit steps. The simplest example is backward Euler:

u ( t n + 1 ) = u ( t n ) + Δ t f ( u ( t n + 1 ) )

Fully implicit methods are relatively expensive because of the need to solve (non)linear systems at each step; however, they have favorable stability properties (Shampine and Gear, 1979; Higham and Trefethen, 1993). Inclusion of explicit elements gives rise to semi-implicit methods (Zhong, 1996; Giraldo, 2005).

Spatial discretization of multiphysics problems (2a)–(2b) may result in problems in which the tendency f ( u ) has components with widely different dynamics, often classified as stiff and non-stiff (Ascher et al., 1997; Verwer and Sommeijer, 2004; Shampine et al., 2006; Hundsdorfer and Ruuth, 2007; Constantinescu and Sandu, 2010b; Cionnors and Miloua, 2011). This case, called additive splitting, generally corresponds to problems in which the physics components act in the same spatial domain, such as reactive transport or RHD (Ruuth, 1995; Giraldo et al., 2003; Estep et al., 2008a; Giraldo and Restelli, 2009; Giraldo et al., 2010; Durran and Blossey, 2012). When a partitioned time-integration method is used with physics components defined on the same spatial region, some work must be done to synchronize the timestepping in each physics component. When one component has a much faster timescale relative to another, a simple strategy is to use an integer multiple of the faster timescale for the slow timescale. In component partitioning some solution components correspond to a stiff subproblem, and the others, to a non-stiff one; examples are FSI and flow with different speeds (Constantinescu and Sandu, 2007; Constantinescu et al., 2008; Sandu and Constantinescu, 2009). Component partitioning can be used, for instance, if physics components are defined on different spatial domains; the same subcycling idea can work, but some protocol for communicating information across the interface is needed to ensure accuracy and stability. This leads to the following two conceptual types of partitioning:

f ( u ) = f 1 ( u ) + f 2 ( u ) a d d i t i v e p a r t i t i o n i n g

13

f ( u ) = f 1 ( u 1 , u 2 ) f 2 ( u 1 , u 2 ) c o m p o n e n t p a r t i t i o n i n g

14

Simple transformations can be used to cast one type of splitting into the other (Ascher et al., 1997; Constantinescu and Sandu, 2010b). We will henceforth use form (13) in which f 1 is relatively non-stiff and its presence points one to use an explicit integrator for efficiency, whereas f 2 is relatively stiff and requires an implicit integrator for stability reasons. For instance, in (13), f 1 can incorporate all slow transport components and f 2 can incorporate all stiff kinetic terms, and with (14), u 1 may represent coarse-grid points and u 2 , finely refined states or structure and fluid states, respectively. Partitioned integrators were introduced to solve such problems efficiently.

Partitioned methods use time integrators with different properties or timesteps for each (physics) component. They allow black-box implementation of the physics components and different timesteps in each physics component. Nevertheless, this formalism allows one to rigorously analyze the stability and accuracy of the compound scheme. Partitioned methods can be tightly coupled (in time); one simple example is the forward–backward Euler method:

u ( t n + 1 ) = u ( t n ) + Δ t f 1 ( u ( t n ) ) + Δ t f 2 ( u ( t n + 1 ) )

Strang splitting methods (Strang, 1968) are also particular cases of (loosely coupled) partitioned schemes:

u ⋆ = u ( t n ) + Δ t 2 f 2 ( u ⋆ ) u ⋆ ⋆ = u ⋆ + Δ t f 1 ( u ⋆ ) u ( t n + 1 ) = u ⋆ ⋆ + Δ t 2 f 2 ( u ⋆ ⋆ )

The caveats of such a strategy include order reduction and low accuracy as illustrated in Knoll et al. (2003). Also, reduced order or even loss of consistency may occur at boundaries for additive splitting methods. One solution to restore consistency in this case is to apply intermediate boundary conditions that are specific to each substep in the method (LeVeque and Oliger, 1983; Connors et al., 2012). The temporal splitting error can be fairly easily avoided in tightly coupled schemes (Kennedy and Carpenter, 2003; Constantinescu and Sandu, 2010b). Other tightly coupled strategies involve embedding explicit evaluations inside an implicit solver (Kadioglu et al., 2010b).

We differentiate two cases that have a more practical flavor: multirate and IMEX. Multirate methods are partitioned methods that use the same scheme in each partition but with different timesteps (Gear and Wells, 1984; Bartel and Günther, 2002; Logg, 2003a,b, 2004, 2006; Constantinescu and Sandu, 2007; Savcenco et al., 2007; Constantinescu et al., 2008; Sandu and Constantinescu, 2009); this goes back to subcycling. Multirate schemes are also used to address accuracy considerations. One instance is applications that employ temporal and spatial adaptive grids. Another instance is shown in Figure 14, which illustrates a multirate scheme that uses smaller timesteps for the reaction component in a reaction–diffusion problem.

OPEN IN VIEWER

Figure 14.

Illustration of multirate operator-splitting time discretization.

An IMEX method is a partitioned method that uses an implicit integrator on one partition and an explicit integrator on the other. Therefore, an optimal trade-off between efficiency and stability becomes possible (Crouzeix, 1980; Kennedy and Carpenter, 2003; Constantinescu and Sandu, 2010b; Filbet and Jin, 2010). Moreover, it is generally desirable to use an explicit integrator when nonlinear or non-differentiable spatial discretizations are used, such as shock resolving schemes (LeVeque, 2002). In this case IMEX alleviates the difficulties posed by such a problem on the nonlinear solver, which is applied only on the stiff component. Furthermore, IMEX methods generalize fully implicit [ f 1 ≡ 0 , f → f 2 ], and fully explicit [ f 2 ≡ 0 , f → f 1 ] methods, and can provide the best balance between accuracy and efficiency. In both cases the temporal splitting error can be made of the same order as the schemes used separately; however, stiffness leakage can be a concern (Kennedy and Carpenter, 2003). Partitioned methods based on the classical linear multistep (Ascher et al., 1995; Ascher and Petzold, 1998; Hundsdorfer and Ruuth, 2007), Runge-Kutta (Ascher et al., 1997; Kennedy and Carpenter, 2003; Verwer and Sommeijer, 2004; Constantinescu and Sandu, 2007) and extrapolation methods (Gragg and Stetter, 1964; Constantinescu and Sandu, 2010b) have been proposed and successfully applied in practice (see e.g., Butcher, 2008). In Figure 15 we illustrate the temporal convergence rate for a two-dimensional Euler rising-bubble problem by using high accuracy in space (Giraldo and Restelli, 2009) and different IMEX time-integration schemes. The splitting is based on linearizing the PDE: extracting the stiff wave components, L , and then integrating (2a)–(2b) as

∂ t u = ( f ( u ) − L u ) + L u

where the first term is treated explicitly and the last linear term, implicitly. The efficiency advantage of IMEX schemes increases with temporal discretization order. Given a timestep selected for accuracy, in terms where the timestep is below the stability limit, explicit integration is employed for efficiency, and in terms where the stability is limiting, implicit integration is employed.

OPEN IN VIEWER

Figure 15.

Global convergence and timing of implicit–explicit (IMEX) schemes for a two-dimensional Euler rising-bubble test problem. Method SI-BDF is second order, and the rest have the order suffixed to their name.

Interest in IMEX schemes is far from limited to multiphysics problems. An early use of IMEX temporal integration, with consideration of stability and numerical efficiency, was in solid mechanics and heat transfer applications (Hughes and Liu, 1978b,a). A common source of localized stiffness apart from multiphysics problems is stiffness induced by adaptive meshes, for example in boundary layers in fluid mechanics or electromagnetism. In Kanevsky et al. (2007), IMEX methods are introduced for this purpose, and implicitness is employed in the diagonal terms only of the highly refined portions of the domain, in pursuit of further efficiency. Truly multiphysics applications of IMEX methods now abound. An algorithmic variant that exploits the JFNK method for the implicit part and fully eliminates splitting error is employed by Kadioglu et al. (2010a,b) in RHD problems. In these studies, a second-order IMEX technique resolves the explicitly treated non-stiff hydrodynamics subsystem inside every nonlinearly implicit Newton step of the radiation subsystem. This causes each Newton correction of the implicit block to be felt immediately by the unknowns governed by the explicit block; however, there is no additional software or data structure complexity for the implicit solver.

Related to the self-consistent IMEX method enabled by the use of JFNK for the outer implicit part is the concept of “implicit balancing” introduced in Knoll et al. (2003), in which in multirate, multicomponent problems, all nonlinear couplings between the components are evaluated at the same point in time, rather than just semi-implicitly. Semi-implicitness is sufficient to preserve stability, but its accuracy properties are limiting. As illustrated for radiation diffusion, MHD, and shallow-water models, implicit balancing permits second-order schemes to attain true second-order temporal accuracy, as opposed to first-order accuracy with semi-implicitness only.

Stability considerations for partitioned methods

Partitioned methods use multiple time integrators with well-known stability properties when applied to the compound problem (2a)–(2b). When applied on partitions, however, the stability properties of individual methods do not transfer directly to the stability of the compound problem. Consider two integrators that are stable when applied separately on each partition but where the resulting compound method may lead to unstable overall computations (Frank et al., 1997; Ropp et al., 2004; Ropp and Shadid, 2009; Constantinescu and Sandu, 2010b). Such a situation is illustrated in Figure 16 for extrapolated IMEX schemes, where the dependence between the implicit and explicit stability regions is clearly visible. Note that larger stability regions for the implicit method correspond to smaller stability regions for the extrapolated explicit methods. This behavior must be taken into account when selecting an extrapolated IMEX scheme for multiphysics problems. If the implicit term is used to integrate symmetric discretizations of hyperbolic operators that lead to pure imaginary eigenvalues, an A-stable (see e.g., Lambert, 1991) implicit method is needed, and one may expect a penalty on the explicit step size because the explicit stability region is smaller (Method A). Discretizations of parabolic operators that lead to negative real eigenvalues can use an A ( α ) -stable implicit method for α < 90 , so the larger stability region for the explicit operator allows larger stable explicit steps (Method B). These examples tacitly assume that the implicit and explicit operators are simultaneously diagonalizable, which happens almost exclusively only for scalar problems; however, in practice, this analysis typically holds for vector-valued problems, especially when stiff components are present (Frank et al., 1997).

OPEN IN VIEWER

Figure 16.

Stability regions of implicit–explicit (IMEX) pairs for extrapolated IMEX schemes. The implicit methods are A ( α ) -stable for different values of α . The explicit stability regions indicate different orders of extrapolation. The stability region of the IMEX scheme contains the cross-product of the regions for these pairs. This figure is adapted from Constantinescu and Sandu (2010a); 2010 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved.

Stability properties of partitioned timestepping methods are also affected by the treatment of interfaces. This can be illustrated by a simple model of heat conduction through two adjacent materials in subdomains Ω 1 and Ω 2 that are coupled across their shared and rigid interface Γ through a jump condition:

∂ u i ∂ t − ν i Δ u i = f i , x ∈ Ω i − ν i ∇ u i ⋅ n i = κ ( u i − u j ) , x ∈ Γ u i ( x , 0 ) = u i 0 ( x ) , x ∈ Ω i u i = g i , x ∈ ∂ Ω i ∖ Γ

15

for i , j = 1 , 2 , i ≠ j . Three time-integration strategies can readily be described. The fully coupled approach solves equation (15) with u 1 and u 2 and all terms evaluated at the same time level. This method is unconditionally stable, conserves energy, and converges optimally under standard smoothness assumptions. A second scheme is to again keep u 1 and u 2 at the same time level but allow the transmission conditions on Γ to lag in time. This IMEX scheme decouples the subdomain solves but is stable only when κ is small compared with ν 1 and ν 2 , together with a CFL-like stability restriction on the timestep. Here the trade-off is between the need for smaller timesteps and the reduced cost of solving easier problems in each subdomain. A third partitioned scheme uses interface data at the previous timestep from the other subdomain. Since each subdomain requires data from the other subdomain, this data-passing approach is less flexible but is stable when κ is large compared with ν 1 and ν 2 . The data-passing approach also does not have a stability constraint on the timestep; see Connors et al. (2009) for details and analysis. Partitioned time-integration methods are often used in a sequential fashion, as in Algorithm 2, first advancing a single physics component in time, followed by updating information on an interface, advancing a second physics component in time, and finally updating interface information from the second component. The strategy is often referred to as staggered computation, and the techniques used to update data on the interface play a crucial role in the stability of this approach, irrespective of the time-integration methods used for the component physics. In particular, there can be a time lag in the interface data that leads to an imbalance in the interface transmission conditions, even for strategies that are formally second-order accurate in time. For example, in FSI, a small error in the structural displacements may result in large errors in the fluid pressure. The stability and accuracy can be improved by using the combined interface boundary condition (CIBC) (Jaiman et al., 2011b), which solves the transformed conditions for the interface quantities based on their space derivatives and time derivatives. The CIBC technique allows for the construction of a staggered coupling procedure with a similar order of accuracy and stability as tightly coupled computations for some classes of fluid–structure coupling problems.

Despite the issues associated with partitioned and mixed approaches to multiphysics time integration, these elements can be combined to produce stable, accurate, and efficient simulations. Such an approach is illustrated by hierarchical adaptive mesh refinement (AMR), which we discuss in the context of lean, premixed turbulent combustion (Day and Bell, 2000). This is a demanding multiscale simulation involving fluid flow, chemical reactions, and heat transfer to resolve spatial scales that span two orders of magnitude and timescales that span three orders of magnitude. A finite-volume formulation is used to enforce crucial conservation properties. The AMR strategy is tailored to the mathematical structure of the problem to exploit the relationship between the scales and physical processes in the problem. The fluid equations are advanced in time by using a fractional-step projection formulation, whose elements include explicit advection (via a Godunov scheme), implicit diffusion (via a Crank–Nicolson predictor–corrector scheme), Strang splitting for the chemical reactions, and two discrete projection steps to enforce a modified, divergence-free constraint.

AMR introduces two complexities. The first is the need to compute discrete differential operators where the resolution changes; this computation is done through a simple direct interpolation of data at a coarse resolution to a finer resolution. Local piecewise quadratic interpolation of coarse data is used to align the coarse data with fine data and to calculate fine-scale fluxes at the interface (Almgren et al., 1998). (Local piecewise linear interpolation can also be used (Ewing et al., 1991). Failure to adequately address this issue can lead to non-physical artifacts in the solution (Choi et al., 2004).) In addition, a re-fluxing procedure (Berger and Colella, 1989) is employed to enforce conservation and free-stream preservation at these changes in resolution. The second complexity is local subcycling in time: a global timestep on the entire locally refined grid can be inefficient and will lead to a wide range of CFL numbers across the grid, compromising the accuracy of the Godunov advection scheme. In order to handle this, timesteps are chosen in each refinement region to maintain a CFL number close to one. Practically speaking, this means that a grid with resolution h requires twice as many timesteps as the next-coarser region with resolution 2 h (assuming the refinement ratio is 2 ). Overall, this AMR scheme is second-order accurate in time and space; see Pember et al. (1998) and Almgren et al. (1998) for details.

We note that it is difficult to devise partitioned time-integration schemes that are better than second-order accurate in time. Spectral deferred correction (SDC) (Frank and Ueberhuber, 1997; Dutt et al., 2000) is one potential way to improve time accuracy. Starting with a provisional solution obtained with a simple, low-order method, SDC generates a sequence of corrections that can lead to higher-order accuracy in the solution. Specifically, one writes the solution to the problem in terms of a Picard integral equation and evaluates the integral using interpolated values of the provisional solution. Spectral accuracy of the integral can be obtained by using Gaussian quadrature nodes in the quadrature. In principle, the order of accuracy can be increased up to the order of accuracy of the quadrature used to evaluate the Picard integral. Variations of SDC have also been developed for semi-implicit time-integration methods (Minion, 2003, 2004).

3.4.1 Operator splitting for reaction–diffusion equations

A classic example of M 3 OD is operator splitting for a reaction–diffusion equation. Accuracy considerations dictate the use of relatively small steps to integrate a fast reaction component. On the other hand, stability considerations over moderate to long time intervals suggest the use of implicit, dissipative numerical methods for integrating diffusion problems. Such methods are expensive to use per step, but relatively large steps can be used on a purely dissipative problem. If the reaction and diffusion components are integrated together, then the small steps required for accurate resolution of the reaction lead to an expensive computation.

In multirate operator splitting, the reaction and diffusion components are integrated independently inside each time interval of a discretization of time and “synchronized” in some fashion only at the nodes of the interval. The reaction component is often integrated by using significantly smaller substeps (e.g., 10 − 5 times smaller is not uncommon) than those used to integrate the diffusion component, which can lead to tremendous computational savings.

We illustrate some possible effects for multirate operator splitting using the well-known Brusselator problem

∂ u 1 ∂ t − ν 1 Δ u 1 = 0.6 − 2 u 1 + u 1 2 u 2 ∂ u 2 ∂ t − ν 2 Δ u 2 = 2 u 1 − u 1 2 u 2 s u i t a b l e i n i t i a l a n d b o u n d a r y c o n d i t i o n s

with ν 1 = ν 2 = 0.025 , over the one-dimensional spatial domain [ 0 , 1 ] . We discretize using:

a second-order, piecewise linear finite-element method in space with 500 elements;

We show the results in Figure 20. On the left, we plot the error versus timestep for different final times. We see that on sufficiently long time intervals, there is a critical timestep above which convergence fails. On the right, we plot a typical long time unstable solution. Note that the oscillations are not on the scale of the spatial discretization.

OPEN IN VIEWER

Figure 20.

Left: On moderate to long time intervals, there is a critical timestep above which convergence fails. Right: A typical unstable solution. See Estep et al. (2008a) for details; 2008 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved.

This example illustrates well the fundamental issue arising in solutions using M 3 OD. When divorced from the stabilizing effect of diffusion, the solution of an unstable reaction component is free to grow too rapidly. In this case, the reaction component is not constrained by boundary conditions imposed on the solution of the diffusion component, and consequently its values near the boundaries grow far too rapidly. This coupling either needs to be compensated for by applying an L-stable (see e.g., Lambert, 1991; Ascher and Petzold, 1998) diffusion solver in the diffusion solve (Ropp et al., 2004; Ropp and Shadid, 2009), or the coupling needs to be implicit for large timesteps.

3.4.2 Operator splitting for advection–diffusion equations

When splitting an advection operator from a diffusion operator, a subtle inconsistency may occur at boundaries. The ensuing error will pollute both the advection and diffusion components of the solution. In some cases, the inconsistency may occur only if advection is subcycled relative to diffusion.

Consider the advection–diffusion problem

∂ u ∂ t + a ∂ u ∂ x = ν ∂ 2 u ∂ x 2 , x ∈ ( x L , x R ) × ( 0 , T ] u ( x L , t ) = g L ( t ) , t ∈ ( 0 , T ] u ( x R , t ) = g R ( t ) , t ∈ ( 0 , T ] s u i t a b l e i n i t i a l c o n d i t i o n s

where a > 0 and ν > 0 are constants. An example of an operator-split calculation is to resolve advection by using a step of Lax–Wendroff between times t n and t n + Δ t = t n + 1 , then using a centered scheme with a Crank–Nicholson timestep for diffusion. We illustrate this using finite difference methods at grid points x i = x L + i Δ x for i = 0 , 1 , … , M + 1 with x M + 1 = x R . At a grid point x i and time level t n , denote the advection degrees of freedom as v i n and the diffusion degrees of freedom as w i n , with the operator-split solution degrees of freedom being denoted by u i n . The discrete algorithm is defined over a timestep by setting v i n = u i n for i = 0 , … , M + 1 , then updating the advection degrees of freedom:

v i n + 1 = v i n − a Δ t 2 Δ x ( v i + 1 n − v i − 1 n ) + a 2 Δ t 2 2 Δ x 2 ( v i + 1 n − 2 v i n + v i − 1 n ) v 0 n + 1 = g L ( t n + 1 ) i = 1 , … , M v M + 1 n + 1 = g R ( t n + 1 )

Alternatively, the right boundary value could be extrapolated since a > 0 . The updated advection data is used to initialize the diffusion step by setting w i n = v i n + 1 for i = 1 , … , M , with boundary values w 0 n = g L ( t n ) and w M + 1 n = g R ( t n ) . The diffusion update is

w i n + 1 = w i n + ν Δ t Δ x 2 ( w i + 1 n + 1 / 2 − 2 w i n + 1 / 2 + w i − 1 n + 1 / 2 ) , i = 1 , … , M w 0 n + 1 = g L ( t n + 1 ) w M + 1 n + 1 = g R ( t n + 1 )

where w i n + 1 / 2 = ( w i n + 1 + w i n ) / 2 . The operator-split approximation is then u i n + 1 = w i n + 1 , for i = 0 , … , M + 1 .

The truncation error for this method is found in the usual way by inserting u ( x i , t j ) for u i j and j = n , n + 1 . Away from boundaries ( 1 < i < M ), the method is consistent, with order Δ t 2 + Δ x 2 . Near the boundaries, however, an inconsistency has been introduced because the boundary data used in the scheme does not generally evolve in a way that is compatible with the split operators. Indeed, one may verify directly that the local truncation error is on the order of Δ t Δ x 2 at x i in the case where i = 1 . Special intermediate boundary values may be derived to amend this issue. See Connors et al. (2012) for examples, a posteriori error analysis, and other details.

3.4.3 Iterative solution of parabolic problems coupled through a common boundary

Operator splitting is an extreme form of M 3 OD because it is a blunt factorization of the diffusion–reaction or advection–reaction differential operator whose effects cannot be repaired by, for example, iteration. We next discuss a less severe example of M 3 OD that nonetheless presents some difficulties and puzzling behavior.

We consider a system of parabolic problems, where each equation is posed on a different domain but the domains share a common boundary Γ . In contrast to equation (15), we impose continuity of state and normal flux across Γ :

∂ u 1 ∂ t − ∇ ⋅ ( ν 1 ∇ u 1 ) = f 1 , x ∈ Ω 1 ∂ u 2 ∂ t − ∇ ⋅ ( ν 2 ∇ u 2 ) = f 2 , x ∈ Ω 2 u 1 = u 2 , x ∈ Γ ν 1 ∇ u 1 ⋅ n 1 = ν 2 ∇ u 2 ⋅ n 2 , x ∈ Γ s u i t a b l e b o u n d a r y c o n d i t i o n s o n t h e o t h e r b o u n d a r i e s s u i t a b l e i n i t i a l c o n d i t i o n s

where n i is the outward unit normal to Ω i along Γ for i = 1 , 2 . This example arises in conjugate heat transfer, representing the exchange of heat between two abutting solids with different temperatures. It is a subsystem in the exchange of heat between a solid and an encompassing fluid. This system also arises in other contexts, for example the core–edge coupling in models of a fusion reactor.

Such a system may be posed as a large, fully implicit system of equations. However, it is common practice to use some form of iterative timestep method. The different components often use different spatial and temporal discretizations and may even be solved with different codes.

A common approach combines Algorithm 1, as an inner iteration, with Algorithm 2 and is specified in Algorithm 4. The solves that are performed on each subdomain may use subcycling in time; recall that we use u 1 ( t n ) and u 2 ( t n ) to denote the nodal values of the numerical solutions at N t time nodes t 0 = 0 < t 1 < t 2 < … < t N t . Of course, we must discretize in space to create full discretizations. Many possible variations of this approach exist.

OPEN IN VIEWER

Algorithm 4.

Given initial values u i ( t 0 ) , i = 1 , 2 :

for n = 1 , … , N t do

Set u i 0 ( t n ) = u i ( t n − 1 ) , i = 1 , 2

for k = 1 , 2 , … , (until convergence) do

Solve for u 1 k ( t n ) with u 1 k ( t n ) = u 2 k − 1 ( t n ) o n Γ

Solve for u 2 k ( t n ) with ν 2 ∇ u 2 k ( t n ) ⋅ n 2 = ν 1 ∇ u 1 k ( t n ) ⋅ n 1 o n Γ

end for

Set u 1 ( t n ) = u 1 k ( t n ) a n d u 2 ( t n ) = u 2 k ( t n )

end for

The number of specified iterations in each step is an important factor. This approach can be interpreted as a fixed-point iteration, and if the fixed-point iteration converges, then by iterating sufficiently, this discretization converges to the fully implicit discretization. On the other hand, this approach often is implemented with only one iteration. Doing so may have significant effects on accuracy. More worrisome is that if the fixed-point iteration actually fails to converge, using only one or a few iterations will fail to reveal this fact. The numerical results are questionable in this circumstance.

To illustrate some of the behavior that can arise in this approach, we consider simple linear heat equations with diffusion coefficients ν 1 = 1 and ν 2 = 1 / 2 on domains Ω 1 = [ 0 , 1 ] and Ω 2 = [ 1 , 2 ] , respectively, coupled by imposing continuity of state and normal flux at Γ 12 = { 1 } . The solution should converge to a steady state. We solve using Algorithm 4 with continuous second-order finite-element methods in space with Δ x = 0.01 and implicit Euler with Δ t = 10 − 3 for each component. Each component is solved sufficiently accurately for the component errors to be negligible.

We plot the total error of the approximation in Figure 21. We see the expected increase during the initial transient. However, the increase in error that occurs over longer time intervals is not expected, since the solution should be converging to a steady state. We can estimate the error that arises from incomplete iteration at each step (we use only one iteration) and the error in the numerical flux that is used in the solution of u 2 . Note that u 1 produces a numerical flux that is only first order. We see that the error from incomplete iteration is significant during the initial rapid transient but becomes unimportant over long time intervals. However, the error due to transferring a numerical flux increases as time passes, and becomes the dominant feature. In fact, the pollution from transferring the numerical flux causes a loss of order in spatial accuracy. Analyses of the errors, as well as examples in two dimensions with both matching and non-matching grids along Γ , are given in Estep et al. (2008b, 2009).

OPEN IN VIEWER

Figure 21.

Results of the Gauss–Seidel algorithm for coupled parabolic problems. The component discretization errors are negligible. Left: the total error. Second from left: the contribution from incomplete iteration on each step. Second from right: the contribution from the passing numerical fluxes. Right: approximation showing loss of an order of accuracy in space as a result of the error in the numerical flux.

3.4.4 Solution of systems of elliptic equations with independent discretizations

In the next example, drawn from Carey et al. (2009), we consider an even more nominally benign form of M 3 OD, which nonetheless produces interesting behavior. We consider a two-component, “one-way” coupled system.

− Δ u 1 = sin ( 4 π x ) sin ( π y ) x ∈ Ω − Δ u 2 = b ⋅ ∇ u 1 x ∈ Ω u 1 = u 2 = 0 x ∈ ∂ Ω b = 2 π 25 sin ( 4 π x ) sin ( π x )

Here Ω = [ 0 , 1 ] × [ 0 , 1 ] . The “lower triangular” form of this system means that we can solve it as a coupled system or we can solve the first equation and then use the solution to generate the parameters for the second problem, as in Algorithm 5. In either case, we use independent discretizations of the two problems.

OPEN IN VIEWER

Algorithm 5.

Construct discretizations for the two equations

Compute the numerical solution u 1 of the first equation

Project u 1 onto the mesh for the second equation

Compute the numerical solution u 2 of the second equation

Solutions obtained by using a standard piecewise-linear, continuous finite-element method on uniform, highly refined meshes and Algorithm 5 are shown in Figure 22.

OPEN IN VIEWER

Figure 22.

Solutions of the component problems of the elliptic system computed on uniform meshes. See Carey et al. (2009) for details; 2009 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved.

Using uniform meshes, we evaluate the standard a posteriori error estimate for the second component problem and the quantity of interest u 2 ( 0.25 , 0.25 ) , ignoring any effect arising from error in the solution of the first component. We estimate the solution error to be ≈ 0.0042 . However, the true error is ≈ 0.0048 , giving a discrepancy of ≈ 0.0006 ( ≈ 13 % ) in the estimate. This is significant for these kinds of accurate a posteriori error estimates since the solutions are smooth.

If we adapt the mesh for the solution of the second component based on the a posteriori error estimate of the error in that component while neglecting the effects of the decomposition (Figure 23), the discrepancy becomes alarmingly worse. For example, we can refine the mesh until the estimate of the error in the second component is ≈ 0.0001 . However, we find that the true error is ≈ 0.2244 .

OPEN IN VIEWER

Figure 23:.

Results obtained after refining the mesh for the second component so that the a posteriori error estimate of the error only in the second component is less than 0.001 . The mesh for the first component remains coarse; consequently, the error in the first component becomes relatively much larger. See Carey et al. (2009) for details; 2009 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved.

4.1.1 Current practices

Current practices for multiphysics coupling in production codes, such those introduced in Section 2, must strike a balance among performance, software reuse, numerical accuracy, and robustness, while also addressing ease-of-use for developers at all levels of the multiphysics software stack. Exactly which practices and approaches are favorable to these goals is a subject of controversy and continued research. However, some common choices exist.

The importance of software reuse has strongly favored operator-splitting coupling methodologies, possibly at the expense of numerical accuracy. Typical workhorse codes in key physics areas such as CFD, particle transport, atomistic methods, and structural mechanics can easily have hundreds of developer-years behind them, including extensive test suites, comprehensive experience with the strengths and limitations of the underlying numerical methods, active user communities, and highly optimized kernels. Since most projects using these codes are driven by application goals (scientific or engineering discovery), any significant restructuring of the codes, let alone a complete rewrite, is extremely risky for application scientists, and often completely infeasible with current software technology. Moreover, many of these applications were not originally written with extensibility in mind and depend on outdated constructs and practices in software design. Thus, their integration into a more modern, modular, component-based software architecture is often more problematic than may initially appear to be the case.

Therefore, it is fair to say that the lowest threshold to early multiphysics results involves treating the coupling variables as I/O to these legacy codes. This approach allows some initial exploration into basic coupled phenomena with very limited modification of the original codes. We also note that higher than first-order accuracy may be achieved within this framework through various iterative techniques, as introduced in Section 3, albeit possibly at a higher performance cost than other choices.

While this approach sounds relatively straightforward, non-trivial implementation issues remain, particularly on HPC architectures, which beget another set of fundamental choices. These typically fall into two general classes.

The service-oriented architecture approach (e.g., the fusion IPS framework (Elwasif et al., 2010)), ensures the maximum independence of each code component by basing coupling on I/O rather than variable/grid-based data formats. The key idea is that interacting with each code at the I/O data level is more flexible and is independent of the internals of the codes to be coupled. Currently file I/O is commonly used for data transfer, but it will be necessary to move toward in-memory transfers for scalability.

A standard for the meaning of the I/O format (e.g., which variables are located where on which mesh type) is still required for the coupling to be numerically consistent, but this approach avoids any direct “contact” with the internal data structures of each component. Note that, while not required, this approach can be naturally carried out in a parallel environment when each component runs as its own executable on a disjoint set of MPI processes. This approach can be implemented on top of a higher-level I/O library.

In the remote procedure call approach, components share information by direct method calls based on a defined abstraction of the information that the components must share. This approach does not imply that the internal data representations of each code are directly modified to accommodate the coupling; rather, more typically, interface routines are used to transmit information between the components for loading into the internal representations of the communicating components. Clearly, the more common the internal representations used by the components, the more efficient, in terms of memory and translation, the information communication process. This approach is more commonly followed with existing HPC codes and is more naturally implemented as a single MPI executable with either disjoint or overlapping processor configurations. When linking independently developed components with different internal data representations, this approach must employ a lingua franca perspective in which the information passed in the calls can be easily extracted from, and placed into, the internal data structures of the individual components. Such an approach was adopted by the ITAPS SciDAC center for dealing with unstructured meshes and the field (solution) information defined over those meshes (Chand et al., 2008; Devine et al., 2009; Ollivier-Gooch et al., 2010).

We note that there is no superior multiphysics software solution in general; the best approach for a particular circumstance depends on application-specific requirements. Usually, though, service-oriented architecture approaches will require much less development work but offer less flexibility and opportunity for performance optimization, in terms of both memory and execution time.

Increasingly, the commercial simulation community is endeavoring to provide “coupling toolkits” to ease the programming burden of building multiphysics simulations. Some of these toolkits, such as MpCCI (2012), claim broad penetration and adoption in the commercial modeling world. Their basic philosophy is to allow users to specify their coupling problem at a high level of abstraction by providing the I/O variables and a description of the orchestration (time sequence of couplings). Such tools are focused on a loose coupling approach, providing the adapter code that can interpolate variables between different simulation packages. For this approach to work, the list of codes that can be coupled is hardwired into the toolkit, something that probably makes more sense in the commercial world where a larger fraction of the community uses a smaller number of codes, and thus it is feasible to support only the most popular ones. Limitations of this approach include little or no accommodation of parallel computing resources, particularly at the very high end, and lack of flexibility in potential candidate code components.

A number of projects have explored implicit coupling based on the JFNK method, introduced in Section 3. This approach requires each of the N c code components to return a solution residual, F i ( u ) , i = 1 , … , N c . The approach can be taken even one step further, with a broader goal of blending the best of both worlds by building in the option to employ an operator-splitting algorithm in cases where a tightly coupled approach is not needed and to employ an operator-splitting solution as a preconditioner when a fully implicit solution is deemed necessary. Section 4.1.3 discusses software to facilitate this approach.

For all these techniques, it is important at the implementation level to categorize the range of coupling types that one encounters across applications (see e.g., Larson, 2009a,b). The frequency and volume of data movement depend greatly on these categories and go a long way in determining the difficulty of implementation or the other strategies described above.

4.1.2 Common needs

From this brief survey of current practices in production multiphysics codes, we identify common needs for HPC multiphysics software to achieve reasonable performance, robustness, portability, and flexibility, as well as to address the expectations of software developers and end-users. These needs fall into several categories that are introduced briefly below and discussed in more depth in Section 4.2. Many of these needs also exist for good-quality HPC uniphysics software. However, while it may be possible to “get away with” some poor software engineering practices in a uniphysics context, addressing these issues is imperative for effective code sharing for multiphysics applications.

Interfaces with appropriate levels of encapsulation and granularity while permitting adequate performance. In designing a multiphysics software package, it is critical to determine interfaces with appropriate levels of encapsulation. If too much is encapsulated, the interface is inflexible, and users cannot access needed data. If too little, objects can easily be used incorrectly, and functionality will not be great enough to encourage adoption. See Section 4.2.1.

4.1.3 Software successes

The development of multiphysics software is becoming increasingly widespread. Several commercial analysis packages contain the word “multiphysics” in their description, and Section 2 mentions some application-specific research software. Below we list a few frameworks, libraries, and other tools that incorporate general-purpose, reusable algorithmic infrastructure in support of parallel PDE-based multiphysics applications, and aim to address the needs identified above. Of course, many additional projects address important aspects of multiphysics simulations. Such projects include approaches that facilitate working at higher levels of abstraction (Sahu et al., 1999; Bulatewicz and Cuny, 2006; Ford et al., 2006; McInnes et al., 2006; DeVito et al., 2011), as well as software such as Cactus (Goodale et al., 2003), Coupler (Liu and Sosonkina, 2011), DDEMA (Michopoulos et al., 2003), FLASH (Dubey et al., 2009), InterComm (Sussman, 2006), Jade (Hegewald et al., 2008), Karma (Mahadevan et al., 2009), MCT (Larson et al., 2005), MUSE (Zwart et al., 2008), Overlink (Grandy, 2004), Overture (Henshaw, 2002), PreCICE (Gatzhammer et al., 2010), Sierra (Stewart and Edwards, 2004), and XpressSpace (Zhang et al., 2011a).

4.2.1 Enabling introduction of new models, algorithms, and data structures

It is particularly important for multiphysics software to attempt to balance the competing goals of interface stability and software reuse with the ability to innovate algorithmically and develop new physical models. Interface stability is important for several reasons:

A related point is the need over time to re-examine interfaces in light of recent advances in mathematics and programming models. For example, providing (possibly optional) access to additional information (e.g., derivatives, not just fields, or other internal states) can inform library design for both couplers and solvers.

We discuss two important aspects of software design that promote interface stability and facilitate experimentation with existing approaches as well as the introduction of new capabilities over the lifetimes of multiphysics codes:

Library interfaces.An important provision for extensibility is the separation of logically (or mathematically) distinct concerns into modules with library interfaces. Because a library cannot make assumptions about program startup or the use of state, it cannot seize control of “main” or assume MPI_COMM_WORLD; see Section 4.2.6 for additional software engineering issues. When preparing an application for use in a multiphysics context, a clear division should be made between the software components responsible for discretization of a physical model and the analysis being performed with the discrete model. The single-physics “analysis” component will typically be bypassed in a multiphysics context, so the interfaces between these components should be thought out carefully, especially with respect to ownership and hidden state.

A key concept is that any state data must be explicitly exchanged through an interface to maintain consistency; examples include common algebraic interfaces, such as vector and matrix classes available via PETSc and Trilinos. This approach may be used to specify inhomogeneous boundary conditions that would otherwise often be computed in a preprocessing stage. It is usually not critical which classes are used by a coupling interface, especially for vectors, because the memory backing a vector from one package can be rewrapped as a vector from a different package without a copy (though dynamic adaptivity and interactions with fast processors such as graphics processing units (GPUs) can cause complications).

Composite operators

Particular to multiphysics coupling is the ability to build block composite operators that aggregate physics applications for fully coupled solution techniques. For example, one could build a composite abstract vector from two physics applications that each supply their own residual vector. Generic tools for building composite block operations take as input abstract operators and produce new abstract operators that can be handed off to solvers.

Also useful is the approach of combining kernels of code to form terms and equations that contribute to the overall multiphysics system residual. These kernels also can supply an approximate Jacobian contribution to be assembled by the coordinating framework or library to be used for preconditioning operations. For example, in MOOSE, these kernels are usually composed of volume and surface integral terms from a variational form of a governing equation, written in residual form.

Another layer of abstraction is the interface to an application that a multiphysics driver requires to support various coupling techniques. Such interfaces can be difficult to define because of the need for balancing design requirements, including flexibility to support solution algorithm requirements, efficiency, and ease of integration. To this end, LIME (Pawlowski et al., 2011) (see Section 4.1.3) provides a high-level domain model for multiphysics coupling, which supports the CASL project (Kothe et al., 2012). This domain model describes the analysis in terms of residuals,

F i ( u ˙ i , u i , { p i , l } , { z i , k } , t ) = 0

transfer operators,

z i , k = r i , k ( { u m } , { p m , n } )

and responses.

g j ( { u ˙ i } , { u i } , { z i , k } , { p m , n } , t ) = 0

Here, F i is the residual vector for physics i ; u i is the vector of unknowns corresponding to physics i ; u ˙ i = ∂ u / ∂ t is the vector of derivatives of unknowns corresponding to physics i with respect to time t ; { p i , l } is the set of all independent parameters for physics i ; { z i , k } is the set of coupling parameters for physics i ; and r i , k is the transfer function k for physics i , where physics i has k = 1 , … , n i k transfer functions. A complete implementation of the domain model above should support all multiphysics coupling algorithms, ranging from loose to tight coupling. Based on the model, if less information is available from an application, the above interface is contracted, thereby limiting the available application coupling algorithms. Five forms of reduced domain models are identified (Pawlowski et al., 2011): the coupling elimination model p → g ( p ) , the state elimination model p → u ( p ) , the fully implicit timestep model f ( u , p ) = 0 , the transient explicitly defined ODE model u ˙ = f ( u , p , t ) , and the transient fully implicit DAE model f ( u ˙ , u , p , t ) = 0 .

4.2.2 Sharing methods and codes among application fields

One of the primary challenges with software development for multiphysics is sharing and integrating algorithms across coupled applications. As codes are integrated into a multiphysics framework, each application may bring in a number of required library dependencies. The reliance on such a large number of libraries can cause difficulties, since software engineering standards and portability can vary dramatically from code to code. Unifying or standardizing common software around libraries with exceptional software practices (domain-driven design, object-oriented design, test-driven design, coverage testing, etc.) and portability will be critical from a maintenance standpoint. Areas where sharing of tools and methods is important include meshing, solvers, data mapping, visualization, I/O, interpolation/mortar methods, and coupling libraries.

The current practice is to develop common components that are released as a framework or library. Two paths are typically followed. The first path is the monolithic Framework (Framework with a capital “F”), which directs and controls all aspects of the problem, including memory management. This approach tends to allow rapid development, since all code is written to concrete implementations and much of the needed utility code is in place. However, the approach tends to lack flexibility because it forces applications to adopt the framework’s data structures and algorithms. The current trend is to move away from monolithic frameworks, since multiphysics coupling requires that applications “play well together”: they cannot insist on being in control of the overall coupled simulation. The second path is to develop frameworks that collaborate (framework with a lower-case “f”). These typically code to interfaces to allow for maximum flexibility. This approach allows users to provide their own concrete implementations of data structures at the cost of being somewhat more complex. See the glossary for additional discussion.

4.2.3 Multiphysics spatial discretizations

Two general strategies are available for spatial integration within multiphysics applications, where the multiple physics models each apply across the entire domain (Hansen and Owen, 2008; Dubcova et al., 2011). The first strategy is based on the use of a single mesh that spans the entire domain. Each physics model is hosted in terms of the basic operators supported on the mesh. The second strategy uses multiple meshes, in the extreme using a unique mesh for each physical model present in the code. This strategy is common when multiple codes are coupled; each code tends to use a mesh tailored to its requirements.

The use of a single mesh is straightforward in concept. However, this is almost always a compromise, since each physical model rarely has similar refinement requirements at a given location in the domain. If a single mesh is used, the mesh must be fine enough everywhere to support the most demanding physical model at that location. Moreover, the developers of the composite application must agree on a meshing strategy that meets the quality requirements and discretization assumptions of all concerned. On the other hand, if these issues can be overcome, this approach has several advantages. First, a single mesh incurs no error in transferring (or interpolating) data between models, since they all are hosted on a single mesh. Second, a single mesh is easier to implement, decompose for parallel execution, and load balance during a transient simulation.

The use of multiple meshes to host a multiphysics problem has nearly opposite considerations. Creating multiple meshes typically magnifies the cost, in person time, of running the application. In some cases, the problem domain may be sufficiently complex that generation of a single mesh is so expensive that the use of multiple meshes is not cost-effective. In addition, as discussed in Section 3.2.1, care must be taken in communicating (transferring) the data between meshes and models in order to avoid sacrificing the accuracy of the final result (Jiao and Heath, 2004a; Johnson et al., 2011).

One option aimed at maintaining the advantages of both approaches is a multimesh assembly strategy (Dubcova et al., 2011) to achieve the spatial diversity advantages of a multiple mesh approach without sacrificing the accuracy of the single mesh approach. For example, Hermes code (Solin et al., 2012) generates multiple hierarchically related meshes using a general version of the finite element method (FEM) known as h p -FEM, where h p -FEM provides for general adaptation and high-order convergence across the composite multiphysics problem. Each of these meshes hosts a physical model that is appropriately fine (in the h p sense) for the desired accuracy of each model. The composite multiphysics model is also refined to meet a global accuracy requirement. The h p -FEM assembly process then integrates across the different models and meshes to provide a coupled solution without using data transfer. Thus, neither data transfer error nor multiphysics coupling error is present in this approach.

In those cases where the individual physics analysis components being linked have independent meshes and internal data representations, a general remote procedure call approach must be taken. Keys to the success of this approach are the abstractions used to define the mesh, solution parameters, and field information, defined over the entities in the mesh. Since the mesh is a discrete representation of the domain of the simulation, it is also desirable that the mesh entities maintain their relationship to the overall definition of the domain. Such relationships are particularly important when adaptive mesh control methods are used and the overall definition of the domain is a solid model as defined in a system for computer-aided design (CAD). The field information, which defines the parameters in the mathematical equations being solved, is best defined in terms of tensor quantities defined over the appropriate geometric entities. In cases where a high-level domain representation is used, the input parameters (loads, material properties, and boundary conditions) can be defined as distributions over the geometric model entities (O’Bara et al., 2002). As part of the discretization process, these fields are mapped onto the mesh. The parameters that are solved for will be defined over the mesh and have known distributions as defined by the equation discretization processes. Maintaining knowledge of the distribution information is key to being able to properly map solution parameters between meshes in a multiphysics simulation.

The fundamental differences in representations used for structured and unstructured meshes indicate that the abstractions best suited for interacting with them are different. In the case of structured meshes, the two key components that must be addressed are the representation of the individual mesh blocks and the interactions between those blocks. In all cases, a compact representation of the mesh blocks in terms of a limited number of spacing and geometric mapping parameters is needed. The representation of the interactions between blocks varies strongly among the various methods (matched mapped blocks, overlapping Cartesian blocks, cut cells, etc.); thus, a separate abstraction tends to be used for each representation.

With unstructured meshes, topological entities from zero to three dimensions, and the adjacencies of those entities (Beall and Shephard, 1997; Ramacle and Shephard, 2003), provide an ideal abstraction for defining the mesh as well as maintaining the relationship of the mesh entities to the geometric domain that is typically defined as a boundary representation. The various ITAPS iMesh implementations (Chand et al., 2008; Devine et al., 2009; Ollivier-Gooch et al., 2010), and others, are based on such a topological abstraction of the mesh. Although it is possible to base the representations used on only the specific mesh entities and adjacencies used in the linkage of specific components, many implementations employ a complete representation in the sense that any of the desired 12 mesh adjacencies can be determined in O(1) time (i.e. are not a function of the number of mesh entities) (Beall and Shephard, 1997; Remacle and Shephard, 2003; Seol and Shephard, 2006).

The parallel execution of multiphysics solution transfer requires extension of the mesh representations that maintain appropriate linkages between the entities in the mesh as they are distributed over the memories of the parallel processors. With structured meshes, the mechanisms for doing this are a strong function of the methods used to understand the interactions between the mesh blocks. For example, the parallel implementation of overlapping mesh blocks will often have one or more individual mesh blocks on a process, along with information on the mesh blocks on other processes with which those blocks interact, whereas with a tree-based decomposition a parallel tree representation might be used (e.g., a parallel octree). In the case of unstructured meshes, the common method of parallel distribution employs mesh patches of neighboring elements defined such that the numbers of mesh entities that share common boundaries with mesh entities on other processes areas are as small as possible. In this case the key component of parallelization is the addition of copies of the common bounding mesh entities on all neighboring processes with links between them (Seol and Shephard, 2006). A mechanism that can be employed to potentially support efficient parallel solution transfer is to “ghost” mesh entities, and their field data, from one process onto neighboring processes that happen to need that information to support solution transfer operations. Since only one process is allowed to change the field data, specific care must be applied to ghosted information to be sure it is up-to-date on any process that will be using it.

4.2.4 Timestep control

Selection of an appropriate timestep for accuracy in a multiphysics calculation can be challenging. Software features with utility for multiphysics calculations span the needs of uniphysics codes, such as robust capabilities for checkpoint and restart and options for timestep control (e.g., manual, programmed, semi-automatic, and automatic timestep selection). When one adds multiple physics models, perhaps in multiple loosely linked modules or separate codes, timestep control across the aggregate application can become very involved. For example, in fission reactor fuel performance simulations, introduced in Section 2.1.2, the timescales of a macroscale structural mechanics model and a coupled mesoscale material model differ by over six orders of magnitude (Tonks et al., 2010). The mesoscale timescale cannot be stepped over, since the evolution of the material properties definitely affects the macroscale simulation. Taking picosecond-size timesteps is not practical for a fuel performance application. Thus one may need to use high-order time integration to support larger mesoscale timesteps, along with subcycling the mesoscale physics. In general, robust and accurate approaches to handling coupled problems with this level of time disparity are still being sought.

In a composite application, one might consider that a timestep is composed of three phases: initializing each package/module in preparation for taking the step, taking the timestep, and then finalizing the step. With multiple modules, the second phase might entail appreciable subiteration between modules, data and results transfer, internal model state updates, and so forth. The third phase includes checking result consistency across the modules, convergence behavior, error metrics, and module intercommunication. The “failure” of a timestep might indicate the inability of subiteration between closely linked modules to converge, or some other error condition. In this case, logic would exist to reduce the timestep or otherwise recover from the error and retake the step. Related algorithmic issues are discussed in Section 3.2.2.

Unfortunately, automatic timestep selection is not robust for general multiphysics simulations and requires tuning for each problem class of interest. Analyst intuition and experience lead to a good understanding of the dynamical timescale of the problem of interest and are valuable in formulating effective approaches and “tuning” them for the application. This observation leads to one of the more effective timestep selection approaches. It is often effective to provide a general interface to a multiphysics code for the manual specification of a timestep history curve that the code should follow, in conjunction with providing residual information as postexecution feedback to the analyst. The analysis can then tune the timestep strategy to provide the desired accuracy and overall simulation behavior, given a few trial simulations.

4.2.5 Achieving performance

While performance is often a major issue in computational science codes, few applications have a quantitative estimate of the achievable performance; hence, they do not know whether the code is using the computing system efficiently. Gaining this understanding is important in judging the impact of different implementation choices and, in particular, the performance trade-offs with different approaches to connecting multiphysics components together.

Achieving this level of understanding requires constructing a performance model good enough to provide an upper bound on performance. Such models are often a combination of instruction, memory reference, and floating-point performance and can identify the critical performance-limiting resource (e.g., Anderson et al., 1999). The “roofline” model, which simply looks at the performance limits due to a number of constraints (each contributing a line to the “roofline”), is a recent representation of this idea.

4.2.6 Software engineering issues for multiphysics integration

Previous sections have discussed some design issues for introducing new models and algorithms into multiphysics applications. In preparing to incorporate an application into a multiphysics environment, a number of issues must be addressed; the following simple suggestions can ease the integration process:

4.3 Difficulties in collaborative multiphysics software

Building on the previous discussion of challenges in multiphysics software design, we now address practical difficulties in collaborative multiphysics software.

4.3.1 Maintenance of a stable and efficient development environment

Managing a multiphysics capability composed of individual application codes from multiple institutions places a number of additional hurdles on coupling infrastructure and application developers. A central issue is how to maintain a stable and efficient development environment in the face of distributed development teams, who may be working under intellectual property, export control, and cybersecurity constraints.

Since individual application codes are under active development at various institutions, procedures must be in place for synchronizing changes to the multiphysics infrastructure. A common procedure is to wait for official releases of each application and then incorporate the changes into the multiphysics infrastructure. Depending on how frequently releases are made, updates based on application release cycles can allow for large windows over which bugs can be introduced and compound. This circumstance makes updating applications difficult and time-consuming, since developers must track down when and how the codes have changed.

A more recent software engineering approach is to use continuous integration (CI), where an automated system continuously checks for “commits” to individual application codes. When a new commit is found, the code is automatically downloaded and tested against the coupling infrastructure. This approach allows an automatic documentation as to which commit breaks a code. Automated emails can be sent to the developer owning the commit. Many bugs can be caught and corrected the day they are introduced, leading to a much more efficient and stable development environment. Such a system requires a rigorous testing system as well as a disciplined development team willing to maintain 100% passing tests so that new failures are not masked by earlier failures.

An advantage of using an automated server for CI is in handling complexities associated with intellectual property and export control. For example, application developers may not be allowed to access the complete coupled code base, yet changes by any application developer may impact the coupled code base. An automated testing system is the first line of defense for notifying developers that their change broke a downstream dependency in code that they cannot access. Cybersecurity concerns can be addressed by locating the CI server outside firewalls.

A recent success story for the CI model is the CASL Nuclear Energy Hub (Kothe et al., 2012). The consortium is composed of 10 core institutions plus an additional 11 contributing partners, each developing their own application codes that are unified under a virtual reactor environment. Because of intellectual property constraints, application developers across the institutions are not allowed access to all application modules. Instead, a three-tiered approach is taken to protect the code base. The first line of defense is a precommit testing script that runs a fast set of unit tests on all codes to which a developer has access. Once this local testing is complete, the commits are pushed to the appropriate institutional repositories. A CI server for the multiphysics infrastructure polls the institutional application repositories every few minutes for changes and will launch a comprehensive test suite that includes all CASL partner codes whenever changes are detected. Results are submitted to a dashboard, and developers are notified by email if their commit introduced a failure. The execution time for tests run for precommit and CI testing is strictly limited so that users can obtain quick feedback that does not impede development efforts. The final tier of defense is a more comprehensive nightly test suite consisting of all precommit/CI tests plus any tests whose execution time is longer than appropriate for CI testing. This is run overnight and takes on the order of hours to complete. This testing can include features such as memory-leak testing (see e.g., valgrind). The procedures developed for CASL are an extension of the Trilinos build and test system.

4.3.2 Investment in software tools and techniques

A significant issue that is exacerbated by multiphysics coupling is the lack of investment in software engineering tools and strong software quality practices across application development teams. The use of software engineering tools (integrated development environments, version control, compilers, build/test and reporting tools) varies widely across institutions and research groups. Multiphysics additionally burdens the performance of the tools because of scaling from a single application to multiple, linked applications. For example, dependency tracking can cause unacceptable configuration and build times if n 2 complexities exist in the configuration system. Such issues might not be noticed in a single application where a smaller number of files and/or dependencies are present.

The lack of adoption of good software engineering tools impacts the development cycle of applications. One issue is that the tool might not be portable to all platforms of interest. Another issue is that a tool may not be funded or supported in the future, leaving users in a difficult position. Investments in open-source software engineering tools should be a priority in order to support HPC applications. An example of a successful project along these lines is the CMake build system (Kitware, 2011). Tools like this are critical, not only for single applications, but especially for multiphysics coupled applications.

Poor software engineering practices can make integration of application components difficult. Some applications may have been developed to run only on a specific system or with a specific compiler with non-standard conforming extensions. Porting to new platforms can be difficult. Memory management, error handling, and global scoping (see Section 4.2.6) can all be controlled with proper software coding guidelines and toolsets (see e.g., Bartlett, 2010).

Multiphysics coupling invariably requires changes to each application code. The stronger the testing in place, the easier refactoring is. Fragile code bases with minimal testing often result in “code stagnation” as developers are wary of making the large changes necessary for integration. The adoption of test-driven development and agile development strategies can address these issues.

4.3.3 Adoption of third-party software

Despite the many advantages of code reuse and success stories of the sort presented in Section 4.1.3, some application developers resist the adoption of third-party software. As introduced in Section 2.2.2, there are technical, programmatic, and personal reasons for this resistance. One is the risk, or perceived risk, of not being in control of the software that is being leveraged, or its long-term availability. Concerns include the long-term maintenance of the third-party software, the viability of the organization producing it, and the ability to influence the development of the software to enhance or maintain its relevance to the application being developed. Concerns may also exist about the “return on investment” of time spent in learning the software; it is not always clear that the time invested in learning the third-party tools will be recouped by increased development efficiency and by reuse in the longer term. Along these lines are concerns about the “opportunity cost”, or perceived cost, of using such tools. Deciding on a particular third-party capability may be hampered by concerns about making a poor selection, particularly if multiple third-party tools appear to meet the user’s requirements. The difficulty in comparing and assessing complex, high-consequence products can result in decision paralysis.

Another concern is the portability of third-party software. The software must span the hardware of interest to the user and must be both supported and tested on these platforms as part of the long-term maintenance strategy.

Another practical consideration is the quality of the third-party software. What processes will be used by the user and supplier to verify and validate the code base being leveraged? If the user makes changes to the third-party software, how do these changes get propagated back to the main repository? How are intellectual-property and export-control issues addressed?

Emotional considerations may also be a factor, such as the resistance to use software that is “not invented here.” Members of the user’s development team might feel security in developing their internal code and insecurity in using third-party capability. Indeed, a developer’s self-worth may become intertwined with the application over time.

4.3.4 Licensing

Additional impediments to using common multiphysics software are (1) the cost of commercial software, especially on machines with large numbers of nodes, and (2) unclear or unacceptable licensing terms for some open-source software. Publicly available open-source software is particularly attractive, and it is consistent with the long tradition within the HPC community. Not only does open-source software remove a cost and inconvenience barrier, but it also allows users to determine exactly what the software does, promotes evolution according to user needs, encourages porting to many platforms, and ensures the long-term availability of widely used software. It is important to use well-known licenses such as the Berkeley Software Distribution (BSD), which confer broad rights to users. Licenses such as the GNU General Public License (GPL), which restrict what users can do with code they have written, are unacceptable to many users, and unfamiliar licenses, regardless of terms, often require lengthy legal consultations at commercial companies and government laboratories.

6.1.1 Leveraging of interfaces

One successful insertion path is that of engineering current software to leverage interfaces. Software infrastructure that encapsulates functionality in terms of well-defined interfaces is able to easily adopt new technologies by swapping capabilities behind those interfaces; see Section 4.2. This design allows a new technology to “look” like the old from the point of view of the application. In this framework, direct comparisons between the old and new technologies can easily be made. In addition, well-defined interfaces often allow the adoption of external libraries, in turn providing the ability to investigate the performance of many methods before choosing a final one for use. Encapsulation of functionality behind good interfaces provides an effective way to transition between technologies that adhere to those interfaces. As a result, more options are available to a given application.

We note that well-defined interfaces can sometimes prohibit extreme performance tuning. For mature technologies, where one technology shows a clear benefit, the code framework may be more closely tied to the technology implementation in order to benefit more from efficiency gains. However, the trade-off often means that interfacing with any further new technologies is much more difficult.

6.1.2 Small side-effort

In many cases, a transition to a new technology starts with a small side-effort. Often this activity is a splinter from a main, funded effort and is conducted by a student or junior staff member who can devote significant time to the project. Within this effort new technologies can be tried in a specific branch or simplified code. If this proof-of-principle effort proves successful, the technology can then be moved into the main code. The principal risk with this mechanism is that since junior staff (typically, temporary or shorter-term) are often the ones conducting most of the work, the people who are most knowledgeable about the new technology may not be present after the initial testing and use.

6.1.3 Full rewrites

An effective way to port new technologies into multiphysics applications is to fully redesign the software considering the demands of the existing multiphysics setting; unfortunately, this is also the most expensive in terms of time required to begin conducting simulations. Additionally, there is no guarantee that for the present project a full rewrite will produce significantly better results than will reusing the existing code. Also, with a full rewrite, it will likely take longer to begin initial simulations.

Even so, fully rewriting software is the best approach if (1) the existing software is out of date by today’s software engineering standards, (2) the mechanism for coupling this software to the multiphysics application is fragile, and (3) someone is willing to pay for a full rewrite. By totally redesigning the technology with a multiphysics mindset, coupling and scalability issues related to data traffic can be addressed effectively. Also, the software can be designed extensibly to take advantage of the most modern algorithms and libraries, as well as to prepare for future evolutions; this approach will help delay any future rewrites and encourage incorporation of new technologies across the community by demonstrating a successful implementation. The value of this approach is in the long-term viability of the multiphysics application achieved by designing it to adapt to existing and coming software practices while still leveraging the knowledge present in the original code. Having the original designer of the technology available is helpful in understanding the underlying mentality of the original design and ensuring that the most effective aspects still appear in the rewrite. We further note that one cannot assume that the new code is as correct as the original; therefore, a full verification is often required when a rewrite takes place.

6.2.1 Clear motivation

The most successful integrations of new technologies into existing simulation codes occur when the new technology or enhancement is clearly needed. However, understanding the goals and hence needs of applications research codes can sometimes be difficult. In addition, assessing the current and potential performance of the code and any new technologies is not straightforward.

A typical research code is often used to analyze or simulate multiple and very different physical phenomena over its lifetime. As a result, the goals of such codes are fluid, and identifying specific requirements can be challenging.

When specific goals and needs are identified, however, a new factor comes into play: the determination of whether a new technology can help the simulation code meet a goal. In order to address this question, performance analyses are required. Multiphysics applications are often large and sometimes difficult to manage, making performance analysis on multiphysics software complex. This complexity, in turn, can make it challenging to quantify the computational goals of multiphysics applications (e.g., memory demands, data traffic, computational speed, order of accuracy). One approach to conducting this analysis is to produce “skeleton codes” that allow for simpler identification and classification of contributing costs. These skeleton codes have the same semantic structure as the complete multiphysics codes but without the full complexity, both computationally and physically.

In order to use this approach, simplified models of the relevant component physics are developed and then included in the skeleton simulation. These skeleton simulations are then profiled, with the goal of identifying the issues that will be present in the full multiphysics simulation and produce clear targets for further progress. These activities invariably help focus the needs of potential new technologies and expose their motivations for use.

6.2.2 Collaboration

One of the most significant impediments to multiphysics applications also produces the most potential for great progress: collaboration between interdisciplinary groups. While the initial implementation of the technology may have been designed only by an application developer, moving to the multiphysics setting, and hence more complex simulations, will likely require the input of software developers to optimize efficiency and the input of mathematicians to study stability and accuracy.

The value of communication goes both ways, though. Math and computer science experts must understand the problems that are involved in the simulation and the assumptions that went into the underlying models. Without this input there is likely no way that full analysis or optimization can be done, and time will be lost while necessary physical assumptions or conditions may be violated during the development. All members of the interdisciplinary development team must understand both their role and the role of the other participants. The most successful collaborations will produce multiphysics simulations that allow for new physical insight obtained with computational efficiency and mathematical accuracy.

Multidisciplinary collaborations also tend to be most successful when each member takes the time to understand the reward structure of other members’ institutions and fields. Authorship order, choice of journals, choice of conferences to present the work, and software releases can all be important factors, but of different weight in different disciplines. It behooves team members to understand these and, to the extent possible, help the members meet their institutional requirements.

6.2.3 Discretization and meshes

The development of new multiphysics capabilities for which solution parameters must be related between physics components requires understanding and interacting with the meshes, as well as distribution of the field parameters over the mesh, as indicated by the equation discretization process. The complexity of developing the needed coupling strongly depends on the degree to which the types of meshes, physics discretizations, and data structures used are, or are not, matched. The insertion options associated with various approaches are outlined in Section 4.2.3.

6.2.4 Software packages

The approaches for addressing multiphysics software challenges given in Section 4 apply to both application-specific projects and to software packages. The burden on package developers is greater, however, because the success of a package depends on its integration into other applications.

Good software engineering alone does not address the reasons given in Section 4.3.3 for resistance to adopting external packages. The developers of software packages need to minimize the real and perceived costs of using their work. Documentation has always been a necessity; and as projects become more complex, that documentation needs to be not just thorough but also well structured. Equally important are examples of the intended use of the package. Examples, particularly multiphysics examples, can demonstrate relevance, give application developers an idea of how to match their codes to the new technology, and drive home the advantages that new technology has over what it is intended to replace.

The chances of an external package being successfully integrated into an application increase if its developers are able to provide support during its use. The developers of multiphysics software often have experience with more varied HPC resources than do their users; hence, their help with migrating to new platforms, troubleshooting, and profiling can be extremely valuable. Support should also include a willingness to adapt software packages to the needs of applications. Just as applications should adopt software packages with new technologies to stay competitive, software packages should incorporate the feature requests of applications to stay relevant.

6.2.5 Funding

Funding is a key element of the successful integration of new technologies. Often, funding in scientific applications areas is provided for use of the simulation to discover new science or conduct an analysis of a given phenomenon. Funding agencies tend to be reluctant to fund new code development without new science. On the algorithms side, funding is usually available only for development of novel methods or software libraries containing algorithm implementations done in a general way so as to benefit multiple application areas.

Clear motivation for new technologies often leads to greater chances of success for their integration. It also leads to greater success in funding that integration, because justifications are more readily made.

One reason science applications and their funding agencies give less priority to adopting new mathematical and computer science technologies is that new technologies can be less robust and hence less reliable. While a new technology may present an opportunity for greater speed or accuracy, the adopting code will need to go through a new verification study. In addition, input files and models may need to be adjusted in order to produce accurate simulations with the new technologies. Even when funding is available for the initial insertion, it is often not available farther down the road for new input files and verifications. Regularly run, detailed verification suites (not solely unit tests) with documented results significantly help in this endeavor. With these suites, new technologies can be verified more quickly, since the full machinery for doing verification with the code is already developed and functioning.

6.2.6 Personnel

Another important consideration for insertion of new technologies is personnel. Often, integration of new technologies requires some software engineering expertise. However, a small- to moderate-sized application effort may not have a dedicated software engineer. Hiring such expertise can be difficult for a short-term position. As a result, scientists often perform large software engineering tasks themselves. However, reward systems at many institutions do not take such work into account. In addition, the work is often not in the scientists’ main line of interest. As a result, software engineering performed in such an environment may not take advantage of some best practices, causing possible delays downstream. Furthermore, a scientist may be less able to tune an algorithm or, when an algorithm fails, determine whether there is a bug in the implementation or an algorithm failure for a problem.

On the mathematics and computer science side, a similar trend is often seen. Again, because institutional reward systems value publishing in disciplinary fields, staff are most interested in implementing and developing new methods, not in placing known methods into a scientific code. This interest is in direct conflict with the scientists’ need for proven and reliable methods to move the science simulation forward. In addition, while some funding exists to package software for use across applications, often package writers are not funded for longer-term support.

A further hindrance to effective cross-disciplinary work from the software developer side is the issue of recognition. Often, software development is seen as a necessary requirement to move a simulation’s capability forward. However, since the main new science is done in modeling and interpreting code results, the software work is not given credit in publications nor in most criteria used for promotions at most institutions. As a result, recruiting for these tasks is difficult, and only rarely have packages and software been supported long-term.

For many of these reasons, implementations often become the task of junior researchers, since they often are the only ones with large enough blocks of time to do complex implementations. These junior staff tend to be less familiar with effective software practices. In addition, junior staff have a high turnover rate and thus cannot be relied on to provide long-term knowledge of and support for the software implementations.