JSweep: A Patch-centric Data-driven Approach for Parallel Sweeps on Large-scale Meshes
Jie Yan Zhang Yang Aiqing Zhang Zeyao Mo




JSweep: A Patch-centric Data-driven Approach for Parallel Sweeps on Large-scale Meshes
    
Jie Yan^⋆   Zhang Yang^†   Aiqing Zhang^†   Zeyao Mo^†

† Software Center for High Performance Numerical Simulation, CAEP

 Institute of Applied Physics and Computational Mathematics, Beijing, China

⋆ Noah's Ark Lab, Huawei Technologies

Correspondence: {yang_zhang, zhang_aiqing, mo_zeyao}@iapcm.ac.cn

    
=====================================================================================================================================================================================================================================================================================================================







In mesh-based numerical simulations, sweep is an important computation
pattern. During sweeping a mesh, computations on cells are strictly ordered by
data dependencies in given directions. Due to such a serial order,
parallelizing sweep is challenging, especially for unstructured and
deforming structured meshes. Meanwhile, recent high-fidelity
multi-physics simulations of particle transport, including nuclear reactor and
inertial confinement fusion, require sweeps on large scale meshes with
billions of cells and hundreds of directions.

In this paper, we present JSweep, a parallel data-driven computational
framework integrated in the JAxMIN infrastructure. The essential of JSweep is
a general patch-centric data-driven abstraction, coupled with a high
performance runtime system leveraging hybrid parallelism of MPI+threads and
achieving dynamic communication on contemporary multi-core clusters. Built on
JSweep, we implement a representative data-driven algorithm, Sn transport,
featuring optimizations of vertex clustering, multi-level priority strategy
and patch-angle parallelism. Experimental evaluation with two real-world
applications on structured and unstructured meshes respectively, demonstrates
that JSweep can scale to tens of thousands of processor cores with reasonable
parallel efficiency.










§ INTRODUCTION



In mesh-based numerical simulations, sweep is an important computation pattern widely used in solving deterministic Boltzmann Transport Equation (BTE) <cit.>, convection dominated or Navier-Stokes equation <cit.><cit.> and so on. During sweep on a mesh, cells are computed from upwind to downwind in the sweeping direction. One cell can compute only when all of its upwind neighboring cells are computed.


General sweep computation on large-scale meshes is challenging. For rectangular structured meshes where the data dependencies can regular, the well-known Koch-Baker-Alcouffe(KBA) <cit.><cit.> algorithm which uses a pipelining wavefront way to parallelize multiple sweeps has been very successful. However, for the more general deforming structured meshes and unstructured meshes in which data dependencies among cells are irregular, the KBA method doesn't fit and is typically impossible. Instead, a data-driven approach <cit.><cit.><cit.> is often considered. This approach models the cells' data dependencies as a directed acyclic graph (DAG), regardless of the mesh types, then sweep on the mesh is equivalent to a topological traversal on the DAG. Unfortunately, although KBA-based sweep on regular structured meshes has scaled to 10^6 CPU cores and billions of cells in 2015 <cit.>,  sweep on ordinary unstructured meshes still doesn't efficiently scale to 10^5 cores. 

Meanwhile, mesh-based application programming frameworks<cit.><cit.><cit.><cit.> have been increasingly critical to today's complex simulations that require to couple multiple multi-physics procedures. On one hand, multiple simulation procedures developed on a unified framework share the same specification of mesh and data layout, and thus are more consistent to work together. Given the fact that coupling different simulation procedures is difficult for both software development and numerical validation, this really makes sense. On the other hand, by providing users a programming abstraction and ensuring reasonable performance, the framework isolates applications from the evolution of underlying HPC systems, and thus achieves good portability. Recently, in areas related to particle transport in which sweep on mesh is the most time-consuming portion, framework-based coupling of multi-physics simulations have led to several successful cases, including the full core reactor simulation based on MOOSE<cit.> and the ICF (Inertia Confinement Fusion) program LARED-I<cit.> based on JASMIN<cit.>. Nevertheless, so far these cases are still on structured meshes only.

In this paper, we focus on the patch-based mesh application framework, specifically JAxMIN <cit.><cit.> (detailed in section <ref>), where the mesh is divided into patches. Patch is conceptually a subdomain of the mesh. The patch-based approach has advantages on adaptive mesh refinement, mesh management and parallel computation scheduling. Over 50 real-world applications have been implemented on JAxMIN. However, JAxMIN, like most counterparts, adopts BSP (Bulk Synchronous Parallel) <cit.> style of parallel computing, in which all subdomains (patches) first compute using previous data of themselves and other subdomains, and then communicate to update their remote copies. Although being efficient and scalable enough for most numerical solvers, BSP is seriously inefficient for data-driven sweep computation where the parallelism is fine-grained. Furthermore, the fact that a patch often can't finish computation at one time and thus has to compute many times, as well as complex factors impacting performance, makes it hard to realize in JAxMIN's BSP-based abstraction.


We propose JSweep, a patch-centric data-driven approach for parallel sweep computation on both structured and unstructured meshes, embedded in the JAxMIN infrastructure. Fig.<ref> illustrates the JSweep modules in the abstracted layers of JAxMIN. Specifically, our key contributions are as follows:


  * The patch-centric data-driven abstraction, a unified model for general data-driven procedures on both structured and unstructured meshes. The core idea is extending the concept of patch to a logical processing element that supports reentrant computation and communicate with other patches (Sec. <ref>). 

  * The patch-centric data-driven runtime module for contemporary multi-core cluster systems, featuring hybrid parallelism (MPI+threads) and dynamic data delivery (Sec. <ref>).

  * A sweep component based on the above patch-centric approach, enhanced by vertex clustering, multi-level priority strategy, patch-angle parallelism and coarsened graph. (Sec. <ref>).

  * Experimental evaluation with real applications of particle transport on both structured and unstructured meshes, demonstrating JSweep's reasonable performance and scalability on up to 76,800 processor cores (Sec. <ref>).


Besides, we present background and motivation in Sec. <ref>, related work in Sec. <ref>, and finally conclusions in Sec. <ref>.








§ BACKGROUND AND MOTIVATION







 §.§ Preliminaries



Throughout this paper, we use a small set of terminologies as illustrated in Fig. <ref>. Note that we describe them in an abstract view and don't explicitly differentiate the structured and unstructured meshes unless otherwise stated. 



  * mesh/grid: the generic way of describing the discretized domain.

  * cell: the smallest unit of a discretized domain.

  * patch: a collection of contiguous cells.

  * local cells: cells owned by a patch. They are updated by an operator applied to the patch.

  * ghost cells: halo cells surrounding local cells. They are needed for computation but not updated by the local operator. They belong to other patches.







 §.§ JAxMIN: Patch-based Mesh Application Infrastructure






JASMIN <cit.> and JAUMIN <cit.> (denoted as JAxMIN[Abbreviation of J Adaptive Structured/Unstructured Mesh INfrastructure.] for simplicity), are software infrastructures for programming applications of large-scale numerical simulations on structured and unstructured meshes. They share the same design principles, i.e., patch-based mesh management and component-based application programming interface. Although they are different in the way of describing mesh geometry and mesh elements, we omit these details by discussing in an abstract level in this paper.

JAxMIN adopts a patch-based strategy to manage the mesh and data. The computational domain, discretized as mesh, is decomposed into patches. Patch is a well-defined subdomain, that (1) each patch has complete information about its own cells as well as other mesh entities, (2) with ghost cells, each patch can explicitly get all adjacency information about its neighboring patches, and (3) it is abstract enough to hide differences of structured and unstructured meshes.

JAxMIN provides users a rich set of components as the programming interface. Component here is a generic implementation of any computational pattern. To implement a parallel program, users only need to instantiate a component by defining the application-specific computation kernel. So far, general patterns such as initialization, numerical computation, reduction, load rebalance, particle-in-cell communication, are provided. Besides, JAxMIN implements amounts of physics-specific integration components.

Traditionally, JAxMIN adopts the BSP model to organize computations in a component. The computations consist of a sequence of iterations, called super-steps. During a super-step, each patch executes logically in parallel, as follows: (1) does computation independently without data exchange with others, and then (2) does halo exchange communication with neighbors using newly computed data. Since many numerical algorithms fit well in BSP, the patch-based framework has been successful in many areas.

JAxMIN implements a high performance runtime system supporting hybrid MPI+threads parallelism and accelerators, with underlying optimization on memory management, data layout and buffering communication.




 §.§ Data-driven Parallel Sweeps



Without loss of generality, we consider the sweep computation in discrete ordinates (S_n) transport solvers. Sweep is the most computationally intensive portion of source iterative methods solving S_n form of Boltzmann Transport Equation<cit.>. As the name implies, sweep in any ordinate direction requires a computational ordering of cells from upwind to downwind. One cell can begin computing only if all of its upwind neighboring cells are computed.

Parallelizing sweep computation is challenging since it can't be efficiently implemented in a BSP manner. For regular structured meshes, the KBA approach, decomposing 3d meshes in a 2d columnar fashion and pipelining the computation for successive angles, is sufficient with BSP. However, for more general deforming structured meshes and unstructured meshes where data dependencies  among cells are irregular and thus the pipeline can't be easily determined, the KBA approach is almost impossible to implement.

Alternatively, we focus on the data-driven parallelization which is a general approach for sweeps on both structured and unstructured meshes<cit.><cit.>. In this approach, any complex and irregular data dependencies can be explicitly modeled by a directed acyclic graph. As an example, Fig. <ref> illustrates a 2d unstructured mesh and the associative directed graph in a given sweeping direction. Then, the sweep on a mesh is equivalent to a topological traversal on the directed graph, generalized with the user-defined numerical computations on the vertex.



In reality, multiple sweeps in different ordinate directions (angles) that are carried out simultaneously. We can model their induced data dependencies in a single graph and implicitly leverage parallelism of sweeps from all angles.





 §.§ Motivation


Now we consider the data-driven parallel sweep procedure in the context of the patch-based framework, especially JAxMIN. Unlike other numerical algorithms, patch-level data-driven computation can't be naturally supported in BSP, due to difficulties described below. These difficulties motivate us to develop a new patch-centric data-driven abstraction in the next section.



  §.§.§ Partial computation

In data-driven scenarios, partial computation of the patch is essential. As illustrated in Fig. <ref> where one mesh is partitioned into two patches, interleaved data dependencies between the patches means that the patch can't be computed as a whole. In reality, the above zig-zag data dependency can be normal in unstructured meshes. Thus, to be reentrant, partial computation of a patch is necessary.





  §.§.§ Simultaneous sweeps on a patch


In real-world applications, sweeps from multiple directions are often performed in parallel. For example, in the S_2 sweeps example illustrated in Fig. <ref>, one patch would be swept by multiple sweep procedures from 4 different directions. Generally, it is common that some sweeping directions are independent to each other. Thus, to enable such parallelism, simultaneous sweeps on a patch is necessary. In JAxMIN, however, patch is the basic unit of parallel computation, so we need to extend its abstraction.




  §.§.§ Priority strategies

Priority strategies are becoming more important and complex. Previous work <cit.><cit.> have proven that ordering of computing cells (or vertices) is often critical to both parallelism and performance. In their settings, since one process (in MPI) or thread handles only one mesh subdomain, it is sufficient to compute priorities of cells (or vertices in the associated graph). However, in the patch-based framework, one process or thread typically is assigned with arbitrary number of patches, which means patch scheduling is always prioritized than cells within a patch. Thus, we need at least a two-level policy to prioritizing both patches and cells.





§ PATCH-CENTRIC DATA-DRIVEN ABSTRACTION





In this section, we introduce the patch-centric data-driven abstraction for mesh-based parallel computations. Its foundation is the completeness and expressivity of the patch concept in JAxMIN described in section <ref>. In our abstraction, the concept of patch is further extended as a logical processing element being able to compute on itself and communicate with any other patches. Users should follow a think-like-a-patch philosophy to program, and focus on only actions of a single patch, i.e., defining the local computation and inter-patch communication. The abstraction doesn't expose any details of underlying patch execution details. It should be suitable for all patch-based mesh application frameworks including SAMRAI<cit.>, (part-based) PUMI<cit.> and especially JAxMIN.



 §.§ Data-driven patch-programs



Data-driven logics on a patch is encoded as a patch-program. The patch-program is identified by a (patch, task) pair, indicating task is executed on patch. Any data communication between two patches is abstracted as a stream. The stream contains the user-defined data and description of source and dest patch programs. Fig. <ref> presents the interface of patch-program and stream, in which the patch-program is factored into five primitive functions.









We define patch-program fully reentrant to support partial computation (detailed in the next subsection). At the beginning, each patch-programs is set active. And, in later execution, the state of a patch-program transits according to the finite state machine given in Fig. <ref>. If its vote_to_halt function is evaluated true, the patch-program becomes inactive. Once receiving a stream, the patch-program becomes active. Conceptually, if there are no active patch-programs globally, the whole program terminates.

The patch-program, identified by (patch = p, task tag = t), is scheduled to run the semantics in Alg. <ref>, as follows.


  * If runs at the first time, the init function is used to initialize a local context.

  * Receives all streams sent to (p, t) by others, which is processed by the user-defined input function.

  * Calls compute function with user-defined numerical kernels.

  * Sends all output streams to and activates targets.

  * Calls vote_to_halt to evaluate whether there remains ready work to do. If not, deactivates itself. Otherwise, keeps active for being scheduled again.






  §.§.§ Partial computation of patch-program



Partial computation is an essential property of the patch-program. On one hand, generally a patch-program couldn't finish at one time and thus requires many times of scheduling. As illustrated by the S_n sweeps case (Sec <ref>), two patch-programs would depend on data of each other, leading to a dead lock if patch-programs are not reentrant. On the other hand, a patch-program may contain multiple parts of computations that depend on data of different patch-programs, so allowing a patch-program to execute multiple times can benefit from finer grained parallelism.

In our abstraction, partial computation of a patch-program is achieved by the following approaches. First, we allow storing of local context so the state are memorized, as illustrated by the implementation of sweeps in section <ref>. Second, the logics of finite state machine in Fig. <ref> maintains state transition of a patch-program, ensuring the correctness of termination after arbitrary times of partial execution.




  §.§.§ Simultaneous tasks on a patch


Our abstraction supports multiple tasks on the same patch. Since any patch-program is identified by the pair (patch, task), multiple tasks on a patch naturally execute in parallel, even with possible inter-task communications. Whether and how to decompose work on a patch into patch-programs is the programers' decision. For the full S_n sweeps discussed in Sec.<ref>, by defining sweep on any patch p from the angle a as a patch-program (patch=p, task=a), sweeps from all directions execute simultaneously.









 §.§ Scheduling patch-programs


The data-driven engine initializes and continues to schedule active patch-programs to run until program termination.

For general patch-centric data-driven computations, the necessary and sufficient condition of program termination is that globally all (patch, task)s become inactive. To detect the termination condition in distributed situations, general negotiating protocols<cit.> are needed. However, in numerical algorithms requiring the data-driven approach, the workload is known in advance. Thus, we can often detect the termination with little or even no distributed negotiation. For example, in sweeps, the program termination condition is all (cell, angle)s are computed, which is known by every patch before computation, and termination detection only need local information. In JSweep's real implementation, we actually allow the patch-program to commit its remained workload (i.e., number of (cell,angle)s in sweeps) to a data structure shared by the master and worker threads of local runtime system (detailed in next section). The master thread, as representative of the process, participates distributed terminate negotiation only when there are no longer patch-programs with remained workload.



Priority policies are known critical for scheduling computations, yet is tightly coupled with the properties of the problem itself. In section <ref>, we shall discuss several strategies used in parallel S_n sweeps.





§ PATCH-CENTRIC DATA-DRIVEN RUNTIME SYSTEM


In this section, we present the runtime that maps the patch-centric data-driven computation and communication onto underlying resources. Our target platform is the multicore cluster widely adopted in contemporary HPC systems. Fig. <ref> shows an overview of the runtime system. In particular, our design emphasizes fast stream delivery, load balance, fine-grained parallelism and low schedule overhead.






 §.§ Hybrid parallelism


The runtime inherits from JAxMIN a hybrid parallel approach of MPI + threads, in which the program is organized with distributed-memory MPI processes and each process consists of multiple threads. On top of this design, JAxMIN have been highly optimized on both domain decomposition and data management. Besides, to reduce NUMA effect in multi-socket systems, by default our runtime launches one MPI process per processor and bind the process to the processor. Within each MPI process, the master-workers multi-threading mode is adopted. As shown in Fig. <ref>, the master thread is in charge of scheduling patch-programs, communicating streams and detecting global termination, while each worker thread executes patch-programs and communicates only with master thread.




 §.§ Dynamic stream communication


In the patch-centric abstraction (section <ref>), communication conceptually happens between a pair of (patch, task)s and is defined as a stream. For data-driven algorithms, the latency of stream transmission is critical for performance, since only received the dependent data can a patch-program execute. By definition, the communication of stream is asynchronous and dynamic.

The runtime system adopts the routable stream concept and reserves a specific core for master thread to support such timely communication. As defined in section <ref>, the stream itself carries full information of source and target (patch, task)s. By identifying the target patch-program and looking up the route table that maps (patch, task) to (process, thread) , the runtime system can deliver any stream to its target place, either locally or in remote process.


The master thread schedules patch-programs by routing streams. At the beginning, all patch-programs are active and assigned to workers evenly. Later as the execution progresses, some patch-programs become inactive. If the master thread receives a stream whose target patch-program is inactive, it chooses and sets a lightest worker as the patch-program's owner, and then routes it the stream.



 §.§ Distributed progress tracking


The master thread calls the progress tracker to detect global program termination. Once temporarily there is no longer work within the process, the progress tracker is activated. A consensus algorithms is implemented to detect distributed termination<cit.>. Besides the general negotiating protocol, as discussed in section  <ref>, for known data-driven algorithms special condition detection methods would be preferred for efficiency in practice. Currently, we support both.




 §.§ Known Issues

Currently, we adopt an one-time assignment policy on scheduling newly ready (patch, task) and one (patch, task) may require multiple times to finish its work, such that priority inversion among workers is possible. In this situation, one worker may hold many high priority tasks that leads to serious waiting of other workers. In the worst case, although in practice rarely happen, it may make the whole execution serial.





§ NEW PARALLEL ALGORITHM OF SN SWEEPS


Now we describe a new parallel sweeps algorithm based on JSweep, the above patch-centric data-driven approach. Further, we explore four optimizations, including scheduling by (patch,angle), vertex clustering, multi-level priority strategy, and coarsened graph, which are natural and efficient to implement thanks to the expressibility of the patch-centric abstraction.



 §.§ Patch-Program Implementation


We assume that the mesh has been decomposed into patches with general spacial domain decomposition methods (for example, the METIS <cit.> and Chaco <cit.> for unstructured meshes, Morton and Hilbert space filling curves for structured meshes). Each process is assigned with an arbitrary number of patches, shared by all its threads.

Formally, we define the directed graph induced by sweeping meshes as G = (V, E), where each vertex is a (cell, angle) pair, and each edge is directed data dependency between two vertices. An edge (u, v) means vertex v depends on vertex u's data. For any patch p and a sweeping direction t, we denote the induced subgraph as G_p,t = (V_p,t, E_p,t), where V_p,t is the set of vertices (i.e., { (cell, t) }) and E_p,t is the set of edges.

Listing <ref> presents the patch-centric implementation of parallel sweeps. As presented, the patch-program consists of two parts, i.e., local context and interface implementation.

language=c++, basicstyle=, tabsize=2, frame=lines, numbers=left, numberstyle=, numbersep=1pt, numberblanklines=false, belowcaptionskip=4pt, xleftmargin=4pt, mathescape=true
emph=init, input, output, compute, vote_to_halt, emphstyle=, emph=[2]solve, emphstyle=[2]
[caption=Patch-Program of data-driven parallel sweeps,label=patch-sweep, firstnumber=1]
//G_p,t = (V_p,t, E_p,t) is the subgraph of patch p with
//task tag (i.e., sweeping angle) t.
SweepPatchProgram(PatchID p, TaskTag t)

	// Part 1: Local Context
	int counts[|V_p, t|];
	PriorityQueue Q;
	Map<Pair<PatchID,TaskTag>,Stream> outstreams;

	// Part 2: Interface Implementation
	void init() 
		Q.clear();
		for(each Vertex v in V_p, t) 
			counts[v] = #. v's upwind vertices;
			if(counts[v]==0) Q.enqueue(v);
		
	
	void input(Stream s) 
		while((edge(u, v),data(u)) = s.read()) 
			counts[v] = counts[v]-1;
			if(count[v]==0) Q.enqueue(v);
		
	
	void compute() 
		Vector<Vertex> vertices;
		// N is the vertex clustering grain
		while(!Q.empty() and vertices.size()<N) 
			Vertex v = Q.dequeue();
			vertices.push_back(v);
			for(each v's downwind Vertex w) 
				if(w is in V_p, t) 
					counts[w] = counts[w]-1;
					if(counts[w]==0) Q.enqueue(w);
				 else 
					Stream    s = outstreams(patch(w),t);
					s.write(edge(v,w), data(v));
				
			
		
		solve(vertices); //user-defined computation
	
	Stream output() 
		Stream s = fetch(outstreams);
		return s;
	
	bool vote_to_halt() 
		return Q.empty();
	



The local context contains all necessary states required by a reentrant sweep on the patch, including: (line 6) an array of counters that count the number of unfinished neighbors for each local vertex, (line 7) a priority queue storing ready vertices, and (line 8) streams later sent to other patch-programs.

The interface functions implement DAG-based data-driven sweeps on the patch p in the direction t. The init function initializes each vertex's count variable to the number of its upwind neighbors, and collect source vertices into the ready queue Q. The input function receives data of vertices from remote patches, updates counts of related local vertices; once a local vertex's count decreases to zero, put it to the ready queue. The compute function collect a sequence of ready vertices and computes on them with user-defined numerical computation, updates their downwind neighboring vertices. The output function generates streams sent to remote patch-programs. The vote_to_halt function evaluates whether the patch-program should deactivate.



 §.§ Optimization: Patch-Angle Parallelism


JSweep naturally supports simultaneous sweeps on a patch, from different angles (i.e., sweeping directions). As shown in Listing <ref>, we achieve this by setting task tag of the patch-program to the id of the sweeping direction. Consider the example in Fig. <ref> again, in which full S_2 transport sweeps are carried on a 2d unstructured mesh of 2 patches. In this example, sweeps from different angular directions are independent, and thus patch-angle parallelism can be fully enabled. This is especially useful for small meshes with large number of angles.


In particular, lets's focus on the colored patch (say p) and sweeping directions. The sweeps in Ω 1 and Ω 2 directions are from the same corner and induce same data dependencies between cells. In full sweeps, there are no dependencies between sweeping directions, thus patch-program(p, Ω 1) and patch-program(p, Ω 2) can be execute in parallel, which could reduce the idle time of computing resources (processor cores).

For sweeps Schedule by patch is not always efficient, since for a patch there is no schedule locality among different sweeping angles. Actually, the interval between starting time of two angles maybe very long, especially when the grid scale is large.




 §.§ Optimization: Vertex Clustering


We adopt vertex clustering in the patch-program. As shown in Listing <ref>, the compute function collects and computes on multiple ready vertices, rather than a single vertex. For example, in Fig. <ref> vertices in the same dashed rectangular are clustered together.

Benefits of this optimization is two-fold. On one hand, vertex clustering can dramatically reduce the scheduling overhead by reducing execution times of a patch-program. For example, in Fig. <ref> the patch 1 and patch 2 need only two and three executions respectively, compared to eight times of no clustering. On the other hand, vertex clustering aggregates multiple streams to the same target into a single stream and thus reduces communication overhead. For example, in Fig. <ref>, the inter-vertex communications of 8 → 9 and 12 →14 can be combined into one message.

Nevertheless, we need to choose the clustering grain carefully. While reducing overhead of schedule and communication, vertex clustering also has potentially negative effect on parallelism since it may defer the communication thus delay scheduling of other patch-programs. Excessive clustering can lead to long communication delay and thus longer execution time. To illustrate this, consider SnSweep-S, an example in JAxMIN package, which implements a Sn solver for neutron transport equations on 3d structured meshes. The experimental results (mesh cells: 160 × 160 × 180, patch size: 20 × 20 × 20, S_2 ordinates, 8*12 CPU cores) are shown in Fig. <ref>.









 §.§ Optimization: Priority Strategy


We adopt a two-level hierarchical priority strategy, i.e., (patch, angle) priority and vertex priority.

The (patch, angle) priority is used for JSweep runtime to schedule patch-programs. For patch p and angle a, its priority is calculated by the following formula:



prior(p, a) = prior(a)*C + prior(p)


, where C is a constant factor. In S_n sweep, to avoid waiting of downwind patches, we want patch-programs with the same angle are continuously scheduled to execute such that the data streams are delivered to the nearby patches as quickly as possible. Thus we set the importance of prior(a) always higher than prior(p) in the formula, by multiplying a large factor C over prior(a). Meanwhile, with respect to prior(p), however we can't reduce a single objective. On one hand, we hope the upwind patches are computed as earlier as possible such that more parallelisms are available. On the other hand, we also want the patches neignbouring other unfinished patches are computed earlier, but these preferred patches hare typically on the downwind of a sweeping direction. Based on the first objective, we develop two priority strategies: LDCP (Longest Distance on Critical Path) for structured meshes and BFS (Breadth First Search) for unstructured meshes. Based on the second objective, we develop the priority strategy SLBD (Shortest Local Boundary Distance, a DFS variant that prefers vertices most close to patch boundary) for both structured and unstructured meshes.

The vertex priority is used within a patch-program to order local ready vertices in PriorityQueue Q in Listing.<ref>. As prior(p), vertex priority also has to trade off more parallelism and earlier communication. The strategies proposed for prior(p), i.e., BFS, LDCP (for structured meshes only) and SLBD, are also suitable for vertex priority. In practice, however we observed that SLBD performs constantly best for unstructured and especially unstructured meshes, as illustrated by SnSweep experiments in Fig.<ref>.






 §.§ Extra Optimization: Coarsened Graph



Coarsened graph, not presented in Listing.<ref> for length limit, can be treated as an extension to vertex clustering. In reality, the mesh structure and its data dependencies are always constant in most or even all sweeping iterations. Thus we can cache the vertex clustering results to build a reusable coarsened graph. For example, in Fig. <ref>, the directed graph (left) is transformed into a much smaller coarsened graph (right) according to the previous clustering results in Fig. <ref>.



Formally, we define a coarsened graph as the property graph: 𝒞G = (𝒞V, 𝒞E, P(𝒞V), P(𝒞E)) where 𝒞V is set of coarsened vertices derived and 𝒞E the set of coarsened edges. 𝒞G is the topology of vertex clusters and the directed communication relationships of vertex clusters. Property of a coarsened vertex cv ∈𝒞V, P(cv), is the series of corresponding clustered DAG vertices, while property of a coarsened vertex ce ∈𝒞E, P(ce), is the combined edges of source vertices and target vertices in DAG. For example, in Fig. <ref>, P(cv2) = (9,14,10), P(cv4)=(7,8,12,13,15), P(cv4 → cv2)=({8,12}, {9,14}). Since 𝒞G is the task graph generated in the scheduling process, we have the computability theorem:


If a directed graph G is acyclic, its derived coarsened graph 𝒞G is also acyclic.


Besides, 𝒞G is distributed at the beginning of its construction. We implement it in same technologies presented in agent-graph <cit.>. With coarsened graph, sweep is carried on DAG in the first iteration and on 𝒞G in all subsequent iterations until the mesh changes. In our practice with JSNT-S<cit.>, the cost of building 𝒞G is less the one DAG-based sweep iteration itself while the speedup of sweeps on 𝒞G over DAG can be 7-10 folds.






 §.§ Application Interface


Parallel sweeps has been implemented as a component in JAxMIN. Like other JAxMIN components, its application interface includes the following two parts.


  * patch strategy class. Users need to define it following a class template which includes three primary member functions: (1) calculate data dependencies of local cells, (2) register patch data variables,  and (3) sweep on a given set of (cell, angle) pairs. By defining this class, users can implement desired numeric computation, such as flux sweeping in neutron and radiation transport.

  * sweep integrator component class. It requires the above user-defined patch strategy as a parameter. In the application, users build a sweep component, and calls its member function to solve the sweeping computation in the mesh.


The SnSweep examples in JAxMIN release packages illustrate how to define and use sweep and other components to implement source iterative method solving neutron transport equations on 3d structured and unstructured meshes.









§ EVALUATION



Platform  All experiments were carried on Tianhe-II, the world's fastest supercomputer in 2015<cit.>. We use at most 3200 nodes. Each node has two Intel Xeon E5-2692v2 12-core processors, equipped with 64GB memory and Tianhe-Express-II network of 40GB/s bandwidth. The operating system is Kylin Linux. All applications are compiled with Intel C Compilers (icc13) and customized MPICH2.

Applications  We use two real JAxMIN-based Sn applications, JSNT-S <cit.> and JSNT-U <cit.>, to investigate the efficiency of JSweep on structured and unstructured meshes respectively. The used meshes are visualized in Fig. <ref>.











 §.§ Evaluation on structured meshes


JSNT-S <cit.> is a JASMIN-based S_n package for structured meshes, which implements most functionalities of TORT <cit.>. We use the well-known Kobayashi benchmark to evaluate JSweep. In particular, we focus on the strong scalability. In all the following experiments, JSweep is configured as follows: patch size = 20 × 20 × 20, vertex clustering grain = 1000,  and the priority strategy is SLBD+SLBD.

We first evaluate JSweep with the original Kobayashi benchmark (Kobayashi-400). It solves the single energy group S_n transport equations with scattering, on a cubic mesh (Fig.<ref>) of 400 × 400 × 400 cells with 320 angular directions. As presented in Fig.<ref>, with increasing number of cores, JSweep shows reasonable scalability constantly, with a speedup 14.3 (or parallel efficiency 44.7%) on 24,576 cores compared to 768 cores. 







We then evaluate JSweep on more CPU cores with a larger problem by modifying the mesh of Kobayashi input to 800 × 800 × 800 proportionally, namely Kobayashi-800. As shown in Fig. <ref>, JSweep scales to 76,800 cores with a reasonable speedup 7.4 (or parallel efficiency 46.3%), normalized to performance on 4,800 cores.


Note that with respect to weak scalability whose results is not shown here, however, parallel efficiency of JSweep is not as good as KBA-based algorithms as reported in other literatures <cit.><cit.>. We conclude the reasons to two aspects. First, rather than the KBA's 2d columnar decomposition, we use a general domain decomposition methods which minimizes the volume-surface rate rather than optimize for pipelining sweeps. Second, the data-driven algorithm introduces overhead of DAG operations that is less scalable. For the second problem, we have proposed a DAG coarsening algorithm in another work <cit.>, which dramatically reduces the DAG scale and improves the sweeps performance. 





 §.§ Evaluation on unstructured meshes


JSNT-U<cit.> is a JAUMIN-based S_n package for unstructured meshes, primarily used for numerical simulations in high energy physics. We evaluate JSweep on two shapes of unstructured meshes, reactor core (Fig.<ref>) and ball (Fig.<ref>). Unless otherwise stated, default configurations of experiments are as follows: priority strategy SLBD+SLBD, patch size = 500 cells, vertex clustering grain = 64, #angles = 24 (S_4) and #energy groups = 4.




  §.§.§ Hyper-parameters' effect to performance

We change and investigate three hyper parameters respectively in order, i.e., patch size, vertex clustering grain and priority strategy, while keeping others default. As shown in Fig.<ref> (left), with increasing patch sizes (i.e., #cells of a patch), the runtime first decreases quickly since the larger patch size reduces total communication between patches, and then slightly increases since the larger patch size also leads to longer waiting time of downwind patches whose execution is driven by data from this patch. Fig.<ref> (right) shows the effect of maximum vertex clustering grains. With increasing vertex clustering grain, the runtime decreases quickly and then keeps steady. Unlike on structured meshes (Fig.<ref>), however, the runtime no longer increases with a large clustering grain. By profiling, we found that the actual number of available vertices is between 16 and 64 at most time, which means the real clustering grain is limited by parallelism. With respect to priority strategies, as shown in Fig.<ref>, their effect to performance is not so significant as that on structured meshes (Fig. <ref>).








  §.§.§ Strong Scalability

As shown in Fig. <ref>, JSweep performs good strong scalability on both small and large meshes. For the small scale problem (ball of 482,248 cells), JSweep shows a speedup of 11.5 (parallel efficiency 72%) at 384 cores and goes to a speedup of 75.8 (parallel efficiency 30%) at 6,144 cores, normalized to the 24-core base. For the large scale problem (ball of 173,197,768 cells), JSweep shows a speedup of 9.9 (parallel efficiency of 62%) at 49,152 cores, normalized to the 3,072-core base. Given that our tested mesh is a ball constructed with tetrahedrons, the above scalability should be reasonably good.






  §.§.§ Weak Scalability

Fig. <ref> presents the results of weak scalability evaluated on the ball (originally 482,248 cells) and reactor core (originally 64,479 cells) meshes. In particular, mesh size is increased in a normal approximate refinement method. As shown, the weak scalability of JSweep is not good enough, although reasonable. For reactor, the parallel efficiency at 12,288 cores is about 40%, while for ball it is lower than 20%. One possible reason is that in JAxMIN, the original small mesh is first partitioned and distributed to processes, and then each process refines the assigned subdomain, leading to thick subdomains that dramatically increase length of critical path in the sweeping direction.





 §.§ Runtime Overhead Analysis


JSweep employs a runtime-based approach, thus the overhead is essential for performance. To investigate the overhead, we carry out a detailed profiling of JSNT-S on small scale Kobayashi benchmark. In particular, the problem has a 200 × 200 × 200 mesh. All optimizations are enabled and all hyper parameters are the same with that in Sec.<ref>. We present one sweep iteration using coarsened graph.

Fig.<ref> shows the time breakup in a strong-scaling fashion. The overhead introduced by JSweep (i.e., the graph-op and pack/unpack) is moderately low (approx. 23%), and the major performance loss comes from idling of CPU cores (22%-46%). Communication takes 13%-19% the total time. With more deep optimization and advanced priority strategy, we expect to lower both the overhead and the idle time.







 §.§ Performance Comparison with other systems



We first compare JSweep's performance with previous JAxMIN, which already implement efficient algorithmic optimizations  <cit.> and achieve good performance. Nevertheless, we show that with innovations on data-driven abstraction and runtime system design, JSweep outperforms them on both structured and unstructured meshes.

Fig.<ref> presents results of JSweep and JASMIN-based SnSweep program (a data-driven implementation of Sweep3D<cit.>). We choose SnSweep because it has been optimized manually with all techniques introduced in Sec.<ref>, including a coarsened graph variant which caches the vertex clusters and their communication relationships by MPI tags. As shown, JSweep's runtime is constantly less than JASMIN.

Fig.<ref> presents runtime comparison of JSweep and JAUMIN-based JSNT-U. Again, JSweep shows constant runtime reduction, and with increasing number of cores the comparative advantage becomes slightly bigger.



Besides, Table-<ref> compares parallel efficiency of JSweep with other work in literatures. We can see that for Kobayashi problem, JSweep demonstrates comparable scalability with Denovo's KBA-based implementation. For unstructured ball (sphere) mesh of tetrahedrons, JSweep scales worse than the manually implemented data-driven algorithm PSD-b<cit.>. However, note that JSweep is a solution of general framework. Due to the lack of common public problems and availability of the systems, it is difficult to compare frameworks directly. 








§ RELATED WORK



The idea of patch-centric abstraction is partly inspired by the vertex-centric models <cit.> in graph-parallel frameworks <cit.><cit.>. In a vertex-centric model, user defines a vertex-program for a single vertex and the framework lifts the vertex computation to the whole graph, conceptually in parallel. However, unlike vertex in graph, patch is not the basic element of mesh, which means patch-centric data-driven abstraction has fundamental difficulties including partial computation, priority inversion and multi-tasks on a single patch. In this paper we comprehensively addresses these issues and formalize a general patch-centric data-driven approach. In fact, the patch-centric abstraction can be seen as a straightforward extension to existing patch-based frameworks (see a survey in <cit.>), such as SAMRAI <cit.>, (part-based) PUMI <cit.> and especially JAxMIN<cit.><cit.>, in which the mesh is decomposed into and managed by patches. 


Task-based programming models, such as PaRSEC <cit.> and more general Charm++ <cit.>, are also suitable to implement location-based data-driven computations. For example, a recent work <cit.> has implemented a PaRSEC-based S_n sweep solver on 3d cartesian meshes, demonstrating high efficiency with 34% of peak performance at 384 cores. Compared to JSweep, however, these task-based frameworks are not well-abstracted for mesh-specific parallel computation and thus require users to remap many conceptions.



The most important application of mesh sweeps is S_n transport sweeps. It is the critical kernel of source iteration method on solving discrete ordinates form of Boltzmann transport equations. For regular structured meshes, Koch-Baker-Alcouffe(KBA) <cit.><cit.> has been very successful. KBA decomposes the 3d mesh in a 2-D columnar fashion and pipelines computation by angles in a way of wavefront, which synchronizes the inter-processor communication. KBA-based multigroup radiation or neutron transport codes, including ASCI Sweep3D <cit.>, Ardra <cit.> and Denovo <cit.>, has scaled to O(10^5) cores with parallel efficiency of up to >90%. For non-confirming structured meshes and unstructured meshes, however, parallelizing S_n Sweeps is challenging because of the irregular data dependency of computation on cells. Pautz <cit.>, Plimpton <cit.> and Mo <cit.> proposed a graph-based parallel pipeline algorithm of S_n sweeps on unstructured meshes respectively, in which data dependencies between cells and angles are explicitly modeled as a DAG and thus sweeps on mesh is equivalent to topological traversal on the DAG. Recently, further with domain overloading strategy, SCEPTRE transport code <cit.><cit.> can scale to at most 68,120 cores with reasonable parallel efficiency.

What distinguishes our work from previous efforts is that JSweeps is implemented in a general context of patch-based application programming frameworks. For the demand of complex multiphysics simulations, application frameworks have become increasingly important, since numerical methods implemented on the same framework have natural advantages on coupling. Successful cases related to particle transport include the full core reactor simulation based on MOOSE <cit.> and the ICF (Inertia Confinement Fusion <cit.>) program LARED-I <cit.> based on JASMIN. Besides sweeps, our patch-centric data-driven approach can efficiently support other data-driven algorithms such as particle trace, and thus is a necessary supplement to existed BSP approach.





§ CONCLUSIONS AND FUTURE WORK



We have presented JSweep, a generic patch-centric data-driven framework integrated in the JAxMIN infrastructure. In particular, we propose the patch-centric data-driven abstraction whose essential idea is extending the concept patch as a logical processing element that is fully reentrant. Also, our abstraction supports multiple tasks on a single patch and arbitrary patch priority strategies. Further, targeting contemporary HPC systems of multicore cluster architecture, we implemented a high performance runtime system to map the patch-centric data-driven computation to underlying system resources. Based the above approach, we implemented a new parallel sweeps algorithm as a component in JAxMIN, featuring patch-angle parallelism, vertex clustering and hierarchical priority strategy. Evaluation with two real Sn software packages demonstrates that JSweep can scale to at least 49,152 cores for unstructured meshes and 76,800 cores for structured meshes with reasonable parallel efficiency.


What distinguishes JSweep from most counterparts is that we consider sweep computations in the context of general mesh-based application frameworks which would be critical in future coupling of multiple multi-physics simulations. Besides S_n transport sweeps, our abstraction also supports other data-driven algorithms well, e.g., particle trace which we have implemented as another component in JAxMIN. Given the increasing importance of data-driven computation and demands on coupling multi-physics simulations, we believe efficient support to both BSP and asynchronous data-driven models are necessary to construct high performance applications.








§ DISCUSSION



In practice, we find that for many applications of structured meshes, overhead of DAG operations is much larger than numerical computation itself. One natural solution to this problem is organizing the computation in the grain of cell blocks rather than the single cell, just like Koch-Baker-Alcouffe(KBA). However, the previous techniques of static cell blocking is not satisfied to our practical situations where the mesh is non-confirming or mosaic. Alternatively, a more elegant solution has been proposed in recent work  <cit.>, which is completely compatible with the approach given in this paper. Its essential idea is to formulate a coarsened graph by dynamically clustering computational vertices during a normal sweeping procedure, and later just do the parallel sweeps by scheduling the vertex clusters dictated by the formulated coarsened graph.


ICF <cit.><cit.> consists of multiple physics including laser ray tracing, dynamic reaction, radiation diffusion and radiation transport. Based on JASMIN, all the above physics simulations can be coupled and implemented in the LARED-I<cit.> code. Among these procedures, two components, i.e., particle trace in laser ray tracing and sweeps in radiation transport, are implemented in the data-driven approach.


Besides Sn transport solver, the downstream sweeping are also used for realization of many other numerical cores for convection dominated or Navier-Stokes equation on rectangular meshes and on unstructured meshes. All these realizations must be carefully designed according to the characteristics of zone shapes, discrete stencils, downstream directions and so on.


Technically, Charm++ plus another patch-based mesh application framework would implement the same approach described in this paper.








1
jasmin
Z. Mo, A. Zhang, X. Cao, Q. Liu, X. Xu, H. An, W. Pei, S. Zhu, JASMIN: a parallel software infrastructure for scientific computing, Frontiers Computer Science of China, 2010, 4(4):480–488. 

jaumin
Q. Liu, W. Zhao, J. Cheng, Z. Mo, A. Zhang and J. Liu, A programming framework for large scale numerical simulations based on unstructured meshes, in Proc. HPSC, NewYork, Apr. 6–8, 2016.

top500
http://top500.org/2015-nov, Nov., 2015.

boltzmann
R. L. Bowers, J. R. Wilson, Numerical modeling in applied physics and astrophysics, Jones and Bartlett publishers, 1991.













dd-app-1-5
J. Bey, G. Wittum, On the robust and efficient solution of convection diffusion problems on unstructured grids in two and three space dimensions, Applied Numerical Mathematics, 1997, 23(1):177–192.




dd-app-2-2
F. Wang, J. Xu, A cross-wind strip block iterative method for convection-dominated problems, SIAM Journal of Computing, 1999, 21:646–665.




downar
T. Downar, A. Siegel, C. Unal. Science Based Nuclear Energy Systems Enabled by Advanced Modeling and Simulation at the Extreme Scale. White Paper on Integrated Performance and Safety Codes, 2009.

t3d
R. Baker, R. Alcouffe. Parallel 3-D Sn Performance for MPI on Cray-T3D. In Proc. Joint International Conference on Mathematics Methods and Supercomputing for Nuclear Applications, New York, Oct., 1997.

cm200
R. Baker, K. Koch. An Sn algorithm for the massively parallel CM-200 computer. Nuclear Science and Engineering, 1998, 28: 312–320.

ardra
Lawrence Livermore National Laboratory, Ardra: Scalable parallel code system to perform neutron and radiation transport calculations, http://www.llnl.gov/casc/ardra.




kba-improvement-2
W. Hawkins, et al., Efficient Massively Parallel Transport Sweeps, Trans. Am. Nucl. Soc., 107, 477, 2012.

sweep3d
Los Alamos National Laboratory. The ASCI Sweep3d Benchmark. http://www. ccs3.lanl.gov/pal/software/sweep3d.

bsp
G. Valiant. A bridging model for parallel computation. Communications of the ACM, 1990, 33(8):108–111.

consensus
J. Misra. Detecting termination of distributed computations using markers. In Proc. PODC, pages 290–294, 1983.

plimpton-2000
S. Plimpton, B. Hendrickson, S. Burns, W. McLendon. Parallel algorithms for radiation transport on unstructured grids. In Proc. SC, Dallas, Nov., 2000.




mo
Z. Mo, A. Zhang and X. Cao. Towards a parallel framework of grid-based numerical algorithms on DAGs. In Proc. IPDPS, Greece, 2006.

dataflow
Hewitt C, Bishop P, Steiger R. A Universal Modular Actor Formalism for Artificial Intelligence. In Proc. IJCAI, San Francisco, 1973.




metis
G. Karypis, V. Kumar. Multi-level graph partitioning schemes. In Proc. of ICPP, Urbana-Champain, Aug., 1995, pp.113–122.

chaco
B. Hendrickson, R. Leland. A multilevel algorithm for partitioning graph. In Proc. of SC, San Diego, 1995.




pautz
S. Pautz. An Algorithm for Parallel Sn Sweeps on Unstructured Meshes. Nuclear Science and Engineering, 2002, 140(2): 111–136.




pmodel
M. Mathis, D. Kerbyson. A General Performance Model of Structured and Unstructured Mesh Particle Transport Computations. Journal of Supercomputing, 2005, 34:181–199.










parsec-sweep
S. Moustafa, M. Faverge, L. Plagne and P. Ramet. 3D Cartesian Transport Sweep for Massively Parallel Architectures with PaRSEC, In Proc. IPDPS, 2015.




frameworks-survey
A. Dubey, et al., A survey of high level frameworks in block-structured adaptive mesh refinement packages. J. Parallel Distrib. Comput. (2014), http://dx.doi.org/10.1016/j.jpdc.2014.07.001.

tort
W. A. Rhoades and D. Simpson, The TORT Three-Dimensional Discrete Ordinates Neutron/Photon Transport Code, ORNL/TM-13221, Oct., 1997.

jsnt-s
T.P. Cheng and L. Deng, JSNT-S manual, IAPCM, 2015.




3dsn
J.X. Wei,  JSNT-U (3DSn) manual, IAPCM, 2010.

sn-colomer
G. Colomer, R. Borrell, F. Trias, and I. Rodrieguez, Parallel algorithms for Sn transport sweeps on unstructured meshes. Journal of Computational Physics, vol. 232, no. 1, pp. 118–135, 2013.

samrai
SAMRAI. https://computation.llnl.gov/casc/SAMRAI, May 31, 2010

pumi
D. A. Ibanez, E. S. Seol, C. W. Smith and Mark S. Shephard, PUMI: Parallel Unstructured Mesh Infrastructure. ACM Transactions on Mathematical Software, 2015.




denovo
G. Davidson, T. Evans, J. Jarrell, S. Hamilton and T. Pandyam, Massively Parallel Three-dimensional Transport Solutions for the k-eigenvalue Problem, Nuclear Science and Engineering, 2014, 177:111-125.

denovo-old
T. Evans, A. Stafford, R. Slaybaugh and K. Clarno, Denovo: a new three-dimensional parallel discrete ordinates code in scale. Nuclear Technology, 2010, 171(8), 171-200.




mo2014 
Z. Mo, A. Zhang, and Z. Yang, A new parallel algorithm for vertex priorities of data flow acyclic digraphs, The Journal of Supercomputing, vol. 68, no. 1, pp. 49-64, 2014.

pautz-2015
S. Pautz and T. Bailey, Parallel Deterministic Transport Sweeps of Structured and Unstructured Meshes with Overloaded Mesh Decompositions, Proc. Joint International Conference on Mathematics and Computation (M&C), Supercomputing for Nuclear Applications (SNA) and the Monte Carlo (MC) Methods, Nashville, TN, April 19-23, 2015.








lared-all
P. Song, C.L. Zhai, S.G. Lee, et.al, LARED-I: The Integrated Code for Laser-driven Inertial Confinement Fusion. High Power Laser and Particle Beam (in Chinese), 2015, 27(3):54-60.

multiphysics-reactor
D. R. Gaston, C. J. Permann, J. W. Peterson, A. E. Slaughter, D. Andrsie, Y. Wang, M. P. Short, D. M. Perez, M. R. Tonks, J. Ortensi, Ling Zou and R. C. Martineau, Physics-based multi-scale coupling for full core nuclear reactor simulation, Annals of Nuclear Energy, 2015, 84:45–54.





parsec
G. Bosilca, A. Bouteiller, A. Danalis, T. Herault, P. Lemarinier and J. Dongarra, DAGuE: A generic distributed DAG engine for High Performance Computing, Parallel Computing, vol.38, no.1-2, 2012.

charmpp
L.V. Kale, E. Bohm, C. L. Mendes, T. Wilmarth and G. Zheng, Programming peta-scale applications with Charm++ and AMPI, Peta-scale Computation: Algorithms Appl. 1(2007): 421–441.


pregel
G. Malewicz, M. Austern, A. Bik, J. Dehnert, I. Horn, N. Leiser and G. Czajkowski. Pregel: a system for large-scale graph processing. In SIGMOD, 2010.

powergraph
J. Gonzalez, Y. Low, H. Gu, D. Bickson, C. Guestrin. PowerGraph: distributed graph-parallel computation on natural graphs. In OSDI, 2012.

graphine
J. Yan, G. Tan, Z. Mo, N. Sun, Graphine: programming graph-parallel computation of large natural graphs for multicore clusters, In IEEE Transaction on Parallel and Distributed Systems, 27(6):1647-1659, 2016.








[Sn]
Radiation effect are often modeled by the discrete ordinates (S_n) form of Boltzmann transport equation


    Ω_m ·∇ I_mg = R(I, t)
       =S_mg - (σ_A+σ_S)I_mg +∑_g'∑_m' w_m σ_S, g' → g, m' → m I_m'g',


    Ω_m ·   ∇ I_mg^i+1 = R(I^i+1, t^i+1),

where the radiative flux I_mg(x,y,z) is discretized into a set of energy groups g and angular directions (ordinates) m  <cit.>. Equation (<ref>) describes how a single flux I_mg propagates in a particular direction Ω_m. At each time step, the solution for a grid space is got by solving a large set of coupled equations describing flux propagations for all m and g, which iteratively continues until the convergent condition is met. The widely used method for this problem is source iteration, e.g., to treat σ_A and σ_S as constants (in each iteration) so that any individual I_mg can be computed independently. In the same iteration, all I_mg can be computed simultaneously with collisions resolving in runtime. For every I_mg, each iteration can be divided into two phases, that (1) to compute local scattering source, and (2) to compute sweep flux from source across the grid in the downstream direction. We simply denote the second phase as Sweeps, which is the most time-consuming portion in source
iteration method.

Radiation effects are often modeled by the discrete ordinates (S_n) form of Boltzmann transport equation. The standard solving method is source iteration, in which the solution is computed by iteratively repeating two phases. First, compute local scattering source. Second, for each ordinate direction, sweep the radiation flux from source across the grid in the downstream direction. In the second phase, sweeps from all ordinate directions usually can be carried out in parallel. The second phase is also the most time-consuming portion, which we simply denote as Sweep(s). Alg. <ref> is the source iterative algorithm of S_n method on solving BTE. The steps 1, 2.1 and 2.3 are typical BSP procedures that compute on all or part of cells, and can be efficiently implemented by components of JAxMIN. For example, the steps 1 and 2.1 can be implemented by instantiating normal numeric components, while the step 2.3  can be implemented in the reduction component. However, the step 2.2, i.e., sweep, is essentially a data-driven procedure and BSP-inefficient.




language=c++, basicstyle=, tabsize=2, frame=lines, numbers=left, numberstyle=, numbersep=1pt, numberblanklines=false, belowcaptionskip=4pt, xleftmargin=4pt, mathescape=true
emph=compute, sweep, reduce,, emphstyle=
[float, caption=Source Iterative Sn Solver in JAxMIN Components,label=jaxmin-source-iteration, firstnumber=1]
void Integrator::initializeLevelData(...)

	//set data dependencies between cells/angles.
	d_sweep_intc->setDataDependency(level, ...);

void Integrator::advanceLevel(Level    level,...)

	double err;
	do 
		// update scattering source
		d_scatter_intc->compute(level, ...);
		// update flux
		d_sweep_intc->sweep(level, ...);
		// compute the iterative error
		d_err_intc->reduce(   err, 1, level, ...);
	 while(err > err_tolerance);



 <ref> is the complete workflow of a data-driven S_n Sweeps for static grids, including three main procedures. First, it eliminates cycles in the unstructured grid, and then transforms the grid into a computable DAG^1. Plimpton et.al. <cit.> proposed an algorithm to detect and eliminate the cycles. Second, it partitions the DAG into multiple subgraphs that are then distributed to processes (processors). The graph partitioning has a direct impact on load balance and further parallel efficiency. Fortunately, previous works <cit.><cit.><cit.>  on graph theory and parallel computing have already provided solutions to this problem. Third, it iteratively carries out Sweeps computation for required time steps. Each processor performs sweep computation on its local DAG, and communicates with other processors if there are cutting edges between their subgraphs. As a matter of fact, while the former two procedures can be considered as pre-processing steps, the third one — Sweeps, is the main computation part in the S_n Sweeps algorithm and thus is our redesigning target in this paper.






§ COVER LETTER




 §.§ Revision summary

In this submission, we have revised the original paper thoroughly. In particular, significant changes are as follows.



  * Added the optimization of coarsened graph in Sec.<ref>. It is used in the new experiments of Sec.<ref> and <ref>.

  * Added detailed profiling and analysis of JSweep's runtime overhead in Sec.<ref>.

  * Added performance comparison of JSweep against previous JASMIN and JAUMIN in Sec.<ref>.

  * Rewritten paragraphs about progress tracking and program termination in Sec.<ref>  and <ref>, making it more clear.




 §.§ Explanations to synthesized reviews




  §.§.§ Comparison with other systems



Synthesis review #1: A comparison of JSweep performance with that achieved by one or more state-of-the-art systems, for structured/unstructured mesh sweeps, with clear details on the benchmark computations.


To address the reviewers' concern, we have tried our best to compare the performance of JSweep to other available systems, by providing extra experimental results and comparing with data in literatures. These new contents are added as Sec.<ref>.

The first approach we have taken is to compare JSweep with mainline JAxMIN implementations. Since JAxMIN supports sweep for a long time, extensive algorithmic and architecture-oriented performance optimizations are already incorporated (<cit.>). It shall be acceptable to consider it one of the state-of-the-art systems. Another motivation of this choice is that we can implement exactly the same S_n benchmark for a fair comparison. The experiments show that JSweep outperforms the mainline JAxMIN constantly.

The second approach we have taken is to compare the performance we observed on the Kobayashi problem to the data in the latest literature. The Kobayashi problem is a widely used benchmark in reactor neutron transport to assess performance of S_n sweep performance. The data shows that JSweep obtain comparable results with specially optimized implementations in Denovo.



Due to the unavailability of source code and lack of functionality, we did not compare JSweep directly with systems used in <cit.> and <cit.>.





  §.§.§ Clarify Progress Monitoring and Progress Guarantee



Synthesis review #2: A more detailed treatment on monitoring of progress, especially for unstructured meshes: how is progress guaranteed with batched notifications; how does order of local vertex selection affect it


To address the reviewers' concern, we have modified Sec.<ref>  and <ref> to make it more clear, with respect to distributed progress tracking.

Progress is guaranteed for both unstructured and structured meshes. Firstly, the batched notification is not strict. There are no such restrictions as “group these five vertexes and process their output in one stream”. Instead, the runtime system just puts per-vertex streams in a batch queue, which is handled by the master thread in Fig.<ref>. The master thread is always spinning processing available messages, so messages are always delivered, and no dead lock is possible given the graph is acyclic. Secondly, priority strategies of both patch-programs and (cell, angle)s within patch-program never violate the topological order in the original DAG. Thirdly, with respect to distributed implementation, we have handled the inconsistent state when there are no active patch-programs globally but on-the-fly communication that will activates patch-programs by forcing the process to ensure all data streams are finished before negotiating for termination.

Priority strategies and vertex clustering do affect performance, though. We evaluate their effect in Section <ref> and <ref>.

In this paper, as other counterpart such as <cit.> and <cit.>, we exclude the case of cyclic directed cyclic graphs. These cases can exist in real-world applications, especially for deforming structured grids and unstructured grids, but they are often preprocessed by other algorithms to eliminate the cycles. Our used meshes, especially the unstructured reactor and ball, are carefully handled to avoid potential sweeping loops in phases of initial mesh construction and refinement. Along with developers of JSNT-U, we have developed some efficient algorithms for special meshes to do loop detection and eliminations.



  §.§.§ Analysis of JSweep Overhead



Synthesis review #3: In-depth analysis of overheads, with details on compute-cost per vertex, and amount of data communicated per edge.


To address the reviewers' concern, we have carried out a detailed overhead analysis of JSweep under real application scenario and added Sec. <ref>.

A small scale Kobayashi problem with a 200x200x200 mesh is used with the JSNT-S application to study the JSweep overhead. The Kobayashi problem is a standard approach to assess sweep performance, so the choice is reasonable. To gain more insight into the evolution of overhead, we carried out a strong scaling experiments where the portion of computation decrease w.r.t. the growth of CPU cores, so we can study the overhead when there are little computations per CPU core. We instrumented JSweep with RDTSC timer to obtain accurate timing data with the lowest overhead.

We divide the runtime of a sweep into graph-op (DAG operation introduced by DAG scheduling), pack/unpack (stream operation), kernel (actural computation on the vertex), comm (communication of streams), and idle (spinning without meaningful work). The result in Sec. <ref> shows that graph-op and pack/unpack is under control even when scaling to the limit, and the first source of inefficient is always the “idling”. This shall justify the adoption of the data-driven approach since it is general enough to handle complex cases yet introduce only moderate overhead, which itself can be further optimized.


