{'fraction_non_alphanumeric': 0.02791605661298194, 'fraction_numerical': 0.00911943108136373, 'mean_word_length': 4.603250253926088, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CUBE: A scalable framework for large-scale industrial simulations

Writing high performance solvers for engineering applications is a delicate task. These codes are often developed on an application to application basis, highly optimized to solve a certain problem. Here, we present our work on developing a general simulation framework for efficient computation of time resolved approximations of complex industrial flow problems - Complex Unified Building cubE method (Cube). To address the challenges of emerging, modern supercomputers, suitable data structures and communication patterns are developed and incorporated into Cube. We use a Cartesian grid together with various immersed boundary methods to accurately capture moving, complex geometries. The asymmetric workload of the immersed boundary is balanced by a predictive dynamic load balancer, and a multithreaded halo-exchange algorithm is employed to efficiently overlap communication with computations. Our work also concerns efficient methods for handling the large amount of data produced by large-scale flow simulations, such as scalable parallel I/O, data compression and in-situ processing.


Introduction
In the past decades, rapid advances in computer architecture and parallel computing, enabling faster and more accurate numerical simulations, have positioned computational fluid dynamics (CFD) as a standard tool in many areas of science and engineering, capable of competing with full-scale wind tunnel experiments. However, despite the remarkable progress in computational power, accurate engineering simulations are still very time consuming and pose several challenges with respect to mesh generation and numerical approximation, in particular in the context of high performance computing (HPC).
CFD simulations of flow phenomena observed in industrial applications, such as aerodynamics of vehicles and aircrafts, flow around wind turbines, etc., require high computational power. Consequently, large simulations of practical interest are not feasible without parallel computing. Parallelization can significantly alleviate the computational burden of such large simulations by distributing the work across many processors. In parallel computing, parallel efficiency and scalability of a numerical method or a flow solver determine how efficiently the computing resources are being utilized. The chosen meshing technique, representation of geometries and numerical methods not only have an impact on the overall accuracy of the simulation, they greatly affect performance and parallel scalability of a solver. This is even more important in today's extreme scale computing environment where simulation frameworks must be capable of harnessing the power of tens of thousands of cores. At this scale, the parallelization strategy of a solver is of uttermost importance, even the slightest amount of workload imbalance, poor communication patterns or serial sections can greatly affect scalability.
In this paper, we present Cube, a framework for efficient computation of time resolved approximations of complex industrial flow problems with moving geometries. The framework is based on a block structured Cartesian meshing technique called Building Cube Method (BCM) [1]. The BCM framework enables easy parallelization and formulation of efficient numerical kernels. A constraint based immersed boundary (IB) method is used to model the flow around complex geometries encountered in industrial applications [2,3,4]. A hybrid MPI+OpenMP based parallelization is used to hide communication cost by overlapped communication patterns, together with a multi-constraint based load balancing framework to enable scalable simulations. Throughout the paper, we will mainly focus on the design and parallelization aspects of the framework.
There are two main approaches that can be categorized under IB methods, namely continuous forcing approach and the discrete forcing approach [5]. The categorization is based on whether the forcing due to an immersed body is applied before or after discretization of the governing fluid flow equations. As the names imply, in the continuous and discrete forcing approaches the forcing is applied before and after the discretization of the governing equations, respectively. An appealing feature of the continuous forcing methods is that they are independent of the discretization employed [6,2,3,7,4]. The discrete forcing methods on the other hand are discretization dependent, but these methods are capable of sharply representing the immersed body interface [8,9,10,11].
In the present work we use a continuous forcing technique, specifically the constraint based IB method, developed previously by Patankar and co-workers [4,2,3], that enables a solver to handle complex immersed bodies without any special treatment during pre-processing. For simulations involving industrial applications, minimization of the overall pre-processing time is important. The overall reduction in pre-processing time depends on whether the IB method employed allows for "non-water-tight" immersed body geometries. A water-tight geometry is one in which the volume inside the immersed body is separated from the volume outside by a closed surface. Industrial scale or production grade CAD data are not always water-tight. Therefore, methods that require water-tight geometries necessitate pre-processing. Thus, we choose the constraint based IB method as it allows non-water-tight geometries, reducing or eliminating the pre-processing of CAD data. Furthermore, the continuous forcing technique coupled with a Lagrangian representation of the immersed geometry is more robust and versatile in handling complex moving geometries compared to discrete forcing methods. A limitation of continuous forcing approach is its inability to sharply resolve fluid-immersed body interfaces. This limitation can be addressed through higher spatial resolution near the fluid-immersed body interface. Thus, we couple the immersed boundary method with BCM so as to achieve high spatial resolutions close to immersed body.
Load balancing is a key factor that strongly influences the efficiency of a parallel code. In a massively parallel computing environment even the slightest workload imbalance can severely affect the performance of a code. Thus, load balancing is an essential aspect of a large scale simulations. It is closely linked to the fundamental aspect of parallel computing namely data decomposition. Data decomposition is typically done either through a pre-processor before the main simulation, or during the initial steps of the simulation. The decomposition is often carried out on the underlying computational mesh that discretizes the computational domain. The aim of such decomposition is to equally distribute the unit of discretization such as tetrahedral or hexahedral cells, blocks across all the workers. Most decomposition techniques assume that the workload of each discretization unit is equal. Although this is true for many applications, there are a large class of problems for which it is not. Lagrangian-Eulerian based immersed boundary methods are one such class of problems. In Lagrangian-Eulerian approach the immersed body is represented by a discrete set of Lagrangian particles which are free to move relative to the fixed background Eulerian mesh. For such a framework, the computational cost of cells or cubes that overlap with Lagrangian particles is different from that of non-overlapping cells. In this work we develop a general multi constraint load balancing technique based on the intelligent remapping approach PLUM [12], to decompose the Lagrangian-Eulerian system.
The outline of the paper is the following; In Section 2 an overview of the numerical methods is given. Our general framework is presented in Section 5 and paralleization, load balancing and I/O strategies in Section 6 together with performance analysis and parallel scalability in Section 7. We present a typical industrial application in Section 8, give conclusions and outline future work in Section 9.

Mathematical Framework
We consider a viscous, incompressible fluid with immersed bodies. The governing equations for such a setup are given by where u is the velocity field, ρ is the density, p is the pressure, and µ is the dynamic viscosity of the fluid. The physical space occupied by the fluid along with the immersed body is denoted by Ω. The domain Ω is divided into Ω f and Ω s to represent the fluid and the immersed body domains, respectively, Ω = Ω s ∪ Ω f . The immersed body, fluid interface is denoted by ∂Ω sf . Note that bold face variables, such as u, are vector quantities, and non-bold faced quantities, such as p, are scalar quantities. A material description is used for immersed body representation, in which the structure resides in a Lagrangian domain Ψ. The physical coordinates of the immersed body, X(c,t), is updated with the immersed body velocity U(c,t), which is defined on the Lagrangian domain Ψ. Here, c is the material coordinate of the immersed body in the Lagrangian domain, c ∈ Ψ, and the physical coordinate, X, of each material point is defined in Ω, X(c,t) ∈ Ω. A material quantity V is related to its Eulerian counterpart v by V = I[v] The operator I maps v from Ω to V in the Lagrangian domain Ψ. A detailed discussion of the Lagrangian-Eulerian framework is presented in Section. 2.1.
Presence of the immersed body in the fluid is modeled by the body force f in the momentum equation. The body force arises due to the rigidity constraints imposed by the immersed body in Ω s , thus by definition f is non-zero in Ω s and zero in Ω f . The constraint force due to the immersed body is evaluated in the Lagrangian domain, defined by F, before it is mapped onto its Eulerian counterpart f . The mapping operator P maps a Lagrangian variable onto its Eulerian counterpart. The operators P and I are defined and discussed in the following section.
The rigidity constraint is imposed in the immersed body domain Ω s and at fluid immersed body interface ∂Ω sf as where u s is the specified velocity of the immersed body, and D is the deformation rate tensor given by Eq. (5) and Eq. (6) physically imply that the rate of deformation of the rigid, flow induced velocities on the immersed body velocity must be zero. In the case where motion or velocities induced by the fluid on the immersed body is not of interest, the above equations are replaced by the following [2] have shown that the rigidity constraint in Eqs. 5 & 6 give rise to a constraint force that is given by λ r , in the equation, is the Lagrange multiplier that enforces rigid motion constraint in the immersed body domain (Ω s ). λ r is analogous to the mechanical pressure, p, which is also a Lagrange multiplier. p is a Lagrange multiplier that imposes the incompressiblity constraint. It is to be noted that we do not directly compute the Lagrange multiplier, rather an equivalent forcing upon discretization.

Lagrangian-Eulerian approach
Accurate representation of intricate details, down to the Eulerian mesh resolution, of immersed body is key to the success of the numerical method. This is especially true when the immersed body is moving. As already introduced in the previous section, a Lagrangian-Eulerian approach is used in this work because a Lagrangian description is a very accurate method of representing complex, mobile immersed bodies. The physical space Ω, Ω s and boundary ∂Ω sf defined in the previous section are based in the Eulerian description. For the Lagrangian description we define a coordinate system Ψ on which the body resides. Ω s is the projection of Lagrangian domain Ψ on to Eulerian domain Ω through the Dirac delta function operator. The operators P and I are defined by the following equations The variable G, which is defined on the Lagrangian domain (Ψ), is projected on to the variable g, which is defined on the Eulerian mesh's immersed body domain (Ω s ). Similarly, Eulerian variable v is interpolated on to the Lagrangian variable V. The coordinate x is defined on the Eulerian domain Ω, x ∈ Ω. And, c is the discrete Lagrangian coordinate defined on Ψ. δ(x) is the Dirac delta function defined as

Building Cube Method
The Eulerian domain (Ω) is discretized using a structured Cartesian meshing technique called the building cube method (BCM), wherein the Eulerian domain is discretized using cubic units called cubes. The cubes are subdivided into fineer cubes in regions of interest, for example around an immersed body, to generate a set of cubes ranging from coarsest level (l = 0) to finest level (l = m − 1) with a total of l n = m levels, and l r = m − 1 number of refinements. Although it possible in theory to represent the cubes in a heirarchical tree structure, in the present version of BCM we do not use a tree data structure to represent the cubes. A simpler, more efficient, data structure, S m a ll t o L a r g e L a r g e t o S m a ll where only the cube coordinates, size, an arbitrary index, number of cells and neighbor adjacency information is used to represent the mesh using linear arrays in this work. An arbitrary refinement ratio can be used between a coarse and fine cube, for the results presented in this work this ratio is restricted to 2. Cubes are further subdivided into cells such that all cubes have the same number of cells. As every cube has equal cells, all cubes are equal and independent units of work. This enables efficient and easy grid partitioning and parallelization. An example of BCM mesh is shown in Fig. 1, the vertical plane slice shows the cubes and the horizontal slice shows the cube subdivided into cells. A cell centered, collocated arrangement is used to discretize the governing equations on the BCM mesh, wherein the Eulerian variables (u, p, and f ) are defined at the center of a cell. The coordinate of the cell center defined by the index i = (i, j, k), in the t th cube, is given by where x t c = (x 1c , x 2c , x 3c ) is the coordinate of the base corner of a cube indexed t, and ∆x t = (∆x t 1 , ∆x t 2 , ∆x t 3 ) is the mesh spacing. Base corner of a cube is the corner with the minimum coordinate along the principle directions.
x c is the set of base corner coordinates of all the cubes in a mesh so that x t c ∈ x c . Similarly, ∆x t ∈ ∆x, where ∆x is the set of mesh spacing of each cube in the BCM mesh.
Discrete stencil operations at a cube's boundary can span across into the neighbouring cube. A halo of cells from neighbouring cubes are defined around each cube to enable discrete stencil operations at cube faces (Fig. 2). The size of the halo is determined by the extent of the widest stencil used in the solver. For the halo exchange between cubes of same size, data is copied from the boundary interior cells of a neighbouring cubes into adjoining halo cells of a given cube. Halo exchange between coarse and fine cube involves data interpolation of data. A schematic representation of halo exchange is shown in Fig. 2. In case of data exchange from fine to coarse cube, the boundary interior cells are interpolated to the center of the halo cell of the coarse cube (Fig. 2b). The interpolation of a cell centered quantity φ from fine cube t to the halo cell of coarse cube T can be expressed as follows Here, i = (i, j) (i = (i, j, k) in 3D) is the cell index of fine cube and I is that of the coarse cube cell, p = (p, q) (p = (p, q, r) in 3D) is the summation index and w is the interpolation weight. Data exchange from coarse cube interior boundary cell to fine cube halo cells is achieved by simply copying the coarse cell data to all the fine cells that lay within the corse cell. The coarse to fine exchange, as depicted in Fig. 2, may be expressed as, Treatment of face centered, staggered quantities is not necessary because we use a collocated arrangement in the present work.

Lagrangian Data Structure
The immersed body is discretized into material/Lagrangian particles that are attached to its surface which constitute the Lagrangian domain (Ψ). An example of the discretization of the immersed body into Lagrangian particles is shown in Fig. 3, which shows the discretization of the surface of a 2D cylinder into Lagrangian particles. Typically, discretization of immersed body surface results in one Lagrangian particle in the Eulerian cell which immersed body surface intersects. In order to enable simple parallelization we use the underlying cube data structure (BCM cubes) of the Eulerian mesh. The Lagrangian points are grouped into an unordered set belonging to each BCM cube. Lagrangian particles, which are inside cube t, are inserted into set(t). Cubes that do not contain any Lagrangian particles will have empty sets. The grouping of Lagrangian particles into sets is demonstrated in Fig. 4, which shows a set assigned to each cube in the mesh shown. In the figure, empty and non-empty sets are highlighted through the use for different colors. Empty sets are colored in black while each of the non-empty set is colored with a different non-black color. The spacing between Lagrangian particles in a cube indexed t is given by ∆c t = (∆c t 1 , ∆c t 2 , ∆c t 3 ), and a given particle is identified through the index R = (R, S, T ). The Eulerian coordinate of a Lagrangian particle is denoted by X t R,S,T . The main purpose of grouping Lagrangian particles into sets is to enable efficient Lagrangian-Eulerian interpolation. Interpolation from Eulerian mesh to Lagrangian points requires, for a given Lagrangian point, information of surrounding Eulerian mesh points. Without the set based Lagrangian data structure, identifying the surrounding points for interpolation is very expensive. This is because the cube indexing in a BCM mesh is arbitrary. As a result, to find surrounding points for interpolation, the cube containing the Lagrangian point in question needs to be searched. Searching is an inherently slow process which is best avoided when possible. The set based data structure eliminates this search operation. The Lagrangian-Eulerian interpolations boils down to algebraic expression for finding the Eulerian indicies i for a given cube t. Typical Lagragian Eulerian interaction is shown in Algorithm 1.    for each cube in BCM mesh do if set is empty then The Lagrangian data structure is cube based, i.e. a set belonging to cube t is independent of the sets of all the other cubes. Therefore, propogation and motion of the Lagrangian particle requires careful design of the set data structure. As the particles of the set of a given cube t move, some of the particles of the set(t) may cross the bounds of the said cube (runaway particles). The runway particles no longer belong the set(t), but they may be within the bounds of one of the many neighbors of cube t. Hence they have to removed from set(t) and added to the set of the neighboring cube. A schematic representation of propogation, insertion and deletion of particles is shown in Fig. 4c. This deletion and insertion of particles from and into sets of cubes, respectively, has to be carried out efficiently. To enable quick insertion and deletion of particles from sets, the set is built on an integer based hash table. Each Lagrangian particle is assigned a unique global integer identifier, which acts as the key to the hash table. Global implies across cubes as well as MPI partitions. And, the uniqueness of the identifier is necessary to avoid duplicity of particles as they are advected across cubes. The average cost of insertion, deletion, or a look up of hash table entries is O(1), and the cost of the worst case is O(N ), where N is the size of the hash table. The frequency of the running into hash table related worst case operations depends on the choice of hash function. Consequently, the choice of a hash function that minimizes hash conflicts is key to the efficiency of the hash table.

Lagrangian-Eulerian Interaction
With the discretization of the Eulerian and Lagrangian domain defined, the projection and interpolation operations between the two domains will be discussed next. The interpolation from Eulerian to Lagrangian domain in discrete form is given by where U t R,S,T is the interpolated velocity at Lagrangian particle indexed R, S, T , and u t i,j,k is the cell centered velocity in the cell (i, j, k) in cube t. Similarly, the discretized version of projection operator in Eq. (7) is given below Here, δ ∆ is the discrete form of the Dirac delta function introduced in Eqs. 7 &8 and it is given by where φ (x/∆x) /∆x is the one dimensional discrete delta function. The following smoothed 3 point function for the one dimensional function φ(r) [13] is used in this work The choice of the type of Dirac delta function is important from an accuracy and a computational efficiency stand point. Wider delta function are known to give better accuracy then narrower ones [14,6,13], but this better accuracy comes at greater computational cost. An n point wide delta function will require n 3 operations to interpolate at one Lagrangian particle. Therefore, from a computational cost stand point it is ideal to choose a delta function that is the narrowest (a 2-point function). A balance between accuracy and computational cost is necessary. [13] have shown that the 3-point wide function in the above equation provides reasonable accuary at lower computational cost compared wider delta functions. Eqs. 9 & 10 can be expressed in a more concise form as respectively. Here, I ∆ and P ∆ are the discrete versions of the interpolation operator I eq. (8) and P eq. (7).

Numerical algorithm
We use a modified fractional step method [15,2] to solve the system of equations in Eq. (1) to (4). In the first sub-step of the fractional step method, a modified momentum equation is solved to obtain a non-divergence free intermediate velocity. The constraint force due to the immersed body is not imposed at this step. In this first substep, the equations may be discretized using any time stepping algorithm. In Cube, Euler, Adams-Bashforth and implicit Crank-Nicolson schemes are available. Here we describe the numerical algorithm using the Crank-Nicolson scheme.
where a superscript indicates time, n corresponds to previous time step and intermediate variables of the fractional step are identified by a tilde, and n + 1 indicates time at the end of current new time step. A and D are discrete operators of the advection and the diffusion terms, respectively. In the next step, the immersed body force is imposed on the intermediate velocityũ to obtain the second where the immersed body force F n , evaluated as the difference between the immersed body velocity U n s and interpolated intermediate velocity

in the Lagrangian domain Ψ is given by
The non-zero divergence of the intermediate velocity,ũ, and that introduced by the immersed body force into u * is removed through the projection step to yield the divergence free velocity field u n+1 where G is the discrete gradient operator. The pressure p n+1 in the projection step is obtained by solving the Poisson equation Eq. (14) is solved with an efficient geometric multigrid solver. Implementing a multigrid solver on a Cartesian grid with discrete forcing immersed body methods (such as ghost cell IB methods [8,9]) is a complex task. In order to retain good convergence, special treatment of the immersed body at every grid level is necessary [16]. Since we use a continuous forcing to represent the immersed body, the effect of the body is only present in the right hand side of Eq. (14). Thus, it will be unaffected by grid coarsening, and only present at the fine grid level. Hence, implementation of the multigrid solver is straightforward, without the need of any special treatment of the immersed bodies at different grid levels of the multigrid solver. Furthermore, we exploit the underlying cube structure, and only create coarser grid levels on a per cube basis. Thus, we can retain the same communication pattern on all grid levels.
Finally, the position of the Lagrangian points/particles are updated by the following equation

Validation
We carry out simulation of two cases to validate implementation of the numerical method. First, we consider the standard benchmark problem of flow around a stationary sphere. Next, to validate the method for moving immersed bodies we consider flow created by an impulsively started plate in a quiescent fluid.

Flow around a sphere
Flow around a sphere is a widely used case to validate three dimensional flows. There are numerous experimental [17] and numerical studies [18,8,9] that have investigated this flow at various Reynolds numbers against which we can validate our simulations. In the present study we carried out numerical simulations for Reynolds numbers ranging from 100 to 1000 (Reynolds number is based on sphere diameter).
Following are the details of the computational domain used for the simulations: −25D to 25D in all three directions with the sphere placed at the center of the domain. Here, D is the diameter of the sphere. The hierarchical mesh of the BCM has l n = 7 and l r = 6. The mesh spacing on the finest level l = 6, where the sphere is placed, is ∆x | l6 = 0.012D. For all the simulations slip boundary condition is used on domain boundaries in y and z directions, and inflow and outflow boundary condition for x − and x + boundaries, respectively. In Fig. 5 we compare the drag coefficient with established numerical and experimental results. We find that our results are in good agreement with both numerical and experimental data.
Table1 shows a comparison of wake bubble measurements for flow at Re = 100 with experimental and numerical data. The measurements of the wake bubble center agrees well with both numerical and experimental data, where as the bubble length shows better agreement with experiments.

Impulsively started plate
We validate the numerical method for a moving immersed body with the simulation of an impulsively started/accelerated infinitesimally thin plate moving perpendicular to its surface in a quiescent fluid. The experimental work of Taneda and Honji [20] is used to validate our work. A plate of dimensions h × 20h is used in a computational domain of size 40h × 40h × 40h, where h is the plate height. The domain extents are: −20h to 20h along all three  [8,9,18,17] at different Re. directions. At t = 0 the plate is located with its centroid at the origin of the computational domain. On the finest level of the adaptive mesh, the mesh resolution was ∆x | l5 = 0.02h. A large aspect ratio of 20 is chosen for the plate so that flow from the plate's lateral edges do not affect the flow at the center of the plate where the flow is analyzed. The domain size is chosen such that the plate is sufficiently far from the domain boundaries so as to avoid any boundary effects. The Plate is moved with a constant velocity U in the x − direction.
The simulation was carried out at Re = 126 and Re = 896, where plate height and velocity are used as characteristic length and velocity, respectively, to define the Reynolds number. A comparison of wake bubble size from the simulations with Taneda and Honji's [20] experimental data is shown in Fig. 6. The figure shows evolution of wake bubble size (L b ) normalized by plate height h as a function of dimensionless time tU/h. We also compare our results with the 2D simulation results of Koumoutsakos' [21] at Re = 126. Our results are in good agreement with both experimental and simulation data from literature.

Software Environment
Based on the Finite Volume Method, we have developed a unified solver framework Cube (Complex Unified Building cubE) for solving large-scale compressible and incompressible flow problems. The framework has a modular design where Cube provides a core library containing kernel functionalities e.g. mesh representation, numerical kernels for computing flow field and I/O routines. Solvers are then developed on top of the kernel by connecting necessary kernel modules together, forming a "solver" pipeline, describing the necessary steps to solve a particular problem.

Object-Oriented abstraction
Written in modern Fortran 2003, Cube employs a lightweight Object-Oriented abstraction. A set of abstract classes/types define canonical components of a solver, which are later provided by a real solver. With a set of predefined solvers and parameters, Cube can be used as a regular CFD solver. Meanwhile, for advanced users it  is also possible to develop their own tailor-made applications by overloading certain kernel components similar to other C++ based solver frameworks, e.g. FEniCS [22] and OpenFOAM [23].
For example, if a user wishes to write her own application specific flux evaluation routine, she would only need to overload the base definition of a numerical flux flux_t (Fig. 7) inside the kernel. Once done, the framework will execute the user's code (in this case the eval subroutine) instead of the one in the kernel every time fluxes are evaluated. This way we can keep Cube's kernel small and general, without application specific code.  6 Enabling Large Scale Simulations

Data Decomposition
Cube has been parallelized using a hybrid MPI + OpenMP approach, where whole cubes are subdivided between MPI ranks, and thread parallelization of the numerical kernels are performed on a per cube basis or with two dimensional slices in the z-direction of each cube. MPI partitioning is performed following a load-balanced linear data distribution [24]; Let P be the number of MPI ranks, p be an MPI rank, N be the global number of cubes, and n be the local number of cubes assigned to a rank. With N = P L + R and 0 ≤ R < P , we have that To minimize data dependencies, the distribution is calculated based on a space-filling curve, the Z-ordering (cf. [25]) of the cubes given by our mesh generator. Using the linear distribution as defined in (15), we can easily calculate an owner of a cube given its global id. Therefore, each MPI ranks builds a look-up table for local to global index mappings, which can later be used to construct adjacency information on the fly. If the Z-ordering is not sufficient, a new distribution can be calculated using ParMETIS [26] or the built in load balancing framework (see Section 6.5).  For a Lagrangian-Eulerian framework, such as the one employed in the present work for modeling immersed geometries, there are more than one methods of decomposing the system. As there are two separate domains, namely the Lagragian domain and Eulerian domain, combined or separate domain decomposition strategies may be employed. Three of the most commonly used approaches for decomposing a Lagrangian-Eulerian system are task parallelism, atomic decomposition, and spatial decomposition [27]. In task parallelism the Lagrangian and Eulerian workload are assigned to separate set of processors. In atomic decomposition, after the Eulerian domain is decomposed, the Lagrangian domain is divided equally among all the processors. Spatial decomposition decomposes the Lagrangian domain on the basis of the Eulerian decomposition. Both task parallelism and atomic decomposition, due to non locality of the Eulerian and Lagrangian domain, necessitate communication between processors for Lagrangian-Eulerian interaction [28]. This is not the case with spatial decomposition. In the present work we employ the spatial decomposition. A schematic of spatial decomposition applied to combined system of the Lagrangian particles and the cubes of the BCM mesh is shown in Fig. 8. In the schematic the partitioning of the BCM mesh is highlighted by the dashed red lines between the cubes, as shown in Fig. 8b. As described in Sec. 2.3, Lagrangian particles are grouped into a set for each cube of the BCM mesh. Lagrangian particle sets of cubes belonging to partition p are assigned to the same partition p.

Halo-exchange
A key to achieve good performance in block structured codes is an efficient halo-exchange routine, locally (between cubes owned by a MPI rank) as well as globally between adjacent MPI ranks. There are many factors that can impact performance, such as the width of the halo or the number of cells per cube [29]. Increasing the number of cells per cube is a less viable approach for Cube. Since we target industrial applications, complex geometries will require very fine resolution in very localized areas of the domain. For a given mesh resolution around a geometry, cubes with greater number of cells would greatly increase the overall mesh size, and quickly exhaust memory resources. Therefore, the most feasible option for our applications is to develop a highly concurrent halo-exchange algorithm, that can handle situations where large amount of halo data has to be exchanged efficiently.
To mitigate this problem, we have developed a highly concurrent multithreaded exchange algorithm, that will overlap packing, sending and unpacking of data with local halo-exchange and other possible work (if the numerical scheme permits). The pseudo-code of the algorithm is given in Algorithm 2, and is divided into two separate functions. When Cube needs to exchange data, the function Exchange is called. The first thread to arrive starts filling communication buffers and posting non-blocking send/recv requests. While the first thread is busy packing and sending data to MPI neighbors, all other threads will perform the exchange of data between cubes local to an MPI rank.
The second function Finalize waits until halo data has arrived, and adds contribution from received data to the local cubes. To increase the concurrency, the routine will process messages on a first come first served basis, by issuing MPI_Testsome calls to pick the first finished recv requests from the message queue.

Overlapped time-stepping
The halo-exchange wait time between Exchange and Finalize until the data has arrived, creates a serial section and severely limits the potential scalability. This overlap window could be utilized to perform operations that do not require information from other MPI ranks. The easiest option to utilize the time between non-blocking send and recive calls is to impose boundary conditions, which does not require any communcation. However, typically, imposition of boundary conditions is not expensive, and will not cover the entire overlap window.
In order to enable utilization of the overlap window in Alg. 2, we subdivide the local cubes on an MPI rank into two different zones, namely internal zone and external zone. Internal zone is a spatial region containing only those cubes all of whose neighbours are on the same rank. And, we define external zone as a region containing cubes with at least one off rank neighbour as shown in Fig. 9a. With such a zoning of cubes, we can extend the overlapped halo-exchange to all stages of the solution algorithm, filling the wait time window between Exchange and Finalize with as much work as the numerical method permits.
Here we describe an overlapped time-stepping method that can make full utilization of the overlap window between halo exchange routines. In the numerical algorithm of the present work (Sec.  locally to compute Eq. 13 in all cubes belonging to the internal zone. After the internal zone calculation is complete, the Finalize function is called to distribute data that has already been received to the relevant external zone cubes. Once the rank neighbor halo information has arrived, the same calculation is performed in the external zone thereby completing substep one (Eq. 13) of the numerical algorithm. Here, an explicit Adams-Bashforth time-stepping was assumed in the description of overlapped time-stepping of the FracTimeStep. Following the completion of the FracTimeStep, the solver can continue to the next substep of the numerical algorithm, the Poisson solver. In the PoissonSolve function, pressure is solved iteratively until the error in pressure converges to a specified tolerance. Within each of these iteration we perform the overlapped halo exchange as described above for FracTimeStep (Algorithm 3). The potential of utilizing the full overlapped time-stepping scheme is illustrated in Fig. 9b, where the time to solution for an incompressible simulation of a full car model (see Fig. 12b) is reduced by almost a factor of half when running on 128 nodes on the K computer with a mesh of 47811 cubes.

I/O strategies
Efficient I/O strategies have in recent years become a very important component of a scalable code. Not only should the I/O routines be able to dump data as fast as possible, but loading input data and restarting from checkpoint files are as important. For large scale simulations, the ability to load/restart from previous runs on any number of MPI ranks without offline processing is also a necessity.
To accommodate these requirements Cube implements parallel I/O in form of MPI I/O. Any kind of data is written to a shared file as a flat binary stream. A small header in the beginning of each file contains necessary information to compute new load-balanced distributions using equation (15), which allows for the data to be read back on any number of ranks automatically.

Data compression
The high spatial resolution required for accurate simulations can easily exceed hundreds of millions of grid points. For such problems, the data size quickly becomes a problem. While quickly writing data to disk is not the biggest concern, the problem is the psychical disk space, which can easily be saturated for a time dependent problem.
To solve this problem, we use the discrete wavelet based lossy data compression algorithm developed by Sakai and co-workers [30]. In this algorithm, a discrete wavelet transform is computed independently for each cube, which we extend to include the halo region. One level of quantization is applied and the entire signal is encoded using Zlib [31]. The result is a stream of compressed data for each cube of various length. All these streams are combined into one large stream, which together with a header data is written to disk as a single file. Since we record the size of each cube's compressed byte stream it is straightforward to uncompress the data on any number of cores call Update Pressure(external cubes) ; end when loading/restarting from already written data. Also, since we use a hybrid parallelization, the compression algorithm is also multithreaded such that each thread will transform and encode independent cubes, and a final reduction will combine all temporary streams from the threads into the stream that will be written to disk. In total, all necessary operations to compress our data account for about one forth of the total I/O time for writing data to disk.

Algorithm 4: Compression algorithm
for i . . . n cube do j ← i + (thrd id − 1); if j <= n cube then dwt stream ← dwt(cube(j)); quant stream ← quantification(dwt stream ); encoded stream ← encode(quant stream ); end i ← i + n thrds ; end The achieved compression ratio depends on at least two factors, the number of cells per cube and information loss introduced by the quantization. For an error tolerance of O(10 −4 ), which is sufficient for a lossy checkpoint restart, the compression ratio ranges between ≈ 1:4 − 1:15 for 4 3 − 16 3 cells per cube (see Fig. 10a), and for higher errors O(10 −2 ) a ratio of up to 1 : 43 can be achieved.
For a typical application, the I/O throughput is on an average around 10GB/s on the K computer (see Fig  10b), which is far away from the K's theoretical peak bandwidth, but well above the average throughput of most HPC applications [32].
However, if data compression is not sufficient, Cube also includes in-situ visualization through VisIt [33], either in traditional interactive mode or in batch mode, where a user can setup a user-defined pipeline that will be rendered during the simulation directly from the simulations memory without any data copying or I/O operations.

Load Balancing
Load balancing is an essential component in today's large scale multiphysics simulations, and with an ever increasing amount of parallelism in modern computer architecture it is essential to reduce even the slightest workload imbalance. An imbalance could severely impact an application's scalability. Traditionally, load balancing is seen as a static problem, closely related to the fundamental problem of parallel computing, namely data decomposition. For a CFD simulation based on BCM, since each cubes contains the same amount of cells the goal is to evenly distribute the cubes among the available cores. However, such a decomposition assumes that the workload for each cube is uniform. For most cubes this is true, but for cubes which contain Lagrangian particles the workload is slightly higher, which implies a workload imbalance. Therefore, to retain good scalability we need to derive a load balancing method that balances the workload not only considering the Eulerian computational mesh, but also the additional workload from the immersed body.

Static Load Balancing
In parallel computing, the idea of data decomposition or static load balancing is simple, namely divide the workload evenly across all the workers. This can be formulated as the partitioning problem.
Given a set of cells C from a domain T , the partitioning problem for p workers can be expressed as, find p with the constraint that the workload: should be approximately equal for all subsets. Solving Eq. (16) can be done in several ways. The least expensive geometric methods, such as space filling curves [25], only depend on the geometry of the domain. These methods are fast, but don't take into account the topology, hence there is a data dependencies between different cells in the domain. For Cartesian meshes such as BCM, neglecting to consider data dependencies will not lead to severe imbalances. In the scenario where all the cells have the same amount of neighbors, if the decomposition method tries to assign cells, which are close to each other, to one worker (in the geometrical sense) then the data dependencies will "automatically" be approximately balanced. However, if cells have a non uniform workload, or the problem has asymmetric data dependencies between cells, we have to resort to graph methods in order to solve Eq. (16).
Graph methods don't solve Eq. (16) directly, instead they consider the k-way partitioning problem. To understand k-way partitioning, consider an undirected graph G = (V, E), with nodes V and edges E. The nodes V are split into k subsets {Q j } k j=1 with the constraint that the number of nodes should be roughly equal in each subset, and the number of edges cut should be minimized. If we model the computational work by V and the data dependencies in the domain by E, we see that this method will balance both the computational work and the dependencies. Furthermore, if we instead consider a weighted graph G and add the constraint that the sum of all weights should be roughly equal in all subsets Q j , the method can then, by allowing different weights in the graph, handle a non uniform workload.

Dynamic Load Balancing
In order to perform dynamic load balancing two components are needed. First, a way to evaluate the workload and second, a way to decompose the data with the constraint to even out the workload. Using graph based methods from Section 6.5.1 we can compute new constrained partitions of our computational domain. But the challenge is to be able to evaluate the current and future workloads, and decide if load balancing is needed. Figure 11: Example of the dual graph of a Building Cube mesh.

Workload Modeling
We model the workload by a weighted dual graph of the underlying Building Cube mesh (see Fig. 11). Let G = (V, E) be the dual graph of the mesh, with nodes V (one for each cube) and edges E (connecting two nodes if their respective cubes share a common face), q be one of the partitions and let w i be the computational work (weights) assigned to the graph. The workload of a partition q ∈ T is then defined as: Let W avg be the average workload and W max be the maximum, then the graph is considered unbalanced if: where κ is the threshold value, determined depending on the problem at hand and/or machine characteristics.
To model a simulation's workload we finally have to assign appropriate values to the graph's weights w i . In order to have a fine grained control over the workload, we let each node have j weights w vj i , representing the computational work for the given node, and we let each edge have k weights w e k i representing communication costs (data dependencies between graph nodes). The total weight for a given graph node is then given by, Algorithm 5: Dynamic load balancing framework.
for each partition q ∈ T do For a typical simulation, we assign the number of grid points in each cube to w v1 i and the size of the halo (number of grid points to exchange between cubes) to w e k i for each of the graph edges connected to node V i . Additional weights can also be added to the edges, but in the present study we limit ourselves to modeling only the halo exchange cost. We add an additional weight to the graph's node to model the additional computational cost of the immersed bodies. The additional immersed body cost w v2 i is added in Eq. (17), which is modeled as a γ · n particles . Here, γ is a cost parameter for the particles and n particles is the number of particles in each cube.
The operations related to the immersed bodies involve computationally intensive operations of interpolation and projection between the Eulerian and the Lagrangian meshes, and is done once every time-step. Unlike w v1 i , which can be modeled as the number of cells in a cube, w v2 i cannot be modeled directly by the number of particles in cube. This is because interpolation and projection of information involves operations in a small box of Eulerian cells (n × n × n, where n depends on the discrete delta function) for each Lagragian particle. γ = n 3 would be an appropriate choice if interpolation and projection operations were carried out inside the pressure solver's iterations and in the Crank-Nicholson scheme's iterative loop. But, the interpolation is carried out once every time-step outside the pressure solver and the temporal integrator. Consequently, it is not clear what value of γ would be optimal. We choose γ in the range of 1 − 4 to investigate an optimal value through a parametric study in sections to follow.
The graph is finally partitioned by a graph partitioner, with the weights as an additional balancing constraint. Once new partitions have been obtained intelligent remapping is used to assign new partitions in such a way that data movement is kept at a minimum. This is achieved by solving the maximally weighted bipartite graph problem (MWBG). A bipartite graph is one in which edges are weighted by the amount of data to be transferred with the old partitions to the right and new partitions to the left. The MWBG problem is solved using a heuristic with a linear runtime [34,35]. The new load balanced partitions are obtained as illustrated in Fig. 15. The entire load balancing framework is given in Alg. 5.

Performance Analysis
To evaluate the performance of the load balancer, we used Cube to solve two different incompressible flow problems on the K computer. And, the total execution time for performing a fixed number of time steps for both an unbalanced (no load balancing) and a balanced case (using load balancing) on various numbers of cores are compared. For both problems we used the QUICK scheme for the convective terms and a geometric multigrid solver for the pressure. Time integration was performed using a second order Crank-Nicholson method. The constraint based immersed boundary method, presented in this work, was used to represent the complex geometries (Fig. 12b &  Fig. 12a) used in the present analysis.
The geometries for the numerical experiments, namely the landing gear and the vehicle, were chosen to represent different types of immersed bodies over which γ can be studied. Although both the geometries are relatively complex, the overall particle density due of the two geometries are different. The particle density of the vehicle geometry is 2.68 particles for every 100 Eulerian cells, while that for the landing gear is 1.16 particles. Here, particle density is defined as the ratio of the number of Lagrangian particles and the number of Eulerian cells. Carrying out the numerical experiments on these two geometries will be helpful in understanding how the cost parameter varies across a range of immersed bodies.

Nose Landing Gear
The first problem is based on the nose landing gear (Fig. 12a) case from AIAA's BANC series of benchmark problems. Our setup uses a mesh consisting of 48255 cubes, subdivided into 16 cells in each axial direction, and the landing gear consists of ≈ 0.5M surface triangles, resulting in ≈ 2.3M particles.
(a) The landing gear geometry.
(b) The full vehicle model. Figure 12: Geometries used during performance analysis.
In Fig. 13 we present the relative time required to perform one time-step. We define relative runtime per time-step (RRPT) as the runtime per time-step (RPT) of a simulation normalized by runtime per time-step of unbalance base case, for each core count. The RRPT measure is helpful in gauging improvement (or the lack of it) of a load balanced simulation with respect to an unbalanced base case. From the results in Fig. 13 we can observe that by using the load balancer the runtime can be reduced to as much as ≈ 60% of the unbalanced case. As the number of cores are increased, the gains of load balancing is less.
The load balancing, or more specifically the estimated workload is also sensitive to the parameter γ which controls how much more expensive it is to perform Lagrangian particle related operations compared to Eulerian operations. We can see that for γ = 1, the lack of runtime improvement in the balanced case is independent of the core count. Except of the improvement at 4096 and 16384 cores, where the balanced case is slightly faster, the balanced and the unbalanced cases consume similar time per time-step. When γ = 2, the balanced case outperforms the unbalanced case at all the core counts except at 256 and 32768 cores, while at 256 the runtime for balanced and unbalanced is almost same, at 32768 the balanced case is slower than the unbalanced case. Unlike γ = 1 & γ = 2, the remaining two cost parameters, γ = 3 & γ = 4, the balanced cases outperform the unbalanced case for all the core counts.
The estimated workload imbalance before and after the use of the load balancer is shown in Fig. 14. The estimated imbalance before the load balancing is termed predicted imbalance, and the estimated imbalance after load balancing is called actual imbalance. The first observation that can be made is that, with the exception of γ = 1, the predicted imbalance is typically greater than the actual imbalance when the load balancer is invoked. The load balancer is invoked when the predicted imbalance is greater than the threshold value of 4%, i.e. κ = 1.04. When γ = 1, for core count up to 2048, the predicted imbalance is less than 4%. As a result, the load balancer is not invoked, thus the actual imbalance and predicted imbalance are equal.
Next, we find that there are a few cases in which the predicted imbalance is less than the actual imbalance, yet the RRPT of these balanced cases is smaller than the base case. In other words, despite there being no reduction in the actual imbalance compared to the predicted imbalance, we see a speed-up in the runtime of the balanced case. Examples of these cases are (γ = 3, core = 256), (γ = 2, = 512), and (γ = 1, core = 4096) in Figs13 &14. The likely reason for this unusual behavior is our choice of the definition of workload imbalance, W max /W avg > κ. This definition of workload imbalance relies entirely on relative difference between maximum workload (W max ) and average workload (W avg ) and makes no consideration of the standard deviation of the workload of the system. For Figure 13: Realtive runtime per time-step for the load balanced and unbalanced cases for the nose landing gear. Figure 14: Comparison of estimated workload imbalance before and after load balancing for the landing gear geometry. The imbalance data before the load balancing is shown in bars filled with white shade and labeled 'Pre.', and imbalance data after load balancing is shown in bars shaded in non white colors and labeled 'Bal.' distributions with large standard deviation in workload distribution prior to load balancing, there is a large gap between maximum workload and minimum workload. And, there could be a large difference between maximum workload and median workload. For such a distribution there is scope for workload optimization through load balancing and overall speed-up in runtime. The data redistribution algorithm used in the load balancer always results in a 'balanced-distribution' whose standard deviation in workload is small. But, it does not guarantee that maximum workload will be close to the average. As a result, there could be situations where the standard deviation in the workload has been reduced by the load balancer, but the maximum workload has not changed much (or may even have increased). Because the standard deviation has reduced, there could be speed-up in the runtime of the balanced case yet the actual imbalance based on our definition of imbalance may remain unchanged.
Lastly, another interesting observation we can make from the imbalance data in Fig. 14 is the trend of the predicted imbalance. The predicted imbalance progressively increases with core count. This trend implies that there is much scope for improvement the runtime of the base case through load balancing as the core count is increased. Yet, in Fig. 13 we see that the speed-up in runtime for the balanced case progressively decreases with increasing core count. This counter intuitive behavior is likely due to the lack of workload model for Lagrangian domain particle communication for imbalance estimation and for data redistribution by the load balancer. When the core count is low, overall number of partitions cutting the Lagrangian domain and the resulting communication cost would be small compared computational cost of Lagrangian operations. As the core count is increased, progressively more and more partitions cut the Lagrangian domain increasing the communication cost while the cost of Lagrangian computations on each partition decreases, which in turn increases relative cost of the Lagrangian communication. Modeling the workload due to communication of Lagrangian particles is important when core count is high. Thus, when the core count is high, where cost of communication dominates (Lagrangian or otherwise), data redistribution by the current load balancer can result in a distribution whose communication pattern may not be very different from that of the unbalanced distribution. As a result, the runtime, which is now dominated by communication costs, of the balanced case may remain unchanged because the workload of Lagrangian communication was not modeled.

Full Vehicle Model
We consider the flow past a full vehicle model (see Fig. 12b) to evaluate the performance of the load balancer for geometries that result in higher particle density. The numerical methods used for this problem are identical to the ones used for the landing gear benchmark. We use a mesh consisting of 38306 cubes with 16 3 cells per cube, and with a vehicle model consisting of ≈ 12.5M surface triangles, resulting in ≈ 4.2M particles.
Relative runtime per time-step data are presented in Fig. 16. It is seen that for all the chosen values of γ, 1 − 4, the load balanced cases perform better than the unbalanced case. For the cases with 256 cores, the runtime reduced to ≈ 40% or less of the unbalanced case when γ was set to 2, 3 or 4. This improvement in the runtime gradually decreases with increasing core count. Unlike the landing gear case, in which a consistent improvement in runtime was seen only for γ = 3 and γ = 4, an improvement is seen for all values of γ for the vehicle case. Although the trend may be a subtle one, it can seen that γ = 3 and γ = 4 more consistently result in a faster runtime compared to the other values of γ. Thus, with regards to the optimal choice of γ, based on the trend observed for the landing gear and the vehicle model cases, γ = 3 or γ = 4 are reasonable first approximations.
The estimated workload imbalance data for the full vehicle case is presented in Fig. 17. It can be seen that the actual imbalance is consistently smaller than predicted imbalance for every all γ and core count. This unlike the landing gear case in which we found that for a few case, the actual imbalance was greater than the predicted imbalance. This difference in trend, with respect to predicted and actual imbalance, between vehicle and landing gear could be due to particle density and volumetric distribution of particles. The particle density, of landing gear is 0.0116 and that of the full vehicle is 0.0268. The particle density of the vehicle is ≈ 2.5 times that of the landing gear case. Owing to a greater volumetric space occupied by the vehicle model (as opposed to the smaller spatial volume occupied by the landing gear), there is greater volumetric spread of the particles relative to the mesh. A consequence of these two factors could be a workload distribution with a small standard deviation from the average workload which ensures that the actual imbalance is always smaller that the predicted imbalance. The trend of predicted imbalance, for all γ, with core count is identical to that of the landing gear case. The reason for this trend is already discussed in detail in the previous subsection.

Parallel Scalability
The parallel scalability and computational efficiency of a software environment depends on the numerical algorithms involved as well as the design and implementation of each component of the software. In previous sections we have discussed the various components of the software, namely halo exchange, overlapped time-stepping, Lagrangian-Eulerian framework, etc., that we have attempted to optimize. While the gain through the optimization of individual components may not be necessarily high, the small gain in each of the components contributes to the overall computational efficiency and parallel scalability of the software. In this section we present results of strong scaling analysis of Cube with and without immersed geometries. In order to do this, we carried out simulations of flow around a Full-vehicle (Fig. 12b). The vehicle geometry was chosen because it is relatively large, dense and has a complex geometry, representative of geometries likely to be used in large scale simulations. Processor cores ranging from 1024 to 65536 were used to carry out the simulations to analyze the cost of the immersed body, and the scalability of Cube. The mesh used for these simulations had 131, 072 cubes, with 16 3 cells in each cube. QUICK scheme was used for the convection term and a geometric multigrid solver for the pressure. The maximum number iterations for the iterative conjugate gradient solver on the coarsest level was fixed at 10 and the maximum number of v cycles were set to 50. A motion was imposed on the vehicle by setting the axial velocity of the vehicle to 1m/s. The magnitude of the vehicle's velocity is inconsequential because for any non-zero velocity of the vehicle, the Lagrangian particle MPI communicator is initiated. The particle communicator is not necessary when the vehicle is stationary. A non-stationary configuration of the vehicle is chosen for the strong scaling analysis because it is the computationally more expensive and severe IB simulation configuration.
Two sets of simulations were carried out, first without any immersed geometry and the second with the immersed geometry. These two sets of simulations will enable us to compare the strong scaling of the Cube with and without the immersed geometry. The chosen geometry and mesh resulted in 5, 097, 567 Lagrangian particles. In this subsection, we shall refer to simulations without immersed geometry as 'channel' simulations (implying channel flow simulations), and the one with the immersed geometry as the 'IB' simulation. In order to trigger the pressure solver in the channel flow simulation a no-slip boundary condition was imposed on the computational domain walls in the y direction. Slip boundary conditions were used along z-direction, and inflow-outflow boundary condition along x-direction. The same set of boundary conditions were used for the simulations with the immersed geometry.
The strong scaling results of the channel flow simulations and the simulations with moving immersed geometry are shown in Fig. 18. Where, the speed up of cube is compared with the ideal speed up compared to the run-time of base case (1024 cores case). Cube shows good scalability and parallel efficiency all the way up to 65536 cores for the both channel and IB cases. Furthermore, both channel and IB show very similar speed-up trend with the core count, which implies that the underlying scalability of the solver is not strongly influenced by the introduction of a geometry into the solver. This can also be seen in Fig. 18b, which is described next. The relative run-time per time-step is plotted in Fig. 18b. We define relative run-time per time-step as the run-time per time-step (RPT) of a simulation scaled by that of the corresponding channel simulation. For example, for any core count, say 4096, RPT of the IB simulation as well as that of channel simulation is normalized by the RPT of channel simulation of the 4096 cores case. This means that relative RPT will always be 1 for the channel case, as seen in Fig. 18b. We can make two observations from this plot. First, scalability of Cube is not strongly influenced by the inclusion of an immersed geometry in the solver. Second, the plot reveals the relative cost of the inclusion of a complex geometry such as a full vehicle geometry. We find that the additional cost of inclusion of the IB varies between 10% and 25%, which on average is around 15% which makes the simulation of large-scale industrial applications in Cube viable.

.1 Flow around a vehicle in turning maneuver
In this section we proceed to demonstrate the capability of Cube to simulate complex geometries undergoing complex motion. For this demonstration, we consider the flow around a vehicle in a turning maneuver. For this case a full-vehicle "dirty" CAD geometry is used. The CAD geometry is not subjected any special surface treatment to fix surface imperfections. All the intricate details of the vehicle geometry, except passenger seats and steering wheel, are retained. The passenger seats were removed to avoid generation of fine mesh inside the vehicle where the flow is not of interest.
The aerodynamic performance and stability a vehicle is strongly influenced by the crosswinds during cruise  and while in turning maneuvers [36,37]. It is difficult to simulate such real-world flow scenarios in wind tunnel experiments [38,39]. Furthermore, it is also difficult to measure unsteady aerodynamic forces in wind tunnel experiments. Thus, it is desirable for numerical methods and frameworks to be able to efficiently and accurately simulate such flow conditions. To this end, here, we present simulation of a vehicle undergoing a turning motion, including wheel rotation and steer, chassis roll, turn and pitch, in a uniform flow. The geometric extents of the vehicle are as follows: wheel base ∼ 2.4 m, ground height ∼ 1.5 m, width in span-wise direction ∼ 1.7 m. The computational domain extended from −25 m to 75 m in the axial direction (x), −25 m to 25 m in the span-wise direction (y) and 0 to 25 m in the vertical direction (z). The mesh size on the finest level near the vehicle surface was ∆x | l7 = 6.1 mm. The mesh used for the simulation is shown in Fig. 1a. Slip boundary condition was imposed on the boundaries along vertical and span-wise directions, and inflow, outflow boundary conditions at x − and x + , respectively. The viscosity and density of the fluid were 1.82 × 10 −5 kg/(s·m) and 1 kg/m 3 , respectively. The vehicle is assumed to be moving at 13.89 m/s (50 kmph) against an oncoming wind of 28.78 m/s (100 kmph). Thus, a uniform flow of 41.67 m/s is specified as the inflow condition at the x + boundary. An outflow boundary condition used at the x + boundary. A slip boundary condition is used on the boundaries along y & z directions. Reynolds number based on the height of the vehicle is 3.435 × 10 6 .
The vehicle motion is divided in time into two types of motion. First, from t = 0 s to t = 0.3 s the vehicle is assumed to be in a rectilinear motion, with a wheel rotation to match the vehicle speed of 13.89 m/s. Next, for t > 0.3 s the turning maneuver is imposed on the vehicle. A vehicle in a turning maneuver involves the following: Rotation of the chassis and wheels about the vertical axis (z-axis), roll of the chassis about the longitudinal axis (x-axis), additional rotation of the front wheels about the vertical axis which causes the vehicle turn. The wheels undergo rotation about the lateral axis (y-axis) when the vehicle is under motion relative to the fluid. The rotation of the vehicle and the wheels are imposed through angular velocities along respective directions. The angular velocity of the chassis is given by ω ch = (ω ch x , ω ch y , ω ch z ) tanh (α(t − 0.3)), the factor tanh (α(t − 0.3)) is used to gradually impose the rotation instead of an impulsively start and α governs the rate of change of the angular velocity. For our simulation we use ω ch = (0.5, 0, 0.025) tanh (20(t − 0.3)) rad/s. Similarly, ω wh = (ω wh x , ω wh y , ω wh z ) is the angular velocity of the wheels. The wheel rotation about the lateral axis is imposed impulsively at the start of the simulation (t = 0 s) with ω wh y = 51.44 rad/s. The angular velocity of the wheel turn of the front wheels is specified as ω wh z = 0.5 tanh (α(t − 0.3)). And, ω wh x = 0 for all the wheels and ω wh z = 0 for the rear wheels. It is to be noted that no rectilinear motion is directly imposed on the vehicle, linear motion of the vehicle is indirectly specified through the inflow condition which results in a relative velocity between the vehicle and the fluid that is equivalent to the relative velocity due to vehicle motion. Only rotational, angular velocities are directly imposed on the vehicle. This configuration is chosen to reduce size of the numerical mesh. By restricting the linear motion  of the vehicle, the fine mesh is generated only in a small region around the vehicle. As the Cube does not have the capability to adaptively refine mesh as the immersed body moves, generating fine mesh around the vehicle path would result in an excessively large mesh. Visualization of the results from the simulation are presented in Figs. 19 &20. The velocity magnitude on a horizontal plane (xy plane) at z = 0.25 m is plotted in Fig. 19. During visualization, the vehicle chassis is box-sliced at y = 0 such that chassis is visible only for y > 0. This enables the visualization of the front left wheel's turning motion, as shown in Fig. 19. It can be seen in the figure that the vehicle and the wheel's turning begins at t = 0.3 s and proceeds until the end of the simulation. As the wheels turn and lose alignment with the chassis (around t = 0.36 s in Fig. 19), a more pronounced separation flow is caused around the wheels compared to relatively less separated flow around the chassis. This is more evident in Fig. 20a -20d, in which iso-surfaces of the vorticity magnitude is shown at four successive instants of the vehicle's turn. The reason for this more pronounced separation near the front left wheel is due to two likely causes. First, as the wheel turns, it looses alignment with the chassis and projects out. Second, interaction between a rotating wheel and the oncoming fluid could interact turbulently leading to the separation. In the sequence of images in Fig. 19 it can be seen that at all instances of the turn, the separation is actually induced by the wheel, not by front bumper.
The forces on the vehicle during the turning maneuver is plotted in Fig. 20e. Forces normalized by the mean axial force,F = F/F avg x , are plotted in the figure. The lateral force F y oscillates about a zero mean until the beginning of the turning motion after which an approximate linear increase with time is seen. The rate of change of the lateral force is as expected. but the lateral force approximately equal in magnitude to the axial force despite the small angle of turn. Furthermore, it is interesting to note that the rate of change of axial and lateral forces during the turn are very different. The rate of change of the axial force is very small and consequently increase in the force due to the turn is very small compared to the increase in the lateral force. Also, the axial force sees a more rapid increase towards the end of the turn (between t = 0.65 s to t = 0.72 s), which implies that the axial force is severely affected when the turn angle is greater than a certain critical value. If the vehicle continues to turn in an axially dominant wind environment, the axial force could increase more rapidly affecting vehicle's stability.

Conclusions
In this work we have presented Cube, a highly scalable framework for large-scale industrial simulations. A Lagrangian-Eulerian based constraint immersed boundary method adapted to the BCM framework was presented. The immersed boundary method is implemented such that it almost entirely eliminates pre-processing of input CAD geometries and allows that usage of 'dirty' CAD data, which was demonstrated in Section8. An Lagrangian data structure, the building cube Lagrangian, that enables efficient Lagrangian-Eulerian interpolation was presented. A spatial decomposition strategy of parallelization for the Lagrangian-Eulerian system in Cube was used in this work.  Various optimization strategies, overlapped halo exchange, overlapped time stepping, I/O strategies, incorporated in Cube were presented. A multi-threaded halo exchange method was presented and a closely related overlapped time-stepping scheme was also presented. Initial test of overlapped time stepping scheme showed a gain by as much as factor 1.6 in runtime. Through MPI I/O we were able to achieve a throughput of 10 GB/s for read performance on the K computer for file sizes of 251 and 618 GB. Through a wavelet based data compression algorithm, a compression ratio of 15 was achieved on Cube. A multi-constraint based load balancing framework for the Lagrangian-Eulerian system was presented. Upto 50% gain in runtime was achieved through the load balancer. Despite the large gains in runtime for most cases, for some cases when core count was large, lack of workload modeling for the Lagrangian communication lead to a significant drop in the gain in speedup. Lastly, a strong scaling analysis of Cube was presented where we show good scalability upto 65536 cores.
Future work includes further optimization of communication patterns, especially on emerging hardware (e.g. massively parallel multicore processors), tuning and improved workload modeling in the load balancing framework, and further development of the numerical methods to extend the framework's multiphysics capabilities.