Mathilde Boltenhagen
email: boltenhagen@lsv.fr

Thomas Chatain
email: chatain@lsv.fr

Josep Carmona
email: jcarmona@cs.upc.edu

Encoding Conformance Checking Artefacts in SAT

come   L'archive ouverte pluridisciplinaire

Introduction

On its core, conformance checking relies on the computation of artefacts that link observed and modeled behavior, which are used with different purposes: spotting deviations, evaluating quality metrics of a process model, extending a process model with evidence-based information, among others [START_REF] Carmona | Conformance Checking -Relating Processes and Models[END_REF].

Conformance checking is expected to be the fastest growing segment within Process Mining in the years to come 3 . Still, the field is facing several challenges. Among them, we highlight two important ones: techniques for a sound replay of event data on top of process models, and the advent of faithful metrics for evaluating process models with respect to observed behavior.

The former challenge is strongly related to the notion of alignment artefact: given a trace and a process model, find a run in the process model that is as close as possible (in edit distance terms) to the observed trace. The seminal work in [START_REF] Adriansyah | Aligning observed and modeled behavior[END_REF] describes an algorithm for computing alignments based on A * . Alternatives to this algorithm have appeared recently: in [START_REF] De | Aligning real process executions and prescriptive process models through automated planning[END_REF], the alignment problem is mapped as an automated planning instance. Automata-based techniques were proposed in [START_REF] Reißner | Scalable conformance checking of business processes[END_REF][START_REF] Sander | Scalable process discovery and conformance checking[END_REF]. Finally, the work in [START_REF] Vincent Bloemen | Symbolically aligning observed and modelled behaviour[END_REF] proposes using binary decision diagrams to alleviate the computation of alignments. In this paper, we provide an alternative to the aforementioned techniques, that is based on encoding the computation of an alignment as a SAT instance. We also show how to encode also in SAT multi-alignments [START_REF] Chatain | Alignment-based trace clustering[END_REF], that generalize the notion of alignments so that not one but several traces are considered when computing the closest process model run.

The latter challenge is mainly concerned with the proposal of sound and meaningful metrics for precision and generalization, nowadays acknowledged as the quality dimensions with less convincing estimations (e.g., [START_REF] Tax | The imprecisions of precision measures in process mining[END_REF]). Anti-alignments, presented in [START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF], are an effective conformance artefact to foresee those model runs that deviate most (again, in terms of edit distance) with respect to a trace or a complete log. In [START_REF] Boudewijn | Aligning modeled and observed behavior: A compromise between computation complexity and quality[END_REF] it is shown how anti-alignments can be used to provide a more consistent estimation to both precision and generalization. In this paper we provide for the first time an encoding in SAT of anti-alignments that relies on the edit distance between traces, in contrast to the Hamming distance used in [START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF][START_REF] Boudewijn | Aligning modeled and observed behavior: A compromise between computation complexity and quality[END_REF]; to the best of our knowledge this is the first implementation of anti-alignments over the edit distance. The use of Hamming distance for anti-alignment shows important limit highlighted by Fig. 1. In summary, this paper presents SAT as a common means to compute important conformance checking artefacts. The prototype implementation, although not mature and therefore without having any optimization, shows interesting features that can make it a good alternative for the state of the art techniques for the same task in the near future.

The paper is organized as follows: next section provides related techniques for the tasks considered in this paper. Then in Section 2 we provide the background for understanding the paper. Section 3 shows how to encode the computation of a run at a given edit distance of a trace. Then Section 4 shows how to adapt this encoding to particular conformance checking artefacts. Section 5 reports experiments, whilst Section 6 concludes the paper.

Related work. Levenshtein's edit distance is commonly used in Process Mining to define similarities between traces [START_REF] Jagadeesh | Abstractions in process mining: A taxonomy of patterns[END_REF][START_REF] Jagadeesh | Context aware trace clustering: Towards improving process mining results[END_REF][START_REF] Tax | Event abstraction for process mining using supervised learning techniques[END_REF], and to align log traces to model traces [START_REF] Jagadeesh | Trace alignment in process mining: opportunities for process diagnostics[END_REF][START_REF] Adriansyah | Aligning observed and modeled behavior[END_REF][START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF]. Works of the literature on alignments use the definition of edit distance more or less Finally, the d loop is enabled when a user donates to the developer of the app.

implicitly [START_REF] Jagadeesh | Trace alignment in process mining: opportunities for process diagnostics[END_REF][START_REF] De | Data-aware process mining: discovering decisions in processes using alignments[END_REF][START_REF] Adriansyah | Alignment based precision checking[END_REF]. Alternatively, distance between process models have been characterized in several works [START_REF] Bae | Development of distance measures for process mining, discovery and integration[END_REF][START_REF] Dijkman | Aligning business process models[END_REF]. The latter uses a dependency edit distance which is computed in a similar way as our approach.

Recent studies focus on SAT implementation of Data Mining algorithm in order to satisfy all the constrains and get optima [START_REF] Métivier | Constrained clustering using sat[END_REF][START_REF] Davidson | A sat-based framework for efficient constrained clustering[END_REF]. By introducing a SAT implementation of alignments in this paper, we hope to push a new family of algorithmic methods for conformance checking in the line of [START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF][START_REF] Boltenhagen | Generalized alignment-based trace clustering of process behavior[END_REF]. However, these works mostly consider Hamming distance between log traces and process models, which is usually considered less appropriate than edit distance (c.f. Fig. 1).

The SAT encoding of the edit distance between words has already been introduced in previous works [START_REF] Backurs | Edit distance cannot be computed in strongly subquadratic time (unless seth is false)[END_REF][START_REF] Groce | Error explanation with distance metrics[END_REF]. The work in [START_REF] Backurs | Edit distance cannot be computed in strongly subquadratic time (unless seth is false)[END_REF] studies the computation of edit distance for its interest in complexity theory.

Preliminaries

We use labeled Petri nets as process models.

Definition 1 (Process Model (Labeled Petri Net) [START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF]). A Process Model defined by a labeled Petri net system (or simply Petri net) is a tuple N = P, T, F, m 0 , m f , Σ, Λ , where P is the set of places, T is the set of transitions (with P ∩ T = ∅), F ⊆ (P × T ) ∪ (T × P ) is the flow relation, m 0 is the initial marking, m f is the final marking, Σ is an alphabet of actions and Λ : T → Σ ∪ {τ } labels every transition by an action or as silent.

A marking is the set of places that contain tokens for a given instant. A transition x can fire if all the places before x, noted • x def = {y ∈ P | (y, x) ∈ F }, are marked. When a transition fires, all the tokens in • x are removed and all the places in x

• def = {y ∈ P | (y, x) ∈ F } become marked. A marking m is reachable from m if there is a sequence of firings t 1 . . . t n that transforms m into m , denoted by m[t 1 . . . t n m .

Definition 2 (Full runs).

A full run of a model N is a firing sequence t 1 . . . t n of transitions that can transform the initial marking m 0 of N to the final marking m f of N . We note Runs(N ) the full runs of N .

The run s, f, b, a is a full run of the model of Fig. 2. Now we formalize logs: Definition 3 (Log). A log L over an alphabet Σ is a finite set of words σ ∈ Σ * , called log traces.

Alignments

Aligning log traces to model traces has been identified as a central problem in Process Mining [START_REF] Adriansyah | Aligning observed and modeled behavior[END_REF]. The problem is to find a run in the process model that is as close as possible to the observed trace. This closeness is usually defined in terms of the edit distance: Definition 4 (Edit distance). The edit distance dist(u, v) between two words u and v ∈ Σ * is the minimal number of edits needed to transform u to v. In our case, edits can be deletions or additions of a letter in words.

Example 1. Considering words u = s, g, c and v = s, b, c, a the number of edits to transform u to v is indeed 3. The letter g has to be removed and the letters b and a inserted. Then the two words are at distance 3.

Formally, we define alignments in a way that not only achieves the optimal edit distance between the log trace and a model trace, but also makes explicit where they match and mismatch. Definition 5 (Alignment, optimal alignment). An alignment of a log trace s = s 1 , . . . , s m ∈ L to a run u = u 1 , . . . , u n of process model N = P, T, F, m 0 , m f , Σ, Λ is a sequence of moves (s 1 , u 1 ), . . . , (s p , u p ) with p ≤ m + n such that each move (s i , u i ) is either:

• (e i , t i ) with e i = Λ(t i ) for a synchronous move • (e i , ) for a log move ( is a skip symbol which indicates the mismatch), i.e. e i is deleted in the trace • ( , t i ) for a model move, i.e. Λ(t i ) is inserted in the trace; projection of s 1 , . . . , s p to Σ * (which drops the occurrences of ), yields s projection of u 1 , . . . , u p to T * (which drops the occurrences of ), yields u An alignment between u ∈ Runs(N ) and s is optimal if it minimizes the number of occurrences of . This minimal number of mismatches corresponds to the edit distance dist(u, s) between u and s.

Alignments can be represented in a two-row matrix. Next figure shows an alignment between the second log trace ( s, g, c ) and the run s, b, c, a of the model of Fig. 2. Aligning traces to model often means searching the minimal number of moves, i.e. the optimal alignment. Different cost functions are used but the most common one applies the weights 0 for a synchronous move and 1 for a log move or a model move.

Next figure shows an alignment for the trace s, g, c and the run s, g, c, τ of the model of Fig. 2. Since τ moves inccur into no cost, the alignment has cost 0, and hence it is optimal. trace s g c run s g c τ

SAT Encoding of the Edit Distance

We first introduce our SAT encoding of the edit distance between two words, which will serve as a building block for alignment, multi-alignment and anti-alignment.

The Boolean satisfiability (or SAT) problem, is the problem of determining, for a given Boolean formula, if there exists a combination of assignments to the variables that satisfies it. In the case of alignment, a SAT formula would encode the following question: Does it exist an alignment between the trace σ and the model N for a cost d ? In other words, we are looking for a formulas that encodes the edit distance d between a trace and a run of the model. Our encoding is based on the same relations that are used by the classical dynamic programming recursive algorithm for computing the edit distance between two words u = u 1 , . . . , u n and v = v 1 , . . . , v m :

               dist( u 1 , . . . , u i , ) = i dist( , v 1 , . . . , v j ) = j dist( u 1 , . . . , u i+1 , v 1 , . . . , v j+1 ) =    dist( u 1 , . . . , u i , v 1 , . . . , v j ) if u i+1 = v j+1 1 + min(dist( u 1 , . . . , u i+1 , v 1 , . . . , v j ), dist( u 1 , . . . , u i , v 1 , . . . , v j+1 )) if u i+1 = v j+1
We encode this computation in a SAT formula φ over variables δ i,j,d , for i = 0, . . . , n, j = 0, . . . , m and d = 0, . . . , n + m. The formula φ will have exactly one solution, in which each variable

δ i,j,d is true iff dist( u 1 . . . u i , v 1 . . . v j ) ≥ d.
In order to test equality between the u i and v j , we use variables λ i,a and λ j,a , for i = 0, . . . , n, j = 0, . . . , m and a ∈ Σ, and we set their value such that λ i,a is true iff u i = a, and λ j,a is true iff v j = a. Hence, the test u i+1 = v j+1 becomes in our formulas: a∈Σ (λ i+1,a ∧ λ j+1,a ). For readability of the formulas, we refer to this coding by [u i+1 = v j+1 ]. We also write similarly

[u i+1 = v j+1 ].
In the following, we describe the different clauses of the formula φ of our SAT encoding of the edit distance.

δ 0,0,0 ∧ d>0 ¬δ 0,0,d (1)
d n i=0 (δ i+1,0,d+1 ⇔ δ i,0,d ) (2)
d n j=0 (δ 0,j+1,d+1 ⇔ δ 0,j,d ) (3)
d n i=0 n j=0 [u i+1 = v j+1 ] ⇒ (δ i+1,j+1,d ⇔ δ i,j,d ) (4)
d n i=0 n j=0 [u i+1 = v j+1 ] ⇒ (δ i+1,j+1,d+1 ⇔ (δ i+1,j,d ∧ δ i,j+1,d )) (5)
Example 2. At instants i = 1 and j = 1 of words u = s, g, c and v = s, b, c, a , the letters are the same, then, by ( 4), the distance is only higher or equal to 0 : (u 1 = v 1 ) ⇒ (δ 1,1,0 ⇔ δ 0,0,0 ). However at instants i = 2 and j = 2, the letters u 2 and v 2 are different. A step before, δ 1,2,1 and δ 2,1,1 are true because of the length of the subwords. Then, by ( 5), the distance at instants i = 2 and j = 2 is higher or equal to 2 : δ 2,2,2 . The result is understandable because the edit distance costs the deletion of g and the addition of b to transform u to v.

SAT Encoding of Conformance Checking Artefacts

The distance between log traces and a process model is not only a distance between words, since process model describe a (possibly infinite) language. In this part, we recall the SAT encoding of full runs of Petri nets [START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF] and combine the implementation with the edit distance.

SAT implementation of process models. For a Petri net N = P, T, F, m 0 , m f , Σ, Λ and n the size of the full runs, the variables m i,p , with i ∈ {0..n} and p ∈ P , represent the marking at instant i. The variables τ i,a encode a firing transition t ∈ T labelled by a ∈ Σ at instant i ∈ {0..n}4 . The following constraints encode the semantics of the Petri net.

-Initial marking:

( p∈m0 m 0,p ) ∧ ( p∈P \m0 ¬m 0,p ) (6)
-Final marking:

( p∈m f M n,p ) ∧ ( p∈P \m f ¬m n,p ) (7)
-One and only one t i for each i:

n i=1 a∈Σ (τ i,a ∧ a ∈Σ\t ¬τ i,a ) (8)
-The transitions are enabled when they fire:

n i=1 a∈Σ (τ i,a =⇒ p∈ • t m i-1,p ) (9)
-Token game (for safe Petri nets):

n i=1 a∈Σ p∈t • ,Λ(t)=a (τ i,a =⇒ m i,p ) (10)
n i=1 a∈Σ p∈ • t\t • ,Λ(t)=a (τ i,a =⇒ ¬m i,p ) (11)
n i=1 a∈Σ p∈P,p ∈ • t,p ∈t • ,Λ(t)=a (τ i,a =⇒ (m i,p ⇐⇒ m i-1,p )) (12)
Process models are now implemented in a CNF formula and can be combined with log traces for alignments.

SAT Edit Distance for Alignments

Log traces are sequences of activities that can be considered as words and implemented as presented in Section 3. SAT encoding of process models has been recalled in previous section. All the above clauses are considered in the SAT implementation of alignments. A last series of constraints is needed to be appended, to relate the fired transitions, represented by the τ i,a , with the actions in the corresponding model trace, represented by variables λ i,a from the encoding of Section 3:

n i=1 a∈Σ (λ i,a ⇐⇒ τ i,a ) (13)
Example 3. All the full runs of the process model of Fig. 2 contain a s at the first instant. So the variable τ 1,s is true. If the log trace is σ = s, f, g then, λ 1,s is true which implies δ 1,1,0 by (4).

Minimization of the edit distance. The conjunction of the previous clauses for the full runs of the model and the their edit distance to a given log trace σ, gives a formula which has one solution per full run of the model. With each solution, the values of the δ n,m,d determine the edit distance between the corresponding model trace and σ. Our goal for optimal alignments is to minimize this distance, which corresponds to the number of variables assigned to true among the δ n,m,d . Pseudo-Boolean solvers like MINISAT+ deal with minimization objectives under the form of a weighted sum of variables; in our case:

d 1 × δ n,m,d .
How to deal with runs of different length. In order to consider different sizes of traces and different sizes of runs, we added a loop on a wait activity on the final marking of the model. The SAT encoding of the edit distance is adjusted so that skipping a wait activity does not increment the distance between words.

Example 4. Fig. 3 shows optimal alignments of every trace of the log of Fig. 2. The deviating trace s, a, a is then highlighted by the distance to its alignment.

SAT Implementation for Multi-alignments

Multi-alignments were introduced in [START_REF] Chatain | Alignment-based trace clustering[END_REF] as a generalization of alignments. Multialignments were used to define a model-based trace clustering method. Instead of aligning a trace to a run of a process model, they align a set of log traces (typically from the same cluster) to a common run of the model.

Definition 6 (Multi-alignment). Given a finite collection C of log traces and a model N , an (optimal) multi-alignment of C to N is a full run u ∈ Runs(N ) which minimizes the sum σ∈C dist(σ, u).

The SAT implementation of multi-alignment requires us to duplicate the variables λ σ i,a that represent actions in the log traces σ ∈ L and the variables δ σ i,j,d that measure the edit distance to the model trace. Similarly to Section 4.1, the optimal multialignment is found by minimizing the number of variables assigned to true in the following objective:

d σ∈L 1 × δ σ n,|σ|,d .
Example 5. We computed the multi-alignment of the model and the full log of Fig. 2. The optimal multi-alignment is the full run s, f, b, a which is at distance d ≤ 3 to all the log traces.

SAT Implementation for Anti-alignments

Anti-alignment was introduced in [START_REF] Chatain | Anti-alignments in conformance checking-the dark side of process models[END_REF]. Contrary to multi-alignments, the aim of antialignments is to get, for a given log, the run of a model which differs as much as possible to all the traces in the log. The notion of anti-alignment is used in some quality metrics like precision and generalization [START_REF] Boudewijn | A unified approach for measuring precision and generalization based on anti-alignments[END_REF].

Definition 7 (Anti-alignment). Given a finite collection L of log traces and a model N , an anti-alignment is a run u ∈ Runs(N ) which maximizes its distance σ∈L dist(σ, u) to the log.

The encoding is then very similar to the multi-alignment version. Instead of minimizing the distance to the set of log traces, we maximize it using the opposite minimization objective:

d σ∈L -1 × δ σ n,|σ|,d 5 .
Example 6. We computed the anti-alignment of model and the set of log traces of Fig. 2 . Limited by a maximum size of run to 8, the optimal anti-alignment found by our tool is s, b, d, f, d, d, d, τ . The minimal distance between each trace and this full run is 9. We compared our result with the module Anti-Alignment of ProM6 , that computes antialignment over Hamming distance. For the same size of run, the algorithm returned the sequence s, b, d, d, d, c, a, d that is indeed linearly far from the log traces. However as either letters c or a are present in every trace, the run looks more similar to the log than the one found with edit distance.

Experiments

The construction of SAT formulas for alignments, multi-alignments and anti-alignments is implemented in Ocaml in our tool DARKSIDER available on github 7 . The software invokes a SAT solver, by default MINISAT+. Examples of the previous section have been fully computed by the SAT formulas to get optimal solutions. Since the approach is heavily influenced by the size of the formulas constructed (which is large even for small models), in this section we focus on heuristics to simplify the formulas, at the expense of sacrifying optimality eventually. Those are preliminary results oriented towards illustrating how to instruct the current tool.

A first heuristic is to remove δ i,j,d variables, for d large . Limiting the maximal number of edits may not change the result when the traces are close to the model. A second approximation is to get optimal prefix alignment by limiting the size of the run.

Tab. 1 (a)-(c) summarizes the experiments. Notice that while for alignments we show average numbers over one hundred traces, for multi-and anti-alingments (where only a run is computed for the whole log) total numbers are provided. The first two columns (Model and |L|) describe the model and the log size, respectively. Column "Size of run" shows the maximal size allowed for the run in the model (which will be an alignment in (a), a multi-alignment in (b), and an anti-alignment in (c)). Sometimes the length is limited when PRE is specified, as explained above. Then the fourth column reports the maximal number of edits allowed, sometimes with a bound as explained above. When LIM is indicated, the distance between the model and the trace is larger than what was tested. The last three columns show the time to construct the formula and the total execution time for our approach, and the time needed in ProM. Notice that since we are providing results for edit-distance based conformance artefacts, only ProM results for alignments are shown.

The goal of this paper was to demonstrate the interest of our approach based on a SAT encoding of Petri net executions and edit distance between traces. For alignments, indeed, the current SAT encoding shows bad execution times compared to the optimized algorithms implemented in ProM. But we obtain, for the first time, an implementation of the anti-alignment and multi-alignment artefacts defined over the edit distance. We believe that there is a lot of space for improvement, for instance by optimizing the encoding or by using heuristics to very efficiently find approximations of the results. This paper has shown a unified approach to compute important conformance artefacts over SAT. Thanks to its high versatility, the encoding as SAT formula allows one to compute exact solutions to various problems (here alignments, anti-alignments and multialignments) under various optimality criteria. In particular, we show for the first time how anti-alignments can be computed for the edit distance by formulating the problem as a SAT instance. Although technically sound, the encodings proposed in this paper suffer from the explosion of the SAT formulas created. As main research direction, we are working into finding better encodings that can alleviate significantly the size of such formulae, but also incorporating other optimizations and heuristics (possibly inspired from optimization techniques used in automated planning) that can make the approach more efficient in practice.

Fig. 1 :

 1 Fig. 1: How Hamming distance penalizes alignment and anti-alignment : Trace b, a, b, a, b, a is closer to a, a, a, a, a, a than a, b, a, b, a, b according to Hamming distance. However there is only a shift of letter a and b between b, a, b, a, b, a and a, b, a, b, a, b . The model run a, b, a, b, a, b seems to be a better anti-alignment for trace b, a, b, a, b, a , as the edit distance shows.

Fig. 2 :

 2 Fig. 2: Example of process model of the behaviors of users rating an app. s represents the start activity, f and c indicate if the user sent a file or wrote a comment. Transitions g and b separate good and bad marks. Bad ratings get apologies noted by activity a.Finally, the d loop is enabled when a user donates to the developer of the app.

3 Fig. 3 :

 33 Fig. 3: Alignment of each trace and the model of Fig. 2

Table 1 :

 1 Experimental results for the computation of optimum and approximations of alignments and anti-alignments with our tool DARKSIDER, obtained on a virtual machine with CPU Intel Xeon 2.67GHz and 50GB RAM.

	Model			|L|	Size of run	Maximal number of	Formula construction	Total execution	ProM execution
	Reference |T | |P |			edits	time (sec)	time (sec)	time (sec)
	Fig. 2	8	7 100	7	5	0.239	0.349	0.002
	M8 of [25] 15 17 100	PRE: 20	LIM:10	10.139	15.530	0.001
	M1 of [25] 40 39 100	PRE: 7	LIM:10	4.924	7.16	0.005
	Loan [10] 15 16 100	PRE: 19	LIM: 10	14.047	20.915	0.002
					(a) Alignments (showing averages).
	Fig. 2	8	7	10 100	8 8	7 7	10.101 99.602	15.362 200.569
	M8 of [25] 15 17	10 100	18 PRE:15	LIM:6 LIM:6	252.471 516.391	414.174 741.162
	M1 of [25] 40 39	10 100	PRE: 13 PRE: 13	LIM:10 LIM: 5	115.706 681.95	172.500 1066.94
	Loan [10] 15 16	10 100	PRE: 19 PRE: 9	15 LIM:10	252.572 359.982	373.683 508.542
					(b) Multi-alignments.
	Fig. 2	8	7	10 100	8 8	LIM: 10 LIM: 10	13.802 137.213	21.502 243.842
	M8 of [25] 15 17	10 100	18 PRE: 10	LIM:10 LIM: 10	103.812 343.529	148.271 496.733
	M1 of [25] 40 39	10 100	39 PRE:13	LIM:10 LIM:5	1337.806 680.556	2069.505 995.361
	Loan [10] 15 16	10 100	PRE: 19 PRE:19	LIM: 10 LIM: 10	140.840 1526.048	203.257 2185.785
					(c) Anti-alignments.

https://www.marketsandmarkets.com/Market-Reports/process-analytics-market-254139591.html

For the sake of simplicity of the encoding, we are abusing a bit the notation, i.e., assuming that labels identify transitions. This can be generalized easily for the general case when severals transitions exist for the same label.

Anti-alignment Precision/Generalization of package AntiAlignments of ProM software version 6.8, http://www.promtools.org/

[START_REF] Jagadeesh | Trace alignment in process mining: opportunities for process diagnostics[END_REF] https://github.com/BoltMaud/darksider

Acknowledgments. This work has been supported by Farman institute at ENS Paris-Saclay and by MINECO and FEDER funds under grant TIN2017-86727-C2-1-R.


Piotr Łuczyński
email: luczynski@ikdg.rwth-aachen.de

David Erdmann

Dennis Többen

Mathias Diefenthal

Manfred Wirsum

Klaus Helbig

Wolfgang Mohr

Stefan Reitschmidt

Fast Calculation Methods for the Modelling of Transient Temperature Fields in a Steam Turbine in Pre-warming Operation

Keywords:

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The situation on world energy markets is changing dynamically. The rising awareness of the environment and especially of the reality and severity of climate change has marked the beginning of a real breakthrough in modern concepts of energy production. Instead of striving for the cheapest possible solutions, societies value environmentalfriendly, energy-sufficient and sustainable energy generation. Electricity production in big, centralized units is reduced, and increased, locally with wind, solar power and other renewables. In the European energy market, the share of renewable energies is increasing considerably over the last years. Regarding Germany, the politics support a priority feed-in electricity tariffs to encourage the use and development of new energy technologies. Similar policy mechanisms may be encountered in other European countries like France, Denmark or Poland. For almost all renewable energy resources the key issue is variability and a certain unpredictability. Therefore, wind parks or solar power always need back-up energy sources. Although 24-hour weather forecasts are reasonably accurate, wind speeds or radiation intensity can change quickly and unpredictably so that the back-up capacity should have the ability to react as quickly as possible. Fluctuation only in wind power may have a strong influence on balancing electricity supply with demand. An example is the sudden reduction of electricity generation from wind mills in Texas on February 26, 2008 [START_REF] Klimstra | Smart Power Generation[END_REF]. Due to rapid decrease of wind speed the generated wind power dropped from 1.7 GW to only 700 MW. In combination with evening peak in power demand of about 3 GW, the grid operator was forced to use a rapid demand reduction of 1.1 GW for sake of system stability [START_REF] Klimstra | Smart Power Generation[END_REF]. This has been reflected in electricity prices peaks and it offered additional profit to the electricity producers, who were able to cover the unbalanced demand. In the light of the above considerations, the current situation changes the role of the conventional power plants on the energy market and produces new, challenging requirements for modern steam turbines regarding the flexibility and increasing number of startups due to the high share of fluctuating power input of renewable generation. Especially hard coal fired power plants already have to shut down on the weekends, for periods of up to 60 hours [START_REF] Wirth | Recent Facts about Photovoltaics in Germany[END_REF]. As a consequence a cold start-up is required afterwards. The frequent start-up and shut-down processes increase the lifetime consumption of the turbine components as they are not designed for this operational pattern. Furthermore, the conventional plants have to start-up as fast as possible in order to guarantee the security of supply and grid stability. Fast start-up causes high metal temperature gradients. Hence, with decreasing start-up time the thermal stress, particularly of the thick-walled components, increases. This leads to an overly high lifetime consumption. Especially critical are the HP and IP steam turbines due to the high temperature difference between inlet steam and rotor temperature after a long period of idleness. Thermal stress is proportional to the difference between rotor surface and rotor average temperature. Faster loading means higher stress, which is limited not to overuse lifetime. Cold starts have the highest temperature difference from start to base load, thus, a start-up usually requires five to eight hours, which does not longer correspond to recent market requirements [START_REF] Vogt | Flexibility improvement of the steam turbine of Conventional or CCPP[END_REF]. For the sake of shorter start-up procedure, the crucial components of steam turbine have to either be kept warm during the idle times or efficiently pre-warmed before operation. Considering a conventional power plant it is commonly known, that the thick-walled components staying in contact to the hot steam limit the start-up speed [START_REF] Vogt | Flexibility improvement of the steam turbine of Conventional or CCPP[END_REF]. On the side of the boiler the temperature gradients in the heat exchanger headers, valves and drums have to be minimalized during start-up procedure. On the side of the steam turbine the special treatment requires the rotors and split lines. The pre-warming or warm-keeping procedure reduces the temperature difference between the steel and the steam and hence, the evoked stresses within the component decrease. Recent studies have shown different possibilities to heat-up the rotor and the casing. In [START_REF] Spelling | Thermal Modeling of a Solar Steam Turbine With a Forcus on Start-Up Time Reduction[END_REF], three different methods to increase the temperature of the steam turbine in a solar power plant have been investigated: heat blankets on the surface of the housing, gland steam temperature increase and additional insulation. These investigations were conducted on the basis of a simplified calculation model. In [START_REF] Topel | Operational Improvements for Startup Time Reduction in Solar Steam Turbines[END_REF], the work by analyzing three further methods for heating up the turbine with steam has been continued: Gland steam pressure increase, back pressure increase and barring speed increase. Due to its simplicity and possibility of application for different types of turbines the idea of heat blankets attached to the outer surface of the housing gained the interest not only of the university researchers but also of big-league turbomachine producers, like Siemens or GE. To the best of the author's knowledge, both companies have not decided yet to present the final results of their investigation to the international scientific community. Nevertheless it should be mentioned that for two-casing turbine the effectivity of rotor pre-warming may be limited by air heat transfer resistance. Additional problem for the heat blankets method is the pre-warming of internal turbine components like valves. The mentioned methods based on pre-warming with hot working medium very often have to be analyzed under the consideration of windage effects. Windage occurs generally at low volume flow and high rotational speed and causes a temperature increase due to friction and dissipation, as found in literature [START_REF] Evers | Strömungsformen im Ventilationsbetrieb einer ein-und mehrstufig beschaufelten Modellturbine[END_REF], [START_REF] Herzog | Numerische Untersuchung von Schwachlastströmungen in mehrstufigen Axialturbinen[END_REF], [START_REF] Binner | Flow Patterns in High Pressure Steam Turbines During Low-Load Operation[END_REF]. These findings imply that it is possible to heat-up a turbine by windage. A further approach for heating up or keeping warm the turbine with hot air was developed [START_REF] Helbig | A warming arrangement for a steam turbine in a power plant[END_REF]. With an additional hot air warming cycle, warm or hot start-up conditions can be provided even after long idle times. Depending on the expected idle time, hot air is used to keep the HP / IP steam turbine warm or to pre-warm it. This air is heated up by means of an electrical heater and passed with a fan through the turbine, while the rotor is continuously rotating. This approach serves as a basis for the present investigations. Based on the single stage turbine model the heat transfer phenomena in warm-keeping operation have been analyzed in [START_REF] Toebben | Numerical Investigation of the Heat Transfer and Flow Phenomena in an IP Steam Turbine in Warm-Keeping Operation with Hot Air[END_REF]. In the next step the pre-warming operation should be investigated. For this purpose, the additional research into the area of numerical simulations and time-efficient calculation approaches has to be investigated. For the optimization of the startup behavior, calculations of the time-depending temperature distribution within the thick-walled components are required. Due to extremely high computational times necessary for transient Conjugate Heat Transfer (CHT) simulations the fast calculation approaches have to be developed. The numerical methods allowing to calculate heat exchange and transient temperature fields within the turbine in pre-warming operation with regard to thermomechanical fatigue as well as their comparison are described in the following sections.

NUMERICAL MODEL

On the basis of the literature research and in cooperation with the industrial partner General Electric, a numerical model of a single stage turbine is developed to calculate the heat transfer in a steam turbine during pre-warming operation. It is based on the numerical model created for the warm-keeping investigation, which was discussed in [START_REF] Toebben | Numerical Investigation of the Heat Transfer and Flow Phenomena in an IP Steam Turbine in Warm-Keeping Operation with Hot Air[END_REF]. Therefore, in the following sections only the most relevant information concerning model setup as well as the changes introduced to original settings are presented. The numerical model aims to analyze the convective heat transfer and to develop different calculation approaches. As the highest steam temperatures and thermally induced stresses arise particularly within the first steam turbine stage, the heat transfer investigations focus on the first stage. In this Figure 1. Numerical model of the solid body with boundary conditions and surface denotation model, the flow and the heat transfer in the solid state are assumed to be of periodic nature in the circumferential direction. With this assumption, a single blade passage (approx. 2.5°) of the turbine wheel is extracted and meshed. All domains included in the numerical model are discretized by hexahedral cells, which are conformed to each other at the interfaces between fluid state and solid body in the turbine stator and rotor. In the fluid boundary layer a dimensionless wall distance lower than y+ = 1.0 is implemented. In order to exactly calculate the heat transfer, the temperature distribution and the flow field in the boundary layer of the solid body has a high resolution as well. The mesh study conducted in [START_REF] Toebben | Numerical Investigation of the Heat Transfer and Flow Phenomena in an IP Steam Turbine in Warm-Keeping Operation with Hot Air[END_REF] has proven, that spatial discretization has a negligible influence on values of averaged Reynolds-and Nusselt-Number. The fluid state is assumed to be a fully turbulent ideal gas (air) with constant fluid properties. To model the thermal boundary, the Low Reynolds kω-SST turbulence model is used. Constant material properties are applied for the solid state as well.

In general, concerning time dependency of the flow, two types of calculations are required to numerically depict pre-warming process: steady and transient simulations. The main goal of steady state simulations is to generate the boundary conditions for investigated transient calculation approaches. The setup of steady state, as well as of transient models is discussed in the following sections. In Figure 1, the numerical model with its boundary conditions and the used denotations is shown. As presented, the axial walls of the model are considered to be adiabatic. In the present work, the pre-warming procedure has been modelled as thermal shock process, in which the turbine inlet temperature TT,inlet is changed abruptly. The principal behavior of the fluid and solid temperatures in a thermal shock test has been discussed in [START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF] and is shown in Figure 2. As graphically presented, the logarithmic temperature growth of solid temperatures in pre-warming process is on both sides limited by steady state operating points, with constant values of heat fluxes between fluid and solid domains.

STEADY STATE SIMULATIONS

For numerical modelling of pre-warming process the boundary conditions at start t0 and at the end t∞ of procedure have to be determined (Figure 2). In this work, the cold start of prewarming is considered, i.e. at t0 the homogenous temperature T0 in all solid domains is given. Therefore, simple CFD simulation with constant wall temperature and turbine inlet temperature TT,inlet = TFL,∞ is sufficient for subsequent initialization of transient calculation. For the purpose of predefinition of temperature distribution at t0 the CHT simulations can be conducted. In order to generate the missing boundary conditions at t∞ (Figure 2) a mentioned steady state CHT calculation have to be carried out. The setup of this simulation corresponds to the settings of numerical model used for warm-keeping investigations and its description can be found in [START_REF] Toebben | Numerical Investigation of the Heat Transfer and Flow Phenomena in an IP Steam Turbine in Warm-Keeping Operation with Hot Air[END_REF]. In warm-keeping operation (steady state OP2, Figure 2) a constant temperature TR = T(qRotor,OP2) inside the rotor (concentrical hole) and TH = T(qHousing,OP2) at the outer surface of the inner steam turbine casing is defined. The temperature level is approximately 45% of the live steam temperature at nominal load. For all steady state simulations at the inlet of the fluid state, a mass flow rate ṁ and a constant air total temperature TFL,∞ (65% of the live steam temperature at nominal load) are defined as boundary conditions. A constant static pressure is used at the outlet. Due to the applied thermal boundary conditions the impact of radiation is not considered.

TRANSIENT SIMULATIONS

The objective of the work presented within this paper is the development of a fast and generally valid method for the calculation of the transient temperature field in turbine pre-warming with regard to thermomechanical fatigue. As shown in [START_REF] Diefenthal | Speed-up Methods for the Modeling of Transient Temperatures with regard to Thermal and Thermomechanical Fatigue[END_REF] and [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF], the investigated calculation approaches may be used not only in with hot air pre-warmed turbine, but also in numerical depiction of thermal shock tests conducted in other types of turbomachinery like turbochargers. The validated numerical model of turbocharger discussed in [START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF], [START_REF] Diefenthal | Speed-up Methods for the Modeling of Transient Temperatures with regard to Thermal and Thermomechanical Fatigue[END_REF] and [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF] has been used in present investigations in order to verify modifications introduced to extensive reference method in pursuit of calculation time reduction. Generally, even a short literature review reveals the two main simulation approaches of transient fluid/sold heat transfer, which are commonly applied to numerical models. In the first coupled or conjugate heat transfer method (CHT), a CFD calculation is iteratively coupled with a conductive FEA calculation, which results in the solution of the energy equation at the coupled boundaries and provides the continuity of the temperatures as well as the heat fluxes at the fluid/solid interfaces. The practical application of this approach may be found in [START_REF] Bohn | Conjugate Heat Transfer Analysis for Film Cooling Configurations with Different Hole Geometries[END_REF] and [START_REF] Bohn | Conjugate Flow and Heat Transfer Investigation of a Turbocharger[END_REF]. In contradiction to conjugate heat transfer, the uncoupled methods are based on the assumption of a linear relation between the convective heat fluxes across the boundary layer and the driving temperature differences. This assumption may be also understood as the requirement of constant heat transfer coefficients on the boundary interfaces. As a result, the heat transfer coefficients may be determined by means of transient CFD simulation and used as boundary conditions for a subsequent steady state FEA calculation. Examples of this approach are the heat transfer calculations of rotor stator cavities published in [START_REF] Alizadeh | The Use of CFD to Generate Heat Transfer Boundary Conditions for a Rotor-Stator Cavity in a Compressor Drum Thermal Model[END_REF] and [START_REF] Lewis | A Non-Coupled CFD-FE Procedure to Evaluate Windage and Heat Transfer in Rotor-Stator Cavities[END_REF]. The main advantage of the uncoupled approach may be the significantly reduced calculation time, since no iterative solution is necessary at the fluid/solid boundary. Depending on geometry and flow boundary conditions, the uncoupled simulation approaches may lead, however, to high inaccuracies locally, mainly caused by secondary flow phenomena as well as the flow instabilities resulting in nonlinear courses of heat transfer coefficients. In this case, the accuracy of uncoupled approaches may be enhanced by means of CFD and FEA iteration loops, which iteratively increase the quality of results. This solution has been proven, however, to be inextricably linked to two main drawbacks. Firstly, the user input and probability of false setup settings obviously increase depending on the automatization level of numerical models and user experience. Secondly, the computational effort is rapidly changing with each required iteration loop. Owing to the high calculation times, the conjugate heat transfer method is mainly used for steady state calculations. A time accurate calculation based on the coupled method requires an extreme computational effort, since the ratio of convective to conductive heat transfer timescales amount to about 10 -4 , as reported in [START_REF] He | Unsteady Conjugate Heat Transfer Modeling[END_REF]. Consequently, simplifying assumptions have to be made to reach acceptable calculation times for a coupled method simulation of a transient heating or cooling process. In the present paper three numerical methods are applied to a numerical model with hot air pre-warmed steam turbine in order to calculate the transient fluid/solid heat transfer. In a first step extensive simulations are conducted, in which the fluid and the solid calculations are coupled over the whole transient process. Therefore, this simulation approach belongs to a coupled methods family and will be referred to as "Frozen Flow" method (FF), as the fluid equations are partly frozen and not solved over the certain periods of time. A detailed description of this method is given in [START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF] and in the following sections. Additionally, a few modifications are introduced to the original FF method in order to further reduce calculation time with simultaneously almost unaffected solution accuracy. For the purpose of validation, the modified FF method is applied to mentioned numerical model of turbocharger and compared with original data. Based on the results of the FF approach obtained from the numerical model of pre-warmed steam turbine the investigation, points P1 in vane and P2 in blade domains have been identified (Figure 1), which serve as reference for subsequent comparison. The selection criterion for the investigation points is the location of the highest temperature gradients occurring in each analyzed domain in pre-warming operation. In second step an uncoupled approach known as "transient finite element analysis method-exponential" (TFEA-EXPO) method is used in order to model the exponential behavior of the heat transfer coefficients in a transient process. This method is based on commonly applied uncoupled approach and has been developed for purposes of thermomechanical analysis of turbocharger detailed presented in [START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF]. Furthermore, a third approach named Equalized Timescales method (ET) is used for conjugate heat transfer simulation of fluid and solid interaction in pre-warming operation of steam turbine. This calculation time saving approach is based on the modification of solid state properties in pursuit of the reduction of cooling or heating processes. In the following sections, all the three methods are described in detail.

FF method

The commonly used method for an analysis of the fluid/solid conjugate heat transfer is an approach based on [START_REF] Heuer | Numerical and Experimental Analysis of the Thermo-Mechanical Load on Turbine Wheels of Turbochargers[END_REF] and [START_REF] Heuer | Thermomechanical analysis of a turbocharger based on conjugate heat transfer[END_REF]. Due to different time scales of solid and fluid state, and with regard to coupled solution for all domains in each time step, it is assumed that the pressure and the velocity distribution of the fluid remain constant over certain periods of time. With this assumption mass, momentum and turbulence equations of the fluid state do not need to be solved and much larger time steps can be applied. It has been show in [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF] that in the case of thermal shock tests, which are characterized by high heat transfer rates between fluid and solid especially at the very beginning of the transient process, solving only energy equation may lead to higher inaccuracies in solution -6% deviation in temperature gradients courses as reported in [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF].

To consider the changes in the fluid field, the method described in [START_REF] Heuer | Numerical and Experimental Analysis of the Thermo-Mechanical Load on Turbine Wheels of Turbochargers[END_REF], [START_REF] Heuer | Thermomechanical analysis of a turbocharger based on conjugate heat transfer[END_REF] may be extended by updating the fluid field in defined time points or time periods. Two main solution approaches can be distinguished. On one hand, it is possible to conduct a very short fully transient simulation in which all the equations (mass, momentum, energy and turbulence) are solved. This approach was applied to solution procedure described in [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF]. On the other hand, it is also conceivable to replace the fully transient updates simulations with steady state CFD simulations. In the presented paper, the second possibility has been investigated using at first step exactly the same validated model of turbocharger, which served as basis for investigations presented in [START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF], [START_REF] Diefenthal | Speed-up Methods for the Modeling of Transient Temperatures with regard to Thermal and Thermomechanical Fatigue[END_REF] and [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF]. The model of turbocharger used for the FF calculations purposes and for comparison between particular updated approaches comprises the turbine housing, the turbine wheel and the inlet and outlet pipes. As presented in Figure 3, the turbine wheel is modeled as a single rotor passage by means of peripheral averaging (mixing plane) at the rotor stator interfaces. The bearing housing and the compressor are not included in the model. On the boundary plane at the shaft, boundary conditions (α1 + Tref,1, α2 + Tref,2, q) are defined to describe the heat fluxes to the lubricant oil and to the compressor side. These boundary conditions are estimated by using experimental data. The heat shield at the turbine wheel back is assumed to be adiabatic as it is expected that the convective and radiative heat fluxes to the bearing housing are small. The mesh of the whole turbine model consists of approximately 6.2 million nodes. 3.7 million nodes are located in the fluid state and 2.5 million nodes discretize the solid body. In the fluid boundary layer, the dimensionless wall distance y+ is lower than one and the low Reynolds kω-SST turbulence model is applied. A mesh study has been conducted to ensure that the results are unaffected by the mesh. In following investigations the heating up thermal shock process with mass flow m = 0.4 kg/s and fluid temperature TFl,-∞ = 200°C being abruptly changed to TFl,∞ = 600°C (Figure 1) is considered. Further information about numerical setup is provided in [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF]. In pursuit of model validation, the temperatures on the turbine wheel during the thermal shock test are measured. The individual positions of the type K thermocouples are shown in Figure 4. The thermocouples in the turbine wheel are connected to a transmitter, which is located at the shaft end near the compressor inlet. The transmission of amplified signals from the rotating to stationary system is done using an antenna. This antenna also supplies the transmitter with HF-energy. Additional type K thermocouples are attached on the outer surface of the turbine housing to record the temperature distribution. A more detailed description of the test rig is given in [21].

Figure 4. Position of the thermocouples at the turbine wheel

As shown in [START_REF] Diefenthal | Temperature Gradients in a Radial Turbine in Steady State and Transient Operation[END_REF], the original FF method is based on fully transient flow field updates, in which all equations (mass, momentum and energy) are solved for very short periods. These short periods are hereafter referred to as "Transient Update" periods. The simplified periods are named "Only Energy" periods. In the case of turbocharger three "All Equation" periods at 0 s, 6 s, and 60 s are calculated for numerical depiction of the thermal shock process. An alternative possibility for fluid flow updates are steady state CFD simulations, which allow adjustment of pressure and velocity fields to calculated temperature values in domains. This fluid field update method is referred to as "Stationary Update". The schematic representation of FF Stationary Update method is provided in Figure 5. At the beginning of simulation procedure, a steady state CHT All Equation 1 In pursuit of FF Stationary Update method validation, the temperature values obtained from both update approaches are compared with experimental data in measuring point MP3. As depicted in Figure 7, the calculated temperature values show good agreement with experimental values. In order to evaluate the temperature deviation in all measurement points on turbine wheel, the numerical and experimental data are normalized to a time step of 1.0 s by a cubical interpolation. Thereby, a root mean square deviation is determined over a period of 800 s:

       n 1 i 2 calc exp av T T n 1 T (2)
The obtained values of RMS deviations are presented in Figure 8. The FF Transient Update method, as well as the FF Stationary Update method, provide the temperature time plots that correlate well with experimental data. Hence, regarding computational time savings, in steam turbine pre-warming numerical calculations as the reference method the FF Stationary Update method will be used. Thus, for the interpolation of the heat transfer coefficients it is assumed, that those linearly depend on the wall temperatures. With regard to the exponential behavior of the wall temperatures this interpolation method results in an exponential time plot of the heat transfer coefficients over the time. Hence, the method is abbreviated with TFEA-EXPO. The following equations are applied to interpolate the heat transfer coefficients for each cell of the prewarmed steam turbine:

  0 , i , i 0 , i i ) t ( I ) t (         (3)     0 , i , i 0 , i i T T / T ) t ( T ) t ( I     (4)
Equation 4 Ti (t) describes the time-dependent wall temperature. This temperature is directly used from the transient FEA calculation. Hence, the boundary conditions depend on the solution of the FEA simulation and are calculated iteratively in every time step. The variables Ti,0 and Ti,∞ may be determined from steady state CHT calculations conducted respectively at the beginning and at the end of thermal shock process (t → ∞). Similarly, the variable αi,∞ in Equation 3 describes the heat transfer coefficient derived from steady state CHTt→∞. Finally, the variable αi,0 refers to heat transfer coefficient at the beginning of the thermal shock process immediately after switching the inlet temperature of the turbine from TFl,-∞ to TFl,∞ (Figure 1). The further information about TFEA-EXPO method may be found in [START_REF] Diefenthal | Speed-up Methods for the Modeling of Transient Temperatures with regard to Thermal and Thermomechanical Fatigue[END_REF].

ET-method

The Equalized Timescales method may be classified just like FF approach as a conjugate heat transfer method. The ET method is a time saving approach, which similarly to the other presented simulation procedures may be used to calculate the transient temperature fields with regard to thermo-mechanical fatigue. Correspondingly, in this method the fluid and the solid state are calculated coupled with the whole transient process. However, the fluid calculation is not simplified like in the FF approach. Thus, a method characterized by high accuracy and small user effort is expected. The main obstacle to overcome remains unchanged: in case of a coupled calculation the ratio of timescales in fluid to timescales in solid state amounts to about 10 -4 , as given in [START_REF] He | Unsteady Conjugate Heat Transfer Modeling[END_REF]. These different timescales cause long conductional heat transfer processes in the solid state that have to be calculated with small time steps to achieve numerical stability of the fluid calculation. In order to reduce the calculation times in this method, the timescales of fluid and solid are equalized by modifying the material properties of the solid state. Thus, this method is named Equalized Timescales method (ET). To describe the approach of the ET method in more detail, a simple analytical approach is considered in the following sections. For a solid body with a homogenous temperature distribution, the timescale of a convectional heating or cooling process can be determined from a simple analytical equation:

  FL SL SL SL SL SL p SL T T A t T V c        , (5)
With the assumptions of constant material properties, a constant heat transfer coefficient and a constant fluid temperature over the time, the solution of Equation 5 leads to Equation 6:

    FL t FL SL FL t V c A FL SL SL T e T T T e T T T SL SL SL p SL            0 , 0 , , (6)
The exponent of the Euler's Number represents the timescale τ of the heating or the cooling process. It depends on the heat transfer coefficient α, the density ρSl and the specific heat capacity cp,Sl. In the ET method, the timescales of the fluid and solid are equalized by reducing the specific heat capacity of the solid state by the speed up factor SF. The modified specific heat capacity is named cp* and is defined by the following relation:

SF c c SL p SL p , * ,  (7)
From 6 it can be derived that the reduction of the specific heat capacity results in a reduced heating or cooling time of the solid body. The modified time is named t* in the present work and is connected with the real time by the speed up factor SF for this simplified case. The definition of the modified time is given in Equation 8:

SF t t  * (8)
Implementing this approach to a coupled fluid/solid heat transfer simulation results in significantly reduced calculation times as shorter conductional heating or cooling processes are calculated and the time step of the fluid calculation can remain constant. To determine the real heating or cooling time of a transient process from these simulations, the resulting time t* is multiplied by the speed up factor SF according to Equation 8. The Fourier-and the Biot-Number are assessed in order to prove the suitability of the described approach for the determination of the thermal behavior of the solid state in a transient process. As described in [22], transient temperature fields of the same geometry and the same boundary conditions are similar if the Biot-Number and the Fourier-Number are equal. Thus, the Biot-and the Fourier-Number must remain constant when the specific heat capacity is modified. In Equation 9, the modified specific heat capacity cp* and the modified time t* are introduced in the definition of the Fourier-Number:

* 2 * , * 2 * , * 2 , Fo L c t L SF c SF t L c t Fo SL SL p SL SL SL p SL SL SL p SL           (9)
As shown, the unmodified and the modified Fourier-Numbers are identical for the described simplified approach. Unlike the Fourier-Number, the Biot-Number is independent of the specific heat capacity and of the time. However, the Biot-Number depends on the heat transfer coefficient:

* Bi L Bi SL     (10)
With the assumption of a constant heat transfer coefficient, the Biot-Number remains constant. To investigate the assumption of a constant heat transfer coefficient in more detail and to verify the ET method in case of numerical modeling of pre-warming processes, all the methods are applied to the numerical model of pre-warmed steam turbine. For the ET method, a speed up factor of 1000 is used.

COMPARISON OF NUMERICAL METHODS FOR PRE-WARMING CALCULATION

The main goal that the TFEA-EXPO and the ET approaches strive for is the possibility to accurately calculate the temperature gradients resulting from 3D flow structures with significantly reduced computational and user effort in comparison with the FF method.

In pursuit of verification and comparison of the methods, a steam turbine pre-warming process modeled as thermal shock has been calculated with all three simulation approaches. It is assumed, that at the beginning of the pre-warming, all solid domains in investigated 1-stage turbine model have homogenous temperature distribution equal to TSL,-∞. Similarly to thermal shock process in the discussed turbocharger, the fluid inlet temperature is being abruptly changed from TFL,-∞ to TFL,∞ at time point t0. In a second step, the calculated temperatures at the investigation points P1 and P2 are compared with reference approach. As shown in Figure 10 all numerical methods provide temperature plots which show a good agreement with the reference FF method.

In order to quantify the respective error, the root mean square (RMS) deviation between measured and calculated temperatures is determined over a period of the first 1000 s for a time step of 1.0 s. The RMS deviations for both investigation points P1 and P2, comparison of the computing times (calculated for first 1000 s of pre-warming procedure) as well as the temperature distribution in investigated single stage turbine model are presented in Figure 11.

CONCLUSION

In the present work two calculation methods -TFEA-EXPO and ET for modeling of transient temperature fields in prewarmed steam turbine are conducted and compared with an extensive calculation approach for verification purposes. Additionally, some changes resulting in reduced computational time are introduced into numerical procedure of an extensive FF method. Moreover, in a subsequent step the modified FF method is validated against the original FF approach and experimental data using known from several other publications presenting the numerical model of turbocharger.

To conclude, all the methods presented in the present work are suitable to calculate the transient temperature field with regard to thermomechanical fatigue within the timescales of a standard design process. Furthermore, the additional savings in computational time have been achieved -89% in case of TFEA-EXPO approach. The ET method has been analyzed only for one value of speed up factor (SF = 1000) and the good accuracy with simultaneously minimum user effort has been achieved. In future research the more computing timeefficient ET simulations with higher values of speed up factor should be investigated.

Figure 2 .

 2 Figure 2. Qualitative idealized behavior of fluid and solid temperatures in a thermal shock test[START_REF] Diefenthal | Thermomechanical Analysis of Transient Temperatures in a Radial Turbine Wheel[END_REF]

Figure 3 .

 3 Figure 3. Numerical model of Turbine Housing (TH) and Turbine Wheel (TW)

  or CFD All Equation 1 calculation (depending on solid body initialization) is conducted. The outgoing from the generated state of fluid the transient CHT Only Energy 1 simulation (Figure 5) has to be calculated until discrete time point of update simulation has been achieved. The temperature TSolid from the CHT Only Energy 1 simulation are used in the next step as the boundary conditions for subsequent CFD All Equation 2 calculation. Finally, in the CHT Only Energy 2 simulation the fluid state is being initialized by the results of CFD All Equation 2 and the temperature distribution of solid domains is set to be equal to these from CHT Only Energy 1 calculation. The presented CFD update procedure may be repeated depending on required number of fluid flow updates.

Figure 5 .Figure 6 .

 56 Figure 5. Schematic representation of FF Stationary Update methodIn order to evaluate the impact of the updates on the fluid field the averaged relative Reynolds number for turbine wheel is determined as follows:

Figure 7 .

 7 Figure 7. Temperature at measuring point 3 over time for a heating up process (TFl,-∞ = 200°C → 600°C at m = 0.4 kg/s)

Figure 8 .

 8 Figure 8. RMS deviation between measured and calculated temperatures at the measuring points 3.2 TFEA-EXPO As mentioned in the previous sections, the TFEA-EXPO method is an uncoupled approach based on transient FEA simulation. In order to model the transient fluid field, heat transfer coefficients at the start and end of the transient process are derived from steady state CFD or CHT calculations. The obtained heat transfer coefficients are interpolated over the time span of transient process for each cell element. After the switch of the turbine fluid inlet temperature from TFl,-∞ to TFl,∞ the inlet conditions and the outlet pressure remain nearly constant during the transient process. Thereby, the changes of the fluid state result mainly from the changes of heat fluxes between fluid/solid interfaces as well as from the wall temperatures.Thus, for the interpolation of the heat transfer coefficients it is assumed, that those linearly depend on the wall temperatures. With regard to the exponential behavior of the wall temperatures this interpolation method results in an exponential time plot of the heat transfer coefficients over the time. Hence, the method is abbreviated with TFEA-EXPO. The following equations are applied to interpolate the heat transfer coefficients for each cell of the prewarmed steam turbine:

Figure 9 .

 9 Figure 9. Heat transfer coefficients on the surface of vane and bladeIn first step the heat transfer coefficients are derived from the FF calculations and used for subsequent comparison with HTCs obtained from the TFEA-EXPO as well as from the ET simulation. The area averaged heat transfer coefficients on the surface of the vane and blade are depicted in Figure9. As expected, all the plots of the heat transfer coefficients prove an exponential behavior. The oscillating plot of HTCs on blade indicates, that the Speed up factor in the ET method cannot be increased without any further limitations. These oscillations result from unsteady flow effects which can have a large influence on the solid temperatures for higher speed up factors due to the similar timescales of fluid and solid state.

Figure 10 .

 10 Figure 10. Temperatures calculated at characteristic points P1 and P2

Figure 11 .

 11 Figure 11. Temperature distribution, RMS deviation and overview over the calculation times in [core-h]

ACKNOWLEDGEMENTS

The present work is a part of the joint research program COOREFELX-Turbo in the frame of AG Turbo. This work is supported by German ministry of economy and technology under grant number 03ET7041B. The authors would like to gratefully acknowledge GE Electric whose support was essential for the presented investigations. The responsibility for the content of this publication lies with the corresponding authors.


Ümit V C ¸atalyürek

Mehmet Deveci

Kamer Kaya

Bora Uçar

UMPa: A Multi-objective, multi-level partitioner for communication minimization

Keywords:

We propose a directed hypergraph model and a refinement heuristic to distribute communicating tasks among the processing units in a distributed memory setting. The aim is to achieve load balance and minimize the maximum data sent by a processing unit. We also take two other communication metrics into account with a tie-breaking scheme. With this approach, task distributions causing an excessive use of network or a bottleneck processor which participates in almost all of the communication are avoided. We show on a large number of problem instances that our model improves the maximum data sent by a processor up to 34% for parallel environments with 4, 16, 64 and 256 processing units compared to the state of the art which only minimizes the total communication volume.

Introduction

In parallel computing, the problem of distributing communicating tasks among the available processing units is important. To solve this problem, several graph and hypergraph models are proposed [START_REF] Aykanat | A hypergraph model for mapping repeated sparse matrixvector product computations onto multicomputers[END_REF][START_REF]Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF][START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF][START_REF] Hendrickson | Graph partitioning models for parallel computing[END_REF][START_REF] Walshaw | Parallel dynamic graph partitioning for adaptive unstructured meshes[END_REF]. These models transform the problem at hand to a balanced partitioning problem. The balance restriction on part weights in conventional partitioning corresponds to the load balance in the parallel environment, and the minimized objective function corresponds to the total communication volume between processing units. Both criteria are crucial in practice for obtaining short execution times, using less power, and utilizing the computation and communication resources better.

In addition to the total data transfer, there are other communication metrics investigated before, e.g., total number of messages sent [START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF], or maximum volume of messages sent and/or received by a processor [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF]. Even with perfect load balancing and minimized total data transfer, there can be a bottleneck processing unit which participates in most of the data transfers. This can create a problem especially for data intensive applications where reducing the amount of data transferred by the bottleneck processing unit can improve the total execution time significantly.

In this work, given a task graph, our main objective is distributing its tasks evenly and minimizing the maximum amount of data sent by a processing unit. Previous studies addressing different communication cost metrics (such as [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF]) work in two phases. In the first phase, the total volume of communication is reduced, and in the second phase the other metrics are addressed. We propose a directed hypergraph model and partition the related hypergraph with a multi-level approach and a novel K-way refinement heuristic. While minimizing the primary objective function, our refinement heuristic also takes the maximum data sent and received by a processing unit and the total amount of data transfer into account by employing a tie-breaking scheme. Therefore, our approach is different from the existing studies in that the objective functions are minimized all at the same time.

The organization of the paper is as follows. In Section 2, the background material on graph and hypergraph partitioning is given. Section 2.3 shows the differences of the graph and hypergraph models and describes the proposed directed hypergraph model. In Section 3, we present our multi-level, multi-objective partitioning tool UMPa (pronounced as "Oompa"). Section 4 presents the experimental results, and Section 5 concludes the paper.

Background

Hypergraph partitioning.

A hypergraph H = (V, N ) is defined as a set of vertices V and a set of nets (hyperedges) N among those vertices. A net n ∈ N is a subset of vertices and the vertices in n are called its pins. The number of pins of a net is called the size of it, and the degree of a vertex is equal to the number of nets it is connected to. In this paper, we will use pins[n] and nets[v] to represent the pin set of a net n and the set of nets vertex v is connected to, respectively. Vertices can be associated with weights, denoted with w[•], and nets can be associated with costs, denoted with c[•].

A K-way partition of a hypergraph H is denoted as Π = {V 1 , V 2 , . . . , V K } where

• parts are pairwise disjoint, i.e., V k ∩ V = ∅ for all 1 ≤ k < ≤ K, • each part V k is a nonempty subset of V, i.e., V k ⊆ V and V k = ∅ for 1 ≤ k ≤ K, • union of K parts is equal to V, i.e., K k=1 V k = V. Let W k denote the total vertex weight in V k (i.e., W k = v∈V k w[v]
) and W avg denote the weight of each part when the total vertex weight is equally distributed (i.e.,

W avg = ( v∈V w[v])/K). If each part V k ∈ Π satisfies the balance criterion (2.1) W k ≤ W avg (1 + ε), for k = 1, 2, . . . , K
we say that Π is -balanced where ε represents the maximum allowed imbalance ratio.

For a K-way partition Π, a net that has at least one pin (vertex) in a part is said to connect that part. The number of parts connected by a net n, i.e., connectivity, is denoted as λ n . A net n is said to be uncut (internal) if it connects exactly one part (i.e., λ n = 1), and cut (external), otherwise (i.e., λ n > 1).

The set of external nets of a partition Π is denoted as N E . There are various cutsize definitions [START_REF] Lengauer | Combinatorial algorithms for integrated circuit layout[END_REF] for hypergraph partitioning. The one that will be used in this work, which is shown to accurately model the total communication volume [START_REF]Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], is called the connectivity metric and defined as:

χ(Π) = n∈N c[n](λ n -1) . (2.2)
In this metric, each cut net n contributes c[n](λ n -1) to the cutsize. The hypergraph partitioning problem can be defined as the task of finding a balanced partition Π with K parts such that χ(Π) is minimized. This problem is also NP-hard [START_REF] Lengauer | Combinatorial algorithms for integrated circuit layout[END_REF].

2.2. K-way partitioning and multi-level framework. Arguably, the multilevel approach [START_REF] Barnhard | Fast multilevel implementation of recursive spectral bisection for partitioning unstructured problems[END_REF] is the most successful heuristic for the hypergraph partitioning problem. Although, it has been first proposed for recursive-bisection based graph partitioning, it also works well for hypergraphs [START_REF] Cevdet Aykanat | Multi-level direct k-way hypergraph partitioning with multiple constraints and fixed vertices[END_REF][START_REF] Bui | A heuristic for reducing fill-in sparse matrix factorization[END_REF][START_REF]Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF][START_REF] Hendrickson | A multilevel algorithm for partitioning graphs[END_REF][START_REF] Trifunovic | Parkway 2.0: A parallel multilevel hypergraph partitioning tool[END_REF]. In the multi-level approach, a given hypergraph is coarsened to a much smaller one, a partition is obtained on the the smallest hypergraph, and that partition is projected to the original hypergraph. These three phases will be called the coarsening, initial partitioning, and uncoarsening phases, respectively. The coarsening and uncoarsening phases have multiple levels. In a coarsening level, similar vertices are merged to make the hypergraph smaller. In the corresponding uncoarsening level, the merged vertices are split, and the partition of the coarser hypergraph is refined for the finer one.

Most of the multi-level partitioning tools used in practice are based on recursive bisection. In recursive bisection, the multi-level approach is used to partition a given hypergraph into two. Each of these parts is further partitioned into two recursively until K parts are obtained in total. Hence, to partition a hypergraph into K = 2 k , the recursive bisection approach uses K -1 coarsening, initial partitioning, and uncoarsening phases.

Several successful clustering heuristics are proposed to coarsen a hypergraph. Although their similarity metrics aim to reduce the cutsize, they cannot find an optimal solution, since the problem is NP-hard. Hence, an optimal partition of the coarser hypergraph may not be optimal for the finer one. To obtain better partitions, iterative-improvement-based heuristics are used to refine the coarser's partition after projecting it to finer. In practice, Kernighan-Lin (KL) [START_REF] Kernighan | An efficient heuristic procedure for partitioning graphs[END_REF] and Fiduccia-Mattheyses (FM) [START_REF] Fiduccia | A linear-time heuristic for improving network partitions[END_REF] based refinement heuristics that depend on vertex swaps and moves between two parts are used.

Task graph and communication volume metrics.

Let A = (T , C) be a task graph where T is the set of tasks to be executed, and C is the set of communications between pairs of tasks. We assume that the execution time of each task may differ, hence each task t ∈ T is associated with an execution time exec(t). Each task t i ∈ T sends a different amount of data data(t i ) to each t j such that t i t j ∈ C. The communications between tasks may be uni-directional, That is t i t j ∈ C does not imply t j t i ∈ C. In our parallel setting, we assume owner computes rule and hence, each task of A is executed by the processing unit to which it is assigned. Let T i be the set of tasks assigned to processing unit P i . Since it is desirable to distribute the tasks evenly, the computational load t∈Ti exec(t) should be almost the same for each P i . In addition to that two heavily communicating tasks should be assigned to the same processing unit since less data transfer over the network is needed in this case. The total amount of data transfer throughout the execution of the tasks is called the total communication volume (totV ). Note that when a task t ∈ T i needs to send data to a set of tasks in T j , the contribution to totV is data(t), since it is enough to send t's data to P j only once.

Although minimizing the total communication volume is important, it is sometimes preferable to reduce other communication metrics [START_REF] Hendrickson | Graph partitioning models for parallel computing[END_REF]. For example, in the context of one-dimensional partitioning of structurally unsymmetric sparse matrices for parallel matrix-vector multiplies, Uçar and Aykanat used a communication hypergraph model to reduce the maximum of number of messages and the maximum amount of data sent and received by a processor [START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF] (see also [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF] and [START_REF] Uçar | Minimizing communication cost in fine-grain partitioning of sparse matrices[END_REF] for other communication metrics).

Let SV [i] and RV [i] be the volumes of communication sent and received by P i , respectively. Hence, the total communication volume equals to totV

= i SV[i] = i RV[i].
In addition to totV , we are interested in two other communication metrics: maximum send volume (maxSV ), which equals to max i (SV[i]); and maximum send-receive volume (maxSRV ), which is max i (SV[i] + RV[i]).

UMPa: A multi-objective partitioning tool for communication minimization

3.1. Directed hypergraph model. We propose modeling the task graphs with directed hypergraphs. Given a task graph A, we construct the directed hypergraph model H = (V, N ) as follows. For each task t i ∈ T , we have a corresponding vertex v i ∈ V and a net n i ∈ N where pins

[n i ] = {v i } ∪ {v j | t i t j ∈ C}, w[v i ] = exec(t i ), and c[n i ] = data(t i ).
In this directed hypergraph model, the communication represented by a net n is flowing from its source vertex, which will be denoted as s(n), to the target vertices pins[n] \ {s(n)}. Given a partition Π, let δ(n, V i ) = 1 if n ∩ V i = ∅, and 0, otherwise. Then the data sent and received by

P i are equal to SV[i] = n,s(n)∈Vi c[n](λ n -1) and RV[i] = n,s(n) /
∈Vi c[n]δ(n, V i ), respectively. Our primary objective is to minimize maxSV , the maximum send volume. While doing this, we also take the maximum send-receive volume and the total communication volume into account. The total volume of communication corresponds to the cutsize definition (2.2) as in the standard hypergraph model. In other words, the sense of direction is not important for the total communication volume totV . On the other hand, the directions of the flow is crucial while minimizing maxSV and maxSRV .

To optimize its metrics, UMPa follows the multi-level approach. Instead of a recursive bisection, it adopts a direct K-way partitioning. Given the hypergraph, UMPa gradually coarsens it, obtains an initial K-way partition for the coarsest hypergraph, and projects it into the original one by uncoarsening and refinement steps at each level.

Multi-level coarsening phase.

In this phase, the original hypergraph is gradually coarsened in multiple levels by clustering subsets of vertices at each level. There are two types of clustering algorithms: matching-based ones and agglomerative ones. The matching-based ones put at most two similar vertices in a cluster, whereas the agglomerative ones allow any number of similar vertices. There are various similarity metrics-see for example [START_REF] Alpert | Recent directions in netlist partitioning: A survey[END_REF][START_REF]Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF][START_REF] Karypis | Multilevel hypergraph partitioning[END_REF]. All these metrics are defined only on two adjacent vertices (one of them can be a vertex cluster). Two vertices are adjacent if they share a net and they can be in the same cluster if the are adjacent.

In this work, we use an agglomerative algorithm and the absorption clustering metric using pins [START_REF] Alpert | Recent directions in netlist partitioning: A survey[END_REF][START_REF] Aykanat | Patoh: A multilevel hypergraph partitioning tool[END_REF]. For this metric, the similarity between two adjacent vertices u and v is

n∈nets[u]∩nets[v] c[n] |pins[n]| -1
This is also the default metric in PaToH [START_REF] Aykanat | Patoh: A multilevel hypergraph partitioning tool[END_REF], a well-known hypergraph partitioner. In each level , we start with a finer hypergraph H and obtain a coarser one H +1 . If V C ⊂ V is a subset of vertices deemed to be clustered, we create the cluster vertex u ∈ V +1 where nets[u] = ∪ v∈V C nets [v]. We also update the pin sets of the nets in nets[u] accordingly.

Since we need the direction, i.e., source vertex information for each net to minimize maxSV and maxSRV , we always store the source vertex of a net n ∈ N as the first pin in pins[n]. To maintain this information, when a cluster vertex u is formed in the coarsening phase, we put u to the head of pins[n] for each net n whose source vertex is in the cluster.

Initial partitioning phase.

To obtain an initial partition for the coarsest hypergraph, we use PaToH [START_REF] Aykanat | Patoh: A multilevel hypergraph partitioning tool[END_REF], which is proven to produce high quality partitions with respect to total communication volume metric [START_REF]Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF]. We execute PaToH ten times and get the best partition according to the maxSV metric. We have several reasons to use PaToH. First, although our main objective is minimizing maxSV , since we also take totV into account, it is better to start with an initial partition having a good total communication volume. Second, since totV is the sum of the send volumes of all parts, as we observed in our preliminary experiments, minimizing it may also be good for both maxSV and maxSRV . Also, as stated in [START_REF] Cevdet Aykanat | Multi-level direct k-way hypergraph partitioning with multiple constraints and fixed vertices[END_REF], using recursive bisection and FM-based improvement heuristics for partitioning the coarsest hypergraph is favorable due to small net sizes and high vertex degrees.

K-way refinement of communication volume metrics.

In an uncoarsening level, which corresponds to the th coarsening level, we project the partition Π +1 obtained for H +1 to H . Then, we refine it by using a novel K-way refinement heuristic which is described below.

Given a partition Π, let a vertex be a boundary vertex if it is in the pin set of at least one cutnet. Let Λ(n, p) = |pins[n] ∩ V p | be the number of pins of net n in part p, and part[u] be the current part of u. The proposed heuristic runs in multiple passes where in a pass it visits each boundary vertex u and either leaves it in part[u], or moves it to another part according to some move selection policy. Algorithm 1 shows a pass of the proposed refinement heuristic. For each visited boundary vertex u and for each available part p other than part[u], the heuristic computes how the communication metrics are affected when u is moved to p. This is accomplished in three steps. First, u is removed from part[u], and the leave gains on the send/receive volumes of the parts are computed (after line 1). Second, u is put into a candidate part p and the arrival losses on the send/receive volumes are computed (after line 2). Last, the maximum send, maximum send-receive, and total volumes are computed for this move (after line 4).

3.4.1. Move selection policy and tie-breaking scheme. Our move selection policy given in Algorithm 2 favors the moves with the maximum gains on maxSV and never allows a move with negative gain on the same metric. To take other metrics into account, we use a tie-breaking scheme which is enabled when two different moves of a vertex u have the same maxSV gain. In this case, the move with maxSRV gain is selected as the best move. If the gains on maxSRV are also equal then the move with maximum gain on totV is selected. We do not allow a vertex move without a positive gain on any of the communication metrics. As the experimental results show, this move selection policy and tie-breaking scheme have positive impact on all the metrics. Figure 1 shows a sample graph with 8 vertices and 13 edges partitioned into 3 parts. Assume that this is a partial illustration of boundary vertices, and any move will not violate the balance criteria. Each row in the table contains a possible vertex move and the changes on the communication volume metrics. In the initial configuration, maxSV = 6, maxSRV = 9, and totV = 12. If we move v 3 from the partition V 2 to the partition V 3 , we reduce all metrics by 1. On the other hand, if we move v 3 to V 1 , we decrease maxSV and maxSRV , but totV does not change. In this case, since its gain on totV is better, the tie-breaking scheme favors the move v 3 to V 3 . Moreover, the moves v 4 to V 1 , v 6 to V 3 and v 7 to V 3 are other move examples where tie-breaking scheme is used. Note that we allow all the moves in the first 13 rows of the table including these two. However, we do not allow the ones in the last three rows.

3.4.2. Implementation details. During the gain computations, the heuristic uses the connectivity information between nets and parts stored in data structures λ and Λ. These structures are constructed after the initial partitioning phase, and then maintained by the uncoarsening phase. Since the connectivity information changes after each vertex move, when a vertex u is moved, we visit the nets of u and update the data structures accordingly. Also, when new vertices become boundary vertices, they are inserted into boundary array and visited in the same pass.

Vertex Part maxSV maxSRV totV v1 V1 -1 +1 -2 v1 V2 -2 -2 -3 v2 V2 0 -1 -1 v2 V3 -1 +1 0 v3 V1 -1 -1 0 v3 V3 -1 -1 -1 v4 V1 -1 -1 0 v4 V3 -1 +1 +1 v5 V3 0 0 -1 v6 V1 -1 0 +1 v6 V3 -1 0 0 v7 V1 -1 +1 0 v7 V3 -1 -1 0 v5 V2 +2 +2 -1 v8 V1 0 0 0 v8 V2 +2 +2 +1 Figure 1.
A sample partitioning and some potential moves with their effects on the communication volume metrics. The initial values are maxSV = 6, maxSRV = 9 and totV = 12. A negative value in a column indicates a reduction on the corresponding metric.

If at least one move with a positive gain on maxSV is realized in a refinement pass, the heuristic continues with the next pass. Otherwise, it stops. For efficiency purposes, throughout the execution of a pass, we restrict the number of moves for each vertex u. If this number is reached, we lock the vertex and remove it from the boundary. In our experiments, the maximum number of moves per vertex is 4.

Let ρ = n∈N |pins[n]| be the number of pins in a hypergraph. The time complexity of a pass of the proposed refinement heuristic is O(ρK + |V|K 2 ) due to the gain computation loops at lines 3 and 5. To store the numbers of pins per part for each net, Λ, we use a 2-dimensional array. Hence, the space complexity is O(K|N |). This can be improved as shown in [START_REF] Cevdet Aykanat | Multi-level direct k-way hypergraph partitioning with multiple constraints and fixed vertices[END_REF].

Experimental results

UMPa is tested on a computer with 2.27GHz dual quad-core Intel Xeon CPUs and 48GB main memory. It is implemented in C++ and compiled with g++ version 4.5.2.

To obtain our data set, we used several graphs from the testbed of 10th DI-MACS implementation challenge [START_REF]10th DIMACS implementation challenge: Graph partitioning and graph clustering[END_REF]. We remove relatively small graphs containing less than 10 4 vertices, and also extremely large ones. There are 123 graphs in our data set from 10 graph classes. For each graph, we execute UMPa and other algorithms 10 times. The results in the tables are the averages of these 10 executions.

To see the effect of UMPa's K-way partitioning structure and its tie-breaking scheme, we compare it with two different refinement approaches and PaToH. The first approach is partitioning the hypergraph into K with PaToH's recursive bisection scheme and refining it by using the proposed K-way refinement algorithm without employing the tie-breaking scheme. The second approach is using UMPa but again without tie breaking. To remove tie breaking, we remove the else statements at lines labeled with 1 and 2 of Algorithm 2. Table 1 gives the average performance of all these approaches normalized with respect to PaToH's performance. Without tie breaking, refining PaToH's output reduces the maximum send volume by 8%. However, it increases the maximum send-receive and total volumes by 5% and 3%, respectively. Hence, we do not suggest using the refinement heuristic alone and without tie breaking. On the other hand, if it is used in the multi-level structure of UMPa, we obtain better results even without a tie-breaking scheme.

Table 1 shows that UMPa's multi-level structure helps to obtain 17% and 7% less volumes than PaToH's partitions in terms of maxSV and maxSRV , respectively. But since PaToH minimizes the total communication volume, there is a 6% overhead on the totV . Considering 17% reduction on maxSV , this overhead is acceptable. However, we can still reduce all the communication metrics 9%-to-10% more by employing the proposed tie-breaking scheme. For K = 4, this leads us a 34% better maximum send volume, which is impressive since even the total communication volume is 16% less compared with PaToH. Actually, for all K values, UMPa manages to reduce maxSV and maxSRV on the average. The percent of improvement reduces with the increasing K. This may be expected since when K is large, the total volume will be distributed into more parts, and the maximum send or send-receive volume will be less. Still, on the average, the reductions on maxSV , maxSRV , and totV are 26%, 16%, and 4%, respectively.

Tables 2 and3 show performance of PaToH and UMPa in terms of the communication metrics and time. There are 20 graphs in each table selected from 10 graph class in DIMACS testbed. For each graph class, we select the two (displayed consecutively in the tables) for which UMPa obtains the best and worst improvements on maxSV . The numbers given in the tables are averages of 10 different executions. For all experiments with K = 16 parts, as Table 2 shows, UMPa obtains a better maxSV value than PaToH on the average. When K = 4, 64, and 256, PaToH obtains a better average maxSV only for 16, 4, and 1 graphs, out of 123, respectively.

There are some instances in the tables for which UMPa improves maxSV significantly. For example, for graph ut2010 in Table 2, the maxSV value is reduced from 1506 to 330 with approximately 78% improvement. Furthermore, for the same graph, the improvements on maxSRV and totV are 75% and 67%, respectively. When K = 256 (Table 3) for the graph memplus, UMPa obtains approximately 50% improvement on maxSV and maxSRV . Although totV increases 26% at the same time, this is acceptable considering the improvements on the first two metrics. Table 4 shows the relative performance of UMPa in terms of execution time with respect to PaToH. As expected, due to the complexity of K-way refinement heuristic, UMPa is slower than PaToH especially when the number of parts is large.

Conclusions and future work

We proposed a directed hypergraph model and a multi-level partitioner UMPa for obtaining good partitions in terms of multiple communication metrics where the maximum amount of data sent by a processing unit is the main objective function to be minimized. UMPa uses a novel K-way refinement heuristic employing a tiebreaking scheme to handle multiple communication metrics. We obtain significant improvements on a large number of graphs for all K values.

We are planning to speed up UMPa and the proposed refinement approach by implementing them on modern parallel architectures. We are also planning to investigate partitioning for hierarchical memory systems, such as cluster of multisocket, multi-core machines with accelerators.

Table 1 .

 1 The relative performance of UMPa and Pa-ToH+refinement without tie breaking. The performance are computed with respect to that of PaToH.

		PaToH + Refinement		UMPa		UMPa
		No tie breaking	No tie breaking	With tie breaking
	K	maxSV maxSRV totV maxSV maxSRV totV maxSV maxSRV totV
	4	0.93	1.05 1.06	0.73	0.83 0.93	0.66	0.77 0.84
	16	0.93	1.06 1.04	0.84	0.94 1.11	0.73	0.83 0.98
	64	0.91	1.04 1.02	0.86	0.98 1.12	0.76	0.87 1.00
	256	0.91	1.03 1.01	0.89	1.00 1.10	0.81	0.91 1.02
	Avg.	0.92	1.05 1.03	0.83	0.93 1.06	0.74	0.84 0.96

Table 2 .

 2 The maximum send and send-receive volumes, and the total volume for PaToH and UMPa when K = 16. The times are given in seconds. There are 20 graphs in the table where two graphs with the best and the worst improvements on maxSV are selected from each class. Each number is the average of 10 different executions.

			PaToH			UMPa
	Graph	maxSV maxSRV	totV Time maxSV maxSRV	totV Time
	coPapersDBLP	62,174 139,600 673,302 91.45 53,619 117,907 842,954 145.47
	as-22july06	1,506	5,063 12,956 0.63	1,144	3,986 13,162 2.70
	road central	500	999 3,926 112.64	279	576 2,810 27.85
	smallworld	12,043	24,020 188,269 3.09 10,920	21,844 174,645 19.27
	delaunay n14	119	235 1,500 0.19	115	236 1,529 0.88
	delaunay n17	351	706 4,100 1.09	322	655 4,237 2.54
	hugetrace-00010	2,113	4,225 25,809 93.99	2,070	4,144 28,572 43.39
	hugetric-00020	1,660	3,320 20,479 60.96	1,601	3,202 22,019 29.51
	venturiLevel3	1,774	3,548 19,020 27.41	1,640	3,282 20,394 16.01
	adaptive	2,483	4,967 27,715 54.00	2,345	4,692 29,444 29.33
	rgg n 2 15 s0	146	293 1,519 0.34	119	254 1,492 1.03
	rgg n 2 21 s0	1,697	3,387 19,627 37.86	1,560	3,215 20,220 16.66
	tn2010	2,010	3,666 13,473 1.26	1,684	3,895 56,780 1.54
	ut2010	1,506	2,673 3,977 0.43	330	677 1,303 0.82
	af shell9	1,643	3,287 17,306 14.83	1,621	3,242 18,430 8.64
	audikw1	15,119	29,280 145,976 161.23 11,900	24,182 159,640 77.16
	asia.osm	63	125	409 40.43	30	62	323 7.67
	belgium.osm	141	281 1,420 4.80	120.6	243 1,406 1.96
	memplus	986	7,138 7,958 0.23	686	3,726 10,082 0.72
	t60k	155	310 1,792 0.29	148.5	297 1,890 0.99

Table 3 .

 3 The maximum send and send-receive volumes, and the total volume for PaToH and UMPa when K = 256. The times are given in seconds. There are 20 graphs in the table where two graphs with the best and the worst improvements on maxSV are selected from each class. Each number is the average of 10 different executions.

			PaToH			UMPa
	Graph	maxSV maxSRV	totV Time maxSV maxSRV	totV Time
	coPapersCiteseer	7,854	16,765 577,278 224.09	5,448	11,615 579,979 658.21
	coPapersDBLP	14,568	34,381 1,410,966 143.97 10,629	23,740 1,371,425 1038.86
	as-22july06	1,555	7,128	28,246 1.01	617	4,543	33,347 12.62
	smallworld	1,045	2,078 232,255 4.55	877	1,751 208,860 36.24
	delaunay n20	301	600	57,089 17.98	279	566	58,454 68.85
	delaunay n21	420	844	80,603 35.01	398	813	83,234 107.35
	hugetrace-00000	407	814	74,563 55.51	415	831	80,176 123.66
	hugetric-00010	502	1,004	91,318 92.45	477	955	97,263 167.69
	adaptive	753	1,505 143,856 96.60	735	1,472 152,859 224.30
	venturiLevel3	568	1,137 107,920 49.97	564	1,132 114,119 132.02
	rgg n 2 22 s0	799	1,589 145,902 151.30	724	1,495 147,331 249.23
	rgg n 2 23 s0	1,232	2,432 219,404 347.32	1,062	2,168 221,454 446.78
	ri2010	3206	5,989 281,638 0.72	2,777	5,782 279,941	8.66
	tx2010	5,139	9,230 124,033 8.47	3,011	7,534 117,960 15.55
	af shell10	898	1,792 174,624 89.90	885	1,769 184,330 158.04
	audikw1	4,318	8,299 680,590 322.57	3,865	7,607 692,714 822.73
	asia.osm	72	146	4,535 72.37	66	135	4,484 18.79
	great-britain.osm	104	209	11,829 50.52	82	168	11,797 25.51
	finan512	199	420	36,023 2.75	192	437	36,827 27.70
	memplus	1,860	7,982	15,785 0.49	946	4,318	19,945	8.25

Table 4 .

 4 The relative performance of UMPa with respect to Pa-ToH in terms of execution time. The numbers are computed by using the results of 10 executions for each of the 123 graphs in our data set.

	K	4	16	64 256 Avg.
	Relative time 1.02 1.29 2.01 5.76 1.98

This work was supported in parts by the DOE grant DE-FC02-06ER2775 and by the NSF grants CNS-0643969, OCI-0904809, and OCI-0904802. c 0000 (copyright holder) 1

(bestM axSV, bestM axSRV, bestT otV ) ← (maxSV, maxSRV, totV