Five poTAGEs and a COLT for an unrealistic predictor * Pierre Michaud Inria

Summary

The predictor we submit to the competition is for exploring the limits of branch prediction, it is not intended to be implemented in hardware. The proposed predictor consists of 5 huge TAGE predictors which we have tuned for better accuracy under huge storage budget. The overall prediction is obtained by combining the 5 TAGE predictions via COLT fusion. One of the 5 TAGE components takes a conventional global path as input. The 4 other TAGE components take as input different sorts of firstlevel histories. On the CBP-4 traces, the proposed predictor achieves 1.782 MPKI.

The TAGE predictor

The predictor we submit to the unlimited-size track is based on a modified version of the TAGE predictor [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF]. A TAGE predictor consists of several tagged tables and one tagless table. Each entry of the tagless table contains a taken/not-taken counter. Each entry of a tagged table contains a tag, a taken/not-taken counter, and a u bit. The tagless table is indexed with the address of (conditional) branches: it provides a default prediction when none of the tagged table give a hit. The tagged tables are indexed through a hash on the branch address and the global path [START_REF] Nair | Dynamic path-based branch correlation[END_REF]. The path lenghs follow a geometric progression [START_REF] Seznec | Analysis of the O-GEometric History Length branch predictor[END_REF]. TAGE is based on the PPM principle [START_REF] Chen | Analysis of branch prediction via data compression[END_REF]. In a PPM-like predictor [START_REF] Michaud | A PPM-like, tag-based predictor[END_REF], the prediction is always given by the longest hitting path. However, there are some situations where the second longest hitting path is likely to be more accurate than the longest one, especially when the longest hitting path entry is "fresh", i.e., it has just been initialized but not yet updated. TAGE exploits this fact with a single USE ALT ON NA counter (acting like a meta-predictor [START_REF] Mcfarling | Combining branch predictors[END_REF]) for selecting between a fresh longest hitting prediction and the second longest hitting prediction [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF]. The u bit in each entry indicates if that entry is likely to be useful or not. When the u bit is set, it pro- tects the entry from being "stolen" [START_REF] Michaud | A PPM-like, tag-based predictor[END_REF]. All the u bits are cleared when it is detected that the predictor is "jammed" [START_REF] Seznec | A new case for the TAGE branch predictor[END_REF]. The u bit is also used to detect if an entry is fresh: an entry is considered fresh if the u bit is not set and the taken/not-taken counter value is weak [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF].

3 The poTAGE predictor

Post predictor

The poTAGE predictor, shown in Figure 1, is similar to TAGE, except that we replace the USE ALT ON NA mechanism with a post-predictor (poTAGE = postpredicted TAGE). In the submitted predictor, the taken/not-taken counter in each TAGE entry is 3-bit wide.

Our post-predictor is a 1024-entry table, each table entry holding a 5-bit taken/not-taken counter. The postpredictor is indexed with the u bit and counter value of the longest hitting TAGE entry, and the counter values of the second and third longest hitting TAGE entries (10 index bits total). The 5-bit counter is used and updated like a conventional taken/not-taken counter [START_REF] Smith | A study of branch prediction strategies[END_REF].

Ramp-up

In a TAGE predictor, the update is crucial, it must be done very carefully: only the longest hitting counter is updated [START_REF] Michaud | A PPM-like, tag-based predictor[END_REF], and only a few new entries are allocated upon a misprediction for path lengths greater than the longest hitting length [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF]. However, because we assume a huge predictor size, it is possible to use a more aggressive update policy to decrease mispredictions due to cold-start effects: Instead of updating only the longest hitting counter, we update all the hitting counters, whether or not the branch was correctly predicted.

Instead of allocating a few new entries, and only upon a misprediction, we allocate entries for all the path lengths greater than the longest hitting length, whether or not the branch was correctly predicted, and provided the entries can be stolen (u bit not set). We stop doing aggressive allocation for path lengths greater than 200 branches when all the hitting counters are saturated.

We use this aggressive update policy during what we the ramp-up period. When the ramp-up period is over, we switch to the careful update policy implemented in the ISL-TAGE predictor [START_REF] Seznec | A new case for the TAGE branch predictor[END_REF].

The ramp-up period length is roughly proportional to the predictor size. For the submitted predictor, we set the ramp-up period to one million mispredictions 1 . A rampup period might also be useful in (large) real branch predictors. However the problem in a real predictor is how to detect when to start the ramp-up.

Beyond TAGE

So far, the most accurate branch predictors, including TAGE, belong to the family of two-level branch predictors [START_REF] Yeh | Two-level adaptive training branch prediction[END_REF][START_REF] Pan | Improving the accuracy of dynamic branch prediction using branch correlation[END_REF][START_REF] Yeh | Alternative implementations of two-level adaptive branch prediction[END_REF]. Yeh and Patt provided a taxonomy of two-level branch predictors [START_REF] Yeh | A comparison of dynamic branch predictors that use two levels of branch history[END_REF]. In particular, they distinguished predictors according to their first history level. The 2006 and 2011 TAGE predictors [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF] are "global" schemes, as their first-level history consists of a single global path shared by all the branches.

Global schemes are the most accurate two-level predictors that we know. However they have some limits. The footprint of a static branch in a particular TAGE table grows exponentially with the path length used for that table. The footprint also grows very quickly with branch entropy: the more unpredictable the branches, the larger the footprint [START_REF] Michaud | A comprehensive study of dynamic global history branch prediction[END_REF]. A large footprint incurs cold-start and capacity misses in a TAGE predictor. Applications exhibiting randomness in their control flow behavior are difficult 1 The careful update policy is almost never used in our simulations.

to predict by a global scheme such as TAGE not only because of their intrinsic randomness, but also because of the large footprint. Per-address schemes, i.e., predictors that use perbranch first-level histories [START_REF] Yeh | A comparison of dynamic branch predictors that use two levels of branch history[END_REF], suffer less from this problem. Indeed, a global path contains some "noise", i.e., some entropy that does not bring any correlation information but just grows the footprint [START_REF] Evers | An analysis of correlation and predictability: what makes two-level branch predictors work[END_REF]. In an ideal (i.e, large enough) per-address scheme, each static branch is predicted with its own subpath, which it does not share with other static branches. Unlike global schemes, per-address schemes cannot exploit global branch correlations, only local correlations. However, in a per-address scheme, hard-to-predict branches do not "pollute" the subpath of easy-to-predict branches.

Global and per-address schemes are two extreme points in the set of all possible two-level predictors. Yeh and Patt also introduced per-set schemes in their taxonomy [START_REF] Yeh | A comparison of dynamic branch predictors that use two levels of branch history[END_REF]. Per-set schemes have been largely ignored since their introduction, until recently when Ishii et al. introduced two per-set schemes in their FTL++ predictor [START_REF] Ishii | Revisiting local history for improving fused two-level branch prediction[END_REF] In a per-set scheme, there are several subpaths like in a per-address scheme, but a subpath is shared by several static branches. Subpath sharing generally exists in practical implementations of per-address schemes because of the limited storage for first-level histories. However, in a per-set scheme, subpath sharing is intentional: the goal is to exploit both local and (to some extent) global correlations, but with a smaller footprint than global schemes.

In a global scheme such as TAGE, hash functions are applied to the global path for indexing the tagged tables. Path aliasing can occur if two different global paths map to the same table entry and tag. In a global scheme, we generally try to minimize path aliasing by using "good" hashing functions. Using a per-set scheme can be seen as intentionally introducing path aliasing for decreasing the footprint. Blind path aliasing is likely to degrade the predictor's ability to exploit global correlation. Still, the hope is that the benefit of a smaller footprint outweighs this degradation. This is generally not true for all branches but for a few of them. This is why per-set schemes should not be used alone but as components of a hybrid predictor, as Ishii et al. did in the FTL++ predictor. The type of path aliasing that per-set schemes exploit is based on the assumption that branches which are statically close to each other are more likely to be correlated than distant branches. Other forms of intentional path aliasing can be imagined (cf. Section 6

The multi-poTAGE predictor

The proposed predictor, called multi-poTAGE, is depicted in Figure 2. It is a hybrid predictor [START_REF] Mcfarling | Combining branch predictors[END_REF] combining 5 different poTAGE predictors P0, P1, P2, P3 and P4. P0 is the global path poTAGE described in Section 3. The other predictors do not use a global path:

• P1 uses 32 per-address subpaths • P2 uses 16 per-set subpaths with 128-byte sets • P3 uses 4 per-set subpaths with 2-byte sets • P4 uses 8 frequency-based subpaths Predictor P4 uses a new sort of first-level history that we describe later in Section 6.

The 5 poTAGEs are combined using COLT fusion, a method invented by Gabriel Loh and Dana Henry [START_REF] Loh | Predicting conditional branches with fusion-based hybrid predictors[END_REF]. The COLT table is indexed with the branch address. Each COLT entry holds a 5-bit taken/not-taken counter.

Table 5 gives an overview of the 5 poTAGEs implemented in the submitted multi-poTAGE. In a two-level predictor, the first-level history consists of a set of subpaths that we call a path spectrum. At prediction time, a subpath S p is selected from the path spectrum. Subpath S p is used to access the second-level history (here, poTAGE), yielding a prediction. At update time, the branch updates one of the subpath S u . In the proposed multi-poTAGE, S p and S u are always the same.

The spectrum of P0 consists of a single subpath shared by all the branches, i.e., it is the whole global path. The spectrum of P1 consists of 32 subpaths, and S p (= S u ) is selected with the 5 least significant branch address bits 0 to 4 (per-address first level). The spectrum of P2 consists of 16 subpaths, and S p (= S u ) is selected with branch address bits 7 to 10 (per-set first level).

The spectrum of P4 consists of 8 subpaths. The subpath S p (= S u ) is selected as follows. The branch address is used to index a Branch Frequency Table (BFT). If the BFT is large enough, each static branch uses a distinct BFT entry. Each BFT entry holds a counter indicating the current frequency of the static branch. The frequency of a branch is the number of times the branch has been executed until now since the counter was reset 2 .

Predictor P4 seeks to exploit correlations between branches having (roughly) the same frequency. Each of the 8 subpaths S[0] to S [START_REF] Michaud | A comprehensive study of dynamic global history branch prediction[END_REF] corresponds to a distinct frequency bin. Let F max be the maximum branch frequency so far. Branches whose frequency lies in [F max /2, F max ] are predicted with subpath S[0]. Branches whose frequency lies in [F max /4, F max /2[ use subpath S [START_REF] Chen | Analysis of branch prediction via data compression[END_REF]. Branches whose frequency lies in [F max /8, F max /4[ use subpath S [START_REF] Evers | An analysis of correlation and predictability: what makes two-level branch predictors work[END_REF]. And so forth.

Note that, after an initial period, the dynamic instances of a static branch will generally use the same subpath.

7 Experimental analysis P0 is the most accurate of the 5 poTAGEs as a single component. Compared with the USE ALT ON NA mechanism used in previous TAGE predictors, the post-predictor decreases the number of mispredictions by 5% on the CBP-4 traces. Using a ramp-up period yields a further 4% decrease. To evaluate the importance of each of the 5 poTAGEs in the multi-poTAGE predictor, we selectively disabled some of them. P2 is the second most important predictor after P0, followed by P1, P4 and P3 in that order. Adding P2, P1, P4, and P3 (successively and in that order) to P0 decreases the number of mispredictions by respectively 5%, 3%, 2.5% and 1%. In total, the multi-poTAGE predictor has 10% fewer mispredictions than P0 alone.

Conclusion

The prediction accuracy of TAGE is very high. Still, when considering a huge storage budget for limit studies, a more aggressive update policy can be used. Global path predictors such as TAGE are theoretically able to capture all the branch correlations if the global path is long enough. However, a global path suffers from "noise" from branches that bring no correlation information but grow the footprint. To (try to) solve this issue, we combined five poTAGEs with different first-level histories, using COLT fusion. This brute force approach may not be reasonable for a realistic predictor. Future research should try to look for more cost-effective solutions.

Figure 2 :

 2 Figure 2: The multi-poTAGE predictor.

6

  Frequency-based path spectrum P4 uses a new sort of first-level history: frequency-based subpaths. P4 does not belong to the taxonomy of Yeh and Patt but to a more general class of two-level predictors.

Table 1 :

 1 ). The 5 poTAGEs. The poTAGE tables are very large. P0 is roughly twice larger than each of the other predictors. The COLT table (not listed) is also very large.

	poTAGE	spectrum	spectrum size # tagged min subpath max subpath
		type	(# subpaths)	tables	(# branches) (# branches)
	P0	global	1	20	7	5000
	P1	per address	32	19	5	2000
	P2	per 128-byte set	16	19	5	500
	P3	per 2-byte set	4	19	5	500
	P4	frequency	8	19	5	500

In the submitted branch prediction algorithm, we reset the frequency counters only once, when the simulation starts.


Pierre Michaud

André Seznec

Pushing the branch predictability limits with the multi-poTAGE+SC predictor *

Introduction

The multi-poTAGE+SC predictor (Figure 1), which we submit to the unlimited size track, combines the multi-poTAGE and TAGE-SC predictors that have been submitted separately to CBP-4 by the first and second authors respectively [START_REF] Michaud | Five poTAGEs and a COLT for an unrealistic predictor[END_REF][START_REF] Seznec | Exploring the branch prediction accuracy limits for the TAGE-SC branch predictor[END_REF]. Both multi-poTAGE and TAGE-SC are based on the TAGE predictor. We combine them by replacing the TAGE component in TAGE-SC with the multi-poTAGE. On the CBP-4 traces, the proposed predictor achieves 1.691 mispredictions per thousand instructions (MPKI).

Sections 2 and 3 of this paper describe briefly the multi-poTAGE predictor. Section 4 describes the statistical corrector. Section 5 provides some performance analysis. [START_REF] Gao | PMPM: Prediction by combining multiple partial matches[END_REF] The poTAGE predictor TAGE is derived from branch prediction by partial matching (PPM) [START_REF] Chen | Analysis of branch prediction via data compression[END_REF][START_REF] Michaud | A PPM-like, tag-based predictor[END_REF]. The poTAGE predictor, shown in Figure 2, is similar to TAGE [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF], except that we replace the USE ALT ON NA mechanism with a post-predictor (poTAGE = post-predicted TAGE). In the submitted predictor, the taken/not-taken counter in each TAGE entry is 3-bit wide. Our post-predictor is a 1024-entry table, each table entry holding a 5-bit taken/not-taken counter. The post-predictor is indexed with the u bit and counter value of the longest hitting TAGE entry, and the counter values of the second and third longest hitting TAGE entries (10 index bits total). The 5-bit counter is used and updated like a conventional taken/not-taken counter [START_REF] Smith | A study of branch prediction strategies[END_REF].

In a TAGE predictor, the update is crucial, it must be done very carefully: only the longest hitting counter is updated, and only a few new entries are allocated upon a misprediction for path lengths greater than the longest hitting length [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF]. However, because we assume a huge predictor size, it is possible to use a more aggressive update policy to decrease mispredictions due to cold-start effects:

Statistical

• Instead of updating only the longest hitting counter, we update all the hitting counters, whether or not the branch was correctly predicted.

• Instead of allocating a few new entries, and only upon a misprediction, we allocate entries for all the path lengths greater than the longest hitting length, whether or not the branch was correctly predicted, and provided the entries can be stolen (u bit not set). We stop doing aggressive allocation for path lengths greater than 200 branches when all the hitting counters are saturated.

We use this aggressive update policy during what we call the ramp-up period. When the ramp-up period is over, we switch to the careful update policy implemented in the ISL-TAGE predictor [START_REF] Seznec | A new case for the TAGE branch predictor[END_REF].

The ramp-up period length is roughly proportional to the predictor size. For the submitted predictor, we set the ramp-up period to one million mispredictions.

The multi-poTAGE predictor

The multi-poTAGE predictor, depicted in Figure 3, consists of 5 different poTAGE predictors, P0, P1, P2, P3 and P4, combined via COLT fusion [START_REF] Loh | Predicting conditional branches with fusion-based hybrid predictors[END_REF]. P0 is the global path poTAGE described in Section 2. The other poTAGEs do not use a global path:

• P1 uses 32 per-address subpaths

• P2 uses 16 per-set subpaths with 128-byte sets

• P3 uses 4 per-set subpaths with 2-byte sets

• P4 uses 8 frequency-based subpaths P2 and P3 are two per-set two-level predictors [START_REF] Yeh | A comparison of dynamic branch predictors that use two levels of branch history[END_REF]. Ishii et al. also used per-set schemes in their recent FTL++ predictor [START_REF] Ishii | Revisiting local history for improving fused two-level branch prediction[END_REF]. P4 uses a new sort of first-level history, frequencybased subpaths. In a two-level predictor, the first-level history consists of a set of subpath that we call a path spectrum. At prediction time, a subpath is selected from the path spectrum. The selected subpath is used to access the second-level history (here, poTAGE), yielding a prediction. The spectrum of P4 consists of 8 subpaths. The subpath is selected as follows. The branch address is used to index a Branch Frequency Table (BFT). If the BFT is large enough, each static branch uses a distinct BFT entry. Each BFT entry holds a counter indicating the current frequency of the static branch. The frequency of a branch is the number of times the branch has been executed until now since the counter was reset 1 .

Predictor P4 seeks to exploit correlations between branches having (roughly) the same frequency. Each of Branches whose frequency lies in [F max /8, F max /4] use subpath S [START_REF] Gao | PMPM: Prediction by combining multiple partial matches[END_REF]. And so forth.

The multi-poTAGE paper provides a more detailed description [START_REF] Michaud | Five poTAGEs and a COLT for an unrealistic predictor[END_REF].

The Statistical Corrector

The Statistical Corrector predictor was introduced as an adjunct predictor to TAGE in [START_REF] Seznec | A 64 kbytes ISL-TAGE branch predictor[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF]. It allows to better predict the class of branches that are not completely determined by the path, but are globally biased by path. The correction aims at detecting the unlikely predictions from TAGE and to revert them: the prediction flowing from TAGE, its confidence level, as well as information on the branch (address, global history, global path, local history) are presented to the Statistical Corrector predictor which decides whether or not to invert the prediction.

In this submission, we use the output of the multi-poTAGE predictor and the prediction outputs of its internal components as inputs of a Statistical Corrector.

We use a statistical corrector similar to the one used in the companion submission [START_REF] Seznec | Exploring the branch prediction accuracy limits for the TAGE-SC branch predictor[END_REF]. It is a perceptroninspired [START_REF] Jiménez | Neural methods for dynamic branch prediction[END_REF] Statistical Corrector [START_REF] Seznec | A 64 kbytes ISL-TAGE branch predictor[END_REF][START_REF] Seznec | A new case for the TAGE branch predictor[END_REF], that combines multiple components:

• the Bias table: indexed through the PC and the multi-poTAGE predicted direction (after selection through the COLT component).

• A second Bias table indexed with the PC, the multi-poTAGE output and the 5 poTAGE component outputs.

• 4 LGEHL components, 3 using a 16-entry history table, and one using 64K-entry history table. Each one features 15 tables.

• 4 perceptron-derived local history components using similar history tables. In these perceptron-derived components, we use the MAC representation of the counters [START_REF] Seznec | Revisiting the perceptron predictor[END_REF]; a counter is associated with 6 consecutive bits of history. Each of these components features 10 tables.

• 2 perceptron-derived components using respectively global branch history and global path history: 10 tables each.

• a global history GEHL component: 209 tables

• a global history component inspired from the MAC-RHSP predictor [START_REF] Seznec | Revisiting the perceptron predictor[END_REF]; a counter is associated with 6 consecutive bits of history and part of the global branch history (1/3) is hashed with the PC: 80 tables.

• Two path skeleton history GEHL components. The first path skeleton are the taken branches whose targets are not too close to the branch source. By too close, we mean 16 bytes for backward branches and 128 bytes for forward branches. The second path skeleton registers the branch in the path only if it was not among the last 8 encountered branches. These components feature 15 tables each.

• Two path skeleton history perceptron-derived components: 10 tables each.

All the tables hold 8 bit counters. The prediction is computed as the sign of the sum of the (centered) predictions read on all the Statistical Corrector tables: a total of 461 counters are summed. The Statistical Corrector predictor tables are updated using a dynamic threshold policy as suggested for the GEHL predictor [START_REF] Seznec | Analysis of the O-GEHL branch predictor[END_REF]. As suggested in [START_REF] Daniel | OH-SNAP: Optimized hybrid scaled neural analog predictor[END_REF], we use a PC-indexed table of dynamic threshold, which yields marginal benefit.

Any of these components has only a limited accuracy impact, but if one removes all the components exploiting local history or exploiting global branch/path history the impact on accuracy is more significant.

Choosing between multi-poTAGE and the statistical corrector ouputs

The prediction flowing out of the statistical corrector is often more accurate than the poTAGE prediction. However using a chooser results in a higher accuracy than just using the statistical corrector output.

The best chooser we have experimented uses a GEHL + LGEHL structure with limited number of tables (i.e. short histories). It is indexed with the PC, the multi-poTAGE prediction, the confidence (high or not ) of the SC prediction.

Performance analysis

The multi-poTAGE and statistical corrector configurations used here are very close to those described in the companion papers [START_REF] Michaud | Five poTAGEs and a COLT for an unrealistic predictor[END_REF][START_REF] Seznec | Exploring the branch prediction accuracy limits for the TAGE-SC branch predictor[END_REF], where they were analyzed separately 2 . Table 1 provides some performance analysis about the combined multi-poTAGE+SC predictor. We disabled some features of multi-poTAGE+SC to see how they impact the overall performance.

When we apply the statistical corrector on the multi-poTAGE predictor, more than half of the accuracy gain comes from the bias and global-history components of the statistical corrector (cf. Table 1, compare MP, MP+SC and MP+SCg).

The local-history statistical corrector components have more effect when statistical correction is applied on a single global-path poTAGE: P0+SCg has 6.3% (12.2% vs. 5.6%) more mispredictions than P0+SC. That is, the nonglobal components of the multi-poTAGE (P1,P2,P3,P4) and the non-global components of the statistical corrector have partially redundant effects. Still, the non-global components of multi-poTAGE bring a 5.3% prediction accuracy gain (MP+SC vs. P0+SC).

It is interesting to notice that when statistical correction is applied to P0, removing the post-predictor from P0 has little impact: P0+SCg and PPM+SCg have close MPKIs (respectively +12.2% and +13.7% more mispredictions than MP+SC). However, removing the postpredictor from an isolated poTAGE degrades prediction accuracy significantly: PPM has about 10% more mispredictions than P0 (+30.8% vs. +18.7%).

In a PPM-like predictor, the longest matching context is not always the most accurate [START_REF] Seznec | A case for (partially) tagged geometric history length branch prediction[END_REF][START_REF] Gao | PMPM: Prediction by combining multiple partial matches[END_REF]. It is this problem that the post-predictor and the statistical corrector try to solve, but the statistical corrector is a more effective solution.

It is interesting to notice that statistical correction brings the same relative improvement over a single poTAGE (P0+SCg vs. P0) than over the multi-poTAGE (MP+SC vs. MP).

Conclusion

On the CBP-4 traces, the multi-poTAGE+SC predictor achieves about 5 % fewer mispredictions than the multi-poTAGE predictor or the TAGE-SC predictor that are presented in two other submissions by respectively the first and the second author.

The non-global components of the multi-poTAGE are responsible for most of the prediction accuracy gain of multi-poTAGE+SC over TAGE-SC. The effectiveness of the statistical corrector lies in the fact that, in a PPM-like predictor such as TAGE or poTAGE, the longest matching context is not always the most accurate.

Figure 1 :

 1 Figure 1: The multi-poTAGE+SC predictor.

Figure 2 :

 2 Figure2: The poTAGE predictor. The post-predictor takes as input the u bit and taken/not-taken counter of the longest hitting TAGE entry, and the taken/nottaken counters of the second and third longest hitting TAGE entries.

Figure 3 :

 3 Figure 3: The multi-poTAGE predictor.

Table 1 :

 1 Average MPKI increase for various configurations. MP is the multi-poTAGE (5 poTAGEs + COLT); SC is the statistical corrector; SCg is the statistical corrector reduced to its bias and global-history components; P0 is the global-path poTAGE alone; PPM is P0 without the post-predictor. The base MPKI for all configurations is that of MP+SC.

	configuration MPKI increase
	MP+SC (base)	0
	PPM	+30.8%
	P0	+18.7%
	MP	+5.7%
	P0+SC	+5.6%
	MP+SCg	+2.5%
	P0+SCg	+12.2%
	PPM+SCg	+13.7%

In the submitted branch prediction algorithm, we reset the frequency counters only once, when the simulation starts.

The configuration parameters used here are slightly different from those used in[START_REF] Michaud | Five poTAGEs and a COLT for an unrealistic predictor[END_REF][START_REF] Seznec | Exploring the branch prediction accuracy limits for the TAGE-SC branch predictor[END_REF] because of the memory size constraint.

* This work was partially supported by the European Research Council Advanced Grant DAL No 267175


J Briottet

L Pierro
email: pierro@onera.fr

Labarre

N Lalanne

J-C Krapez
email: krapez@

C Le Niliot

X Briottet
email: xavier.briottet@onera.fr

J Pierro

L Labarre
email: luc.labarre@onera.fr

Development and validation of a numerical tool for the simulation of the temperature field and infrared radiance rendering in an urban scene

come

Development and validation of a numerical tool for the simulation of the temperature field and infrared radiance rendering in an urban scene

Introduction

Human activity in the urban context influences the local environment. For instance a practical consequence in terms of urban climatology is the development of urban heat islands. The urban climate is the result of thermal and water flow processes between the urban surfaces and the environment. Infrared remote sensing is a powerful tool for measuring the influence of these interactions on the urban surfaces at different spatial scales. The infrared images obtained with airborne sensor inform about the thermal balance of urban components and, at a large scale, about the urban heat island. The radiance images are however a complex result of various influences both on the thermal level and on the radiative level; therefore a theoretical model is often necessary for interpreting these images. This model should first take into account the energy budgets for the different surfaces (roofs, roads, walls, etc.) including the radiative interactions between them. It should then solve the heat transfer problem through the solid materials of the urban scene at a scale compatible with the sensor resolution, which can be of the order of a tenth of a meter. The computed surface temperature field has finally to be combined with the emissivity field in a radiative transfer module for getting the atsensor infrared radiance image. The energy exchanges between cities and the atmosphere can be simulated at a scale between the mesoscale (whole city) and the building scale; such models as for example TEB [START_REF] Masson | A physically-based scheme for the urban energy balance in atmospheric models[END_REF] and SM2-U [START_REF] De La | Modelling microclimate in urban environments and assessing its influence on the performance of surrounding buildings[END_REF]. A second category of models includes the thermo-radiative models establishing the energy budget at a scale significantly lower than the building size (typically one meter or less). These models estimate all terms of the energy budget, in order to compute the temperature. A list of such models is described in table 1.

Table 1. Main existing thermo-radiative urban models

Model : autors, (NAME), year Groleau et al. SOLENE [START_REF] Groleau | Simulation of the radiative behavior of an urban quarter of Marseille with the SOLENE model[END_REF] (2003) Krayenhoff, Voogt TUF-3D [START_REF] Krayenhoff | A microscale three-dimensional urban energy balance model for studying surface temperatures[END_REF] (2007) Asawa, Hoyano [START_REF] Asawa | Thermal design tool outdoor spaces based on heat balance simulation using a 3D-CAD system[END_REF] (2008) Kastendeuch, Najjar [START_REF] Kastendeuch | Simulation and validation of radiative transfers in urbanized areas[END_REF] (2008)

Yang, Li MUST [START_REF] Yang | Development of a 3D urban energy model for predicting and understanding surface temperature distribution[END_REF] The thermo-radiative model SOLENE [START_REF] Groleau | Simulation of the radiative behavior of an urban quarter of Marseille with the SOLENE model[END_REF], developed at CERMA Nantes, solves the heat conduction problem in walls with a nodal method based on a preliminary evaluation of the radiative heat fluxes on the boundaries. This 5 nodes model can describe a wall with only one or two layers. The simplicity of this model does not correctly capture the behavior of complex walls, heavy inertial materials, or semi-infinite ground [START_REF] Hénon | High resolution thermo-radiative modeling of an urban fragment in Marseilles city center during the UBL-ESCOMPTE campaign[END_REF]. The radiative exchanges between the surfaces and the sky are computed from the view factors and an approximate radiosity method, partly based on an iterative algorithm named "progressive refinement". In the IR domain, a blackbody assumption is used for computing the radiative heat exchanges between surfaces. Experimental validations have shown that this model can induce deviations between the simulated and the measured heat fluxes, and thereafter between the simulated and measured temperatures [START_REF] Idczak | An application of the thermo-radiative model SOLENE for the evaluation of street canyon energy balance[END_REF]. [START_REF] Wang | Transient heat flow calculation for multilayer constructions using a frequency-domain regression method[END_REF] identify three general classes of methods for solving the heat conduction equation: numerical methods, response factor (RF) methods and conduction transfer function (CTF) methods [START_REF] Stephenson | Calculation of heat conduction transfer functions for multilayer slabs[END_REF]. Numerical methods correspond to methods discretizing the studied domain in meshes, and where the solutions of the differential equations are approximated by numerical techniques as FDM, FVM, FEM, etc. The RF method is based on the responses of the thermal system to an elementary stress and they are convoluted by the observed thermal excitations. The CTF takes further into account previous values of the searched variable. These methods are based on the Laplace transformation and are commonly used in building energy software packages. For instance, the RF method is adopted by DOE-2 whereas the CTF method is found in EnergyPlus or TRNSYS. Both methods are used for computing the wall heat flux density from the outside and inside temperatures (actually the outside sol-air temperature and the inside air temperature) which actually can be considered as an inverse problem. However the main unknown variable for simulating an infrared image of a given urban scene is not the wall heat flux density but the building external temperature, therefore we developed our own thermal model called SUSHI (Simulation in Urban Scene of Heat dIffusion).

It is linked to SOLENE from which it gets some radiative parameters for evaluating the radiative energy budget like the sky model, the view factor matrix and the net radiative flux in the solar spectrum. SUSHI solves the heat equation for wall elements where the heat conduction is either 1D or 2D. The purpose is to simulate extended areas where heat flows normally to the surface and localized areas, generally associated to thermal bridges, where the heat flow is 2D, as for L-shapes (building corners) and T-shapes (wall/floor junction).

Once the temperature field is known, the radiance images at sensor level can be realized. A radiative transfer code is required to take into account the photonic interactions in the scene. The atmospheric radiative transfer codes are based on the radiative transfer equation and can calculate spectral characteristics of the atmosphere along the line of sight (transmission, self-emission, molecular scattering …). They are limited to simplified geometry (plane or spherical ground). 3D radiative transfer codes were thus developed in order to take into account radiative interactions between 3D elements of the scene. In the infrared domain, such codes are DART [START_REF] Gastellu-Etchegorry | A 3d model for remote sensing images and radiative budget of earth surfaces[END_REF] and TITAN [START_REF] Fontanilles | Thermal infrared radiance simulation with aggregation modeling TITAN an infrared radiative transfer model for heterogeneous 3D surface -application over urban areas[END_REF]. Main limits DART are in the spatial resolution of its thermal model. On the other hand TITAN computes radiative transfers with a high spatial resolution, but it has no thermal model. Thereafter, we use the new radiative transfer code MOHICANS which is a combination of AMARTIS [START_REF] Thomas | AMARTISv2: 3D radiative transfer code in the [0.4; 2.5 µm] spectral domain dedicated to urban areas[END_REF] and TITAN. The interest of this code is to simulate a hyperspectral sensor in [0.4 -14 µm].

In the following we will describe the main aspects of the thermal code SUSHI and its coupling with SOLENE and MOHICANS. Then we will present the results of various validation steps, in particular the comparison with experimental results obtained during the BATIR 2013 measurement campaign where infrared thermography of two buildings was performed with two ground cameras and an airborne camera.

Models

Outline

In order to obtain a simulated spectral image of an urban scene in the thermal IR spectral domain [3-15 µm], the temperature field over all external surfaces has to be computed. For this purpose, the heat equation describing the thermal transfer through the walls, roofs and soil has to be solved. The high spatial resolution of present IR cameras requires considering surface elements of area less than one square meter. Depending on their location on the wall, the heat transfer through these elements can be considered as one dimensional or not. Generally speaking, 2D transfer is observed at each corner of the building envelope (L-shape), and at thermal bridges like the floor/wall junctions (T-shape). 3D transfer is observed for example at a corner of a floor/wall junction. Finally, 1D transfer is observed over the building envelope far from previous thermal singularities. The numerical tool we developed allows simulating the thermal behaviour of walls presenting areas where heat flow is one of these two types: 1D or 2D. It takes into account the atmospheric conditions as measured with a series of meteorological sensors. It also takes into account the building internal conditions as recorded with temperature sensors. These external/internal conditions are recorded for a sufficiently long period of time for allowing a dynamic thermal analysis prone to capture all thermal inertia effects.

A series of programs are used as described in the chart in figure 1: -SketchUp is used to build-up the 3D-scene and subdivide each wall in areas where heat flow is either 1D or 2D; -SOLENE software which was developed at CERMA laboratory is used for meshing each area, for computing the view factors and then the solar radiation energy budget after multiple reflections modelling and shadow tracing at each time step; -SUSHI thermal software (Simulation in Urban Scene of Heat dIffusion) is used for modeling the infrared radiation energy budget, the heat convection flux and then for solving the 1D and 2D heat transfer problems inside the wall and soil elements.

-MOHICANS is finally used to get the spectral infrared radiance reaching each element of the sensor matrix from the temperature field computed by SUSHI and from the atmosphere radiative parameters computed by MATISSE [START_REF] Labarre | An overview of MATISSE-v2.0[END_REF]. Temperature and fluxes are assumed to be periodic. This allows applying the Fourier transform which, thanks to the Fast Fourier Transform tool, leads to semi-analytical solutions and yields the results in a very short time. However, for the input, we have to consider real meteorological data and these data are in part stochastic. The question was then how to choose the period aimed at representing the periodical phenomena. Of course it should be 24h at least. If it is too short, let us say only 24h, the computed temperature would be in error due to fact that the atmospheric conditions were actually different during any day before. One should thus extend the duration of the period. In this way, the first mismatch between the modelled atmospheric conditions (i.e. through the periodicity assumption) and the real atmospheric conditions would be rejected further back in the past. The error on the computed temperature would then be lowered at the expense of a longer computing time. The optimal choice for the period duration will be discussed in § 3.1.

Boundary conditions

Elements of mesh can exchange energetic fluxes with the environment. In the general case, these fluxes are convective and radiative. The radiative fluxes are subdivided into solar and infrared components. One supposes that solar radiation can be absorbed or reflected by a surface, and that infrared components can be absorbed, reflected and reradiated to others surfaces. The solar radiation may also be split into a direct and a diffuse component, in order to give hourly shadow distribution. The software SOLENE is used to compute these solar components. In this software, the direct component follows the Perrin de Brinchambault formula, and the atmospheric diffuse component follows the "all weather model for sky luminance" of Perez et al. [START_REF] Perez | All weather model for sky luminance distribution_Preliminary configuration and validation[END_REF]. This provides spectral-integrated fluxes depending on place and hour, and on sky clearness and brightness for the diffuse component. The energy budget through multiple reflections is computed by the progressive refinement radiosity method [START_REF] Cohen | A progressive refinement approach to fast radiosity image generation[END_REF] which first requires the calculation of the view factors between all surfaces of the scene including the "sky facets". SOLENE however doesn't apply the progressive refinement radiosity method for evaluating the infrared radiative net fluxes, it rather considers a blackbody assumption. For avoiding this approximation we chose to perform a more rigorous infrared radiative exchange calculation by implementing the same progressive refinement radiosity method as in the visible.

Let us now consider the convective flux on the outdoor building surface. It may be modelled through a convective coefficient which depends on wind velocity and orientation, on air temperature, on local and global geometry, etc. Using a 3D fluid mechanics solver provides a precise modelling for this coefficient [START_REF] Bouyer | Modélisation et simulation des microclimats urbains Étude de l'impact de l'aménagement urbain sur les consommations énergétiques des bâtiments[END_REF]. However, for keeping a reasonable computational time, we preferred to use one among several empirical formulas giving the convective coefficient vs. wind velocity in the building thermal energy context; the coefficient is thus time dependant:   t h . The thermal balance at the outdoor building surface can be expressed as follows:

      i j i net IR net vis a i i T T t R t R t T T t h z T          , , ) ( ) (  (1)
It can be considered as a second order boundary condition where the specified heat flux corresponds to the right-hand side of the equation. would also require iterations even though an alternative exists which avoids performing iterations but requires more involved algebra [START_REF] Krapez | Modélisation analytique de la variabilité spatiale et temporelle de la température de sols nus naturels, Congrès SFT 2004 Presqu'île de Giens[END_REF].

For speeding up the iteration process, the second order boundary condition can be replaced by the following third order boundary condition with a constant arbitrary coefficient fict h :

    t T T h z T sa i fict i        (2)
where  

              i j i net IR net vis a i fict i sa T T t R t R t T T t h h T t T       , , 1 (3)
In other words, the thermal excitation associated with the 3 rd order BC is:

            i j i net IR net vis a i fict T T t R t R t T t h T t h h      , , (4)
A first iterative process identified by (1) in figure 2 is implemented by updating only the convective flux term

    i fict T t h h 
. When the convergence criterion is reached, the radiative flux term

  i j i net IR T T t R  , ,
is updated (step (2) in figure 2) and the global iterative process is repeated. The iterative process is split into two steps because updating the radiative flux terms is the most expensive. For the indoor boundary condition, it's assumed that the environment is homogeneous and characterised by the air temperature and by a global constant radiato-convective transfer coefficient.

Fig. 2. Flow diagram of the computed terms by SOLENE and SUSHI

1D configuration

As explained before, the 1D transfer model is used for those facets of wall that are far from thermal singularities like facade corners, wall/floor junctions, … . Its semi-infinite version is used for soils. The walls may be composed of many layers. It is assumed that each layer is homogeneous, isotropic with constant thermal properties. The thermal problem is solved by applying the quadrupole technique after performing a Fourier transform to the heat equation and to the boundary conditions [START_REF] Maillet | Thermal quadrupoles: solving the heat equation through integral transforms[END_REF]. The following backward transfer relation relates the temperatures and conduction heat flux at front and back sides of the wall:

                         s s e e T D C B Α T   ~ (5)
where the tilde represents the Fourier transform of the corresponding variable, and where

                   layer s N i i i i i D C B A D C B A 1 (6)
is the total quadrupole matrix, which is calculated by multiplication of each quadrupole matrix, with A i , B i , C i and D i the quadrupole coefficients of the layer i. By isolating the front temperature it can be expressed by the indoor and outdoor exchanged net flux densities. The formulation of the terms of the quadrupole matrix proposed by [START_REF] Pailhes | Thermal quadrupole method with internal heat sources[END_REF] is used. The temperature expression is also written:

in ex ex Q Z Q X T   (7) with:   B h A X g in      1
the outdoor transfer function,

) ( ) ( ~air ex ex k ex ex k net IR abs VNIR ex T h TF T h TF Q        the outdoor net flux density,                         j Z layer s N i 1 i 1 exp
the indoor transfer function,

air in g in in T h Q   the indoor net flux density and D h C B h h A g in g in fict ex fict ex          , A

2D configuration

When the one-dimensional assumption is no valid, as observed in the vicinity of thermal bridges, the transfer functions related to those facets are pre-calculated by the finite element method as provided by the CAST3M software. The objective is to get periodic transfer functions corresponding to the response of a thermal bridge when an elementary flux is applied over each one of the characteristic surfaces of the thermal bridge. Thus, when considering the symmetric L-shape at left of figure 3, a uniform flux is applied on S1, then on S3 (due to symmetry, the response is the same when the perturbation is applied on S1 or on S2). Similarly, for the symmetric T-shape at right of figure 3 a uniform flux is applied on S1, then on S2. In all cases, a unitary step heating is applied for a relatively long period of time, typically several days. The transient response of the system is subtracted to itself after moving it one time step for getting the response to an elementary pulse. This response is then added to itself after moving it 24h, then 48h, … N times for getting the periodic response to an elementary pulse (N should be sufficiently long for the true periodic response). The Fourier transform of the periodic responses

   3 1 ) ( ). , ( ) , ( ~l k l k p l k Q x G x T    (8)
where ) ( ~k l Q  is the Fourier transform of the net heat flux density averaged over surface l S and

x is the position on the external surface with respect to the thermal bridge center. For returning into the temporal space an inverse fast Fourier transform is finally performed with the IFFTW algorithm [20]. The software MOHICANS is a numerical tool developed at ONERA which aims at evaluating all radiative components contributing to the optical signal received by a multi or hyperspectral sensor over a 3D scene. It is actually a fusion of previously developed codes AMARTIS [START_REF] Thomas | AMARTISv2: 3D radiative transfer code in the [0.4; 2.5 µm] spectral domain dedicated to urban areas[END_REF] for the visible to near infrared range and TITAN [START_REF] Fontanilles | Thermal infrared radiance simulation with aggregation modeling TITAN an infrared radiative transfer model for heterogeneous 3D surface -application over urban areas[END_REF] for the midwave to longwave infrared range. Both of them rely on atmospheric radiative computations by MATISSE [START_REF] Labarre | An overview of MATISSE-v2.0[END_REF]. The spectral domain finally covers the 0.4 to 15 µm range. MOHICANS is organized in four calculation steps: (i) computation of the sky irradiance hemisphere and the atmosphere transmission terms by calls to MATISSE, (ii) computation of geometrical terms corresponding to sun visibility, facet to facet visibility and facet to pixel visibility, (iii) computation of total irradiance, corresponding to the sum of solar direct irradiance, the atmosphere irradiance and the environment irradiance, (iv) computation of the radiance of the radiative flux reaching each pixel of the imaging sensor. A schematic of the radiative exchanges between two facets of a wall-soil angle is presented in figure 4.

Fig. 3. Drawing of L-shape (at left) and T-shape (at right) (imposed flow located by arrows)

Spectral model: MOHICANS

Fig. 4. Drawing of radiative terms computed by MOHICANS which contribute to the radiance at sensor level

Computation results

Temporal configuration

The first analysis concerns the validation of the time periodic model. For simulating the temperature field over a one day time interval [t 0 -24,t 0 ] we assume that the radiative and convective fluxes are periodic with a period D and that their variations are those observed during the time interval [t 0 -D,t 0 ] with D=kx24; k≥1. The longer the period D, the lower the error on temperature over the interval [t 0 -24,t 0 ] due to the wrong assumption of periodic phenomena. The minimum number k of days to consider for getting an error lower than a specified value was determined by considering both an elementary part of a typical wall (see properties in table 2) and a semi-infinite ground (effusivity: b=707 SI). Both were submitted to the real atmospheric fluxes observed during the BATIR campaign for a horizontal surface. A long period D of 12 days was considered as a reference. Temperature was computed for progressively shorter periods D and compared to this reference. The maximum difference and the RMS error value over the day of interest (i.e. the last one) are plotted in figure 3 for the wall and the ground). These simulations highlight that a period of two days is long enough for the wall (maximum error of 0.1°C; RMS<0.02°C), but that a period of five days is necessary for soils, when setting the maximum error criterion to 0.1°C. Anyway, with a period of two days, the maximum error on temperature for the soil is not more than about 0.2°C. For the remaining we will thus perform the thermal simulations with D=48h.

Fig. 5. RMS and absolute maximal difference temperature according to the number of considered periods for a wall (left)

and a ground (right)

Sensitivity study

A MATLAB program which calculates surface temperature in the very same manner as the 1D SUSHI model is used for local and fast simulations. Analytical validations of the quadrupole model were performed on basic configurations like a semi-infinite layer with or without convective losses, with a sinusoidal imposed boundary flux. Software validations on finite walls were then realized by means of comparison with SOLENE simulations. In the case of a single horizontal two layer slab (with an isolating layer), both models give similar results. However, SUSHI results are much closer to the analytical solution as compared to SOLENE results for very thick materials and for soils. In addition, SUSHI permits to model walls with more than two materials.

In a second phase, the sensitivities of the front temperature with respect to thermal, convective and radiative parameters are analysed. area series of 14 parameters were analyzed: the effusivity and the square root of the diffusion time each of the three layers of the wall (see table 2), the albedo and the emissivity of the external surface, the external and internal convective coefficients, the external and internal air temperature, and the VNIR (Visible and Near Infrared, i.e. solar spectrum) and IR incident flux density.

We define two groups of parameters, according to their time dependency. The external convective coefficient, the air external and internal temperature and the incident flux densities are those which vary with time. For the static parameters, as usual, we introduced a specific variation around their reference value:

  K ref   1  
. For the timedependent parameters, we performed two analyses. The first one was according to     

K t t ref   1   (labeled as
Wall Ground "systematic" in figure 7). For the second one (labeled as "stochastic" in figure 7), we added random variations to the nominal time evolution as observed over a period of two days. We applied a uniform centered distribution according to:

      1 2 1 , 0    u K t t ref ref    (9) with cst t t K ref    ) ( ) (  
and u 0,1 is the uniform distribution on the support [0,1] ;

) (t ref 
in Kelvin if  is the air temperature. In both cases of systematic deviation of the parameter, the reduced sensibility is approximated by the finite difference:

  K T T S i i i i        (10)
For the second analysis of the time-dependent parameters, series of simulations have to be performed with different random drawings. At each time, the temperature scatter is then compared with the nominal temperature at this time. The reduced sensitivity is then given by the ratio of the standard deviation of the simulated temperature with respect to the reference temperature, over the standard deviation of the uniform distribution of the considered parameter:

    K t t S T i 3    (11)
We plotted the reduced sensitivity to the thermal parameters (effusivity and square root of diffusion time) of the three layers of the wall in figure 5. For instance, a reduced sensitivity of 1°C means that a deviation of 10% on a parameter implies a deviation of 0.1°C on front temperature. We can see that the parameters of the third layer have less influence on the temperature than the first layer and the insulating layer. The sensitivities of this two first layer are inferior to 3°C. Maximum absolute values of sensitivity are those with respect to the front layer effusivity; it is negative at day and positive at night. x i 1

x i 2

x i 3

Fig. 6. Relative sensitivity of thermal properties of the three-layer wall

The optical parameters and the incident flux densities describe the radiative front excitation. In figure 6 we plotted the absolute values of the reduced sensitivity; it is actually negative for the albedo. The ratio between the sensitivities to albedo and to the VNIR flux varies from 5 to 12. This ratio should be about 3.6 when taking into account the actual value of the albedo. We observe a higher value because the sensitivity to VNIR variations as obtained through Eq. ( 9) is lower than the one obtained when considering a systematic VNIR relative variation. This deviation is due to the inertial effect of the wall: for the radiative flux, which is randomly scattered at high frequency about the reference values, the front temperature is less influenced than with a uniform deviation of the radiative flux which obviously has the same consequence that as an absorptivity variation.

Fig. 7. Relative sensitivity of convective (a) (at left) and radiative (b) (at right) parameters

We can see on figure 7.a the sensitivity of the front temperature to the convective coefficients and to the external and indoor air temperature. Overall, the front temperature has a low sensitivity to the indoor parameters. The external parameters and a systematic error on indoor temperature are most critical to the results. For instance, during the two days presented here, the sensitivity to the air temperature reaches about 300°C for the systematic variation analysis, respectively 100°C for the stochastic analysis. It means that an error of 0.5% on the air temperature (in Kelvin), i.e. an error of around 1.4°C, then implies an uncertainty of around 1.5°C, resp. 0.5°C, on the front surface temperature.

The front temperature is more sensitive to thermal and geometrical parameters of first two layers, to the outdoor convective parameters and to the IR radiation and emissivity. Therefore these parameters have to be measured with a high accuracy: the uncertainty has to be lower than about 5%, or 0.17% for air temperature (in Kelvin), to get an error lower than 0.5°C. Furthermore indoor parameters would be extremely difficult to evaluate from external surface temperature through data inversion.

Experimental

Description

The experimental measurement campaign BATIR took place in winter 2013 on the military airbase 701 at Salon-de-Provence, France. The thermal behaviour of two buildings was studied for a period of two weeks (see figure 7). Two meteorological stations were used: one station (Campbell CR23X) was placed in an open area next to the runway and the other one (Campbell CR3000) was placed between the two buildings. Meteorological data were collected in order to estimate outdoor air temperature, wind velocity, global solar irradiance and IR irradiance. Air temperature inside the rooms was measured by probes suspended close to back side of the facade of interest (EL-USB-1-PRO). The spectral reflectance of the participating surfaces was measured with a spectroreflectometer, in the [0.4-2.5 µm] domain (ASD) and in the [2.9-13.3 µm] domain (SOC400T). The thermography analysis was realised during two separate periods of 24 hours. Two 512x640 IR cameras were used: a band II camera (FLIR SC6000) and a band III camera (FLIR SC7000). The latter camera has been provided with a narrow filter. The half-height spectral bands are respectively [3-5 µm] and [8.4-8.7 µm]. The cameras were periodically aimed at an extended blackbody for drift correction and for establishing a relation between output signal and IR radiance in the respective bandwidths. The first 24h thermographic measurement was performed by aiming the south façade of building LC3 from the north façade of building LC4. For the second 24h thermographic measurement we performed the opposite. During the second thermographic measurement, five thermocouples and a fluxmeter were fixed on the north façade of building LC4 for independent values of the temperature and the heat flux at six particular locations of the wall. Thermography was also performed with a microbolometer camera onboard of a motor glider. The camera was flown at about 40 m.s

-1 and about 400 m altitude. It was aiming at nadir. During the first thermography period, two flights were performed short after sunrise, at midday, and during the second period, three flights were performed short after sunrise, at midday and shortly before sunset.

Experimental validation of SUSHI

The first experimental validation consisted in comparing our thermal model results with the thermocouple measurements on the north facade of the building LC4. The wall is made of three layers, consisting of a reinforced concrete layer, an isolating PS layer and a BA13 plaster board. The respective layer thicknesses and thermal properties, as obtained from literature are given in table 2. The outside convective coefficient was calculated from wind speed by the empirical formula, with wind speed measured at the building level. Wind speed data were available from two meteorological stations, one in an open area, and the other between the buildings. The first data series was available for whole month and the second one only for the last day of interest (started at 10:00). One day of data was thus missing for allowing the 48h analysis as described in § 2.2.1. Data from the first station were thus used for the missing day by applying a correction factor that was established by comparing the wind data recorded during the second, i.e. common day. This extrapolation is justified by the fact that the wind direction didn't change during this couple of days. The outside air temperature was given by the first station. The inside convective coefficient was arbitrarily fixed to 10 W.m -2 .K -1 , with air temperature measured at an approximate distance of 5 cm from the surface. The albedos were calculated from the measured spectral reflectance by weighting them with the solar spectral radiation provided by the atmosphere radiative transfer model MATISSE. The average emissivities were calculated from the IR measured reflectance by weighting them with the Planck's law corresponding to a surface at 5°C. The MATLAB version of the program was used here again, with boundary conditions fixed from 3D computations by SOLENE and SUSHI. Let us mention here that the net radiative infrared flux was obtained from the temperature 3D field estimated by the SOLENE thermal core; only the first group of iterations aimed at updating the convective heat flux term was thus performed; there was no updating of the radiative infrared flux term (see figure 9). Temperature measured with a given thermocouple is compared with the simulated temperature of the corresponding facet. These thermograms are shown in figure 9.

00:00 03:00 06:00 09:00 12:00 15:00 18:00 21:00 3 These images were recorded at 05:00 PM on January 29th. Two days were used to take into account the thermal dynamic of the scene. The scene was triangulated with 9975 facets. Areas with known thermal bridges (and visible on the IR images) were meshed more densely. For convenience, grounds, stairs, roof, windows and the door contained few facets (and the steps of the stair wasn't modelled). At the moment, half-transparency materials aren't processed by the programs; this is why temperatures and radiances of windows aren't realistic. 2D thermal bridges induced by floors were taken into account by the simulation but not the thermal bridges induced by indoor stairs (on the left of the façade) and the 3D areas induced by the floors near the right corner. The image of the relative difference shows that the radiances of 1D-wall parts are correctly estimated with a relative difference less than 2.5%. Heat losses induced by the thermal bridges at floor junctions are relatively well retrieved by the simulation despite the size of the facets used here. Differences observed on the ground and the roof are explained the significant roughness of these surfaces (they were considered flat in the 3D model). The difference reaches 9% at the soffit (the underside of the roof): the radiance of this part is higher than revealed by the simulation because the internal temperature can be higher than the measured temperature and the material is unknown. From this part some warm air could flow from the attic against the façade, which explains the difference on the top of the 3 rd floor.

Conclusion

In this work a numerical tool was developed for the simulation of the temperature field and the infrared rendering in an urban scene. The radiation energy budget at each facet of the modelled environment takes into account the multiples reflections of the incident solar flux (by SOLENE) and the IR radiative coupling between facets (by the new code SUSHI). The 1D response factors are calculated interactively by the quadrupole method whereas the 2D response factors are extracted form a database that was built from finite element calculation results. Computed results show that two days are necessary for catching the thermal evolution of a typical wall at 0.1°C precision though our periodic approach (the error is about 0.2°C for soils). The sensitivity study indicates the significance of thermal properties of the first two layers of a wall, and the influence of measurement errors about external air temperature, emissivity and infrared incident flux. Validations of this tool were performed by thermocouples and infrared thermography, with data acquired during the experimental campaign BATIR. Final results show little difference between measured and simulated radiance for the wall and the considered thermal bridges. Higher differences are observed for the roof (because the geometry of tiles) and the ground (because of its roughness) and at some locations for which boundary conditions were not well known (soffit and crawl space level). Further studies will be performed for analysing the airborne data.

Fig. 1 .

 1 Fig. 1. Flow chart of programs used in the full simulation

aT

  is air temperature, net vis R is the net radiative flux in the visible domain and net IR R is the net radiative flux in the infrared domain; it depends on the atmospheric IR radiation, on the temperature of the surface of interest and on the temperature of the surfaces exchanging radiation with it. This dependence on temperature is nonlinear which requires solving the problem through iterations. The time dependence of   t h

  t T sa corresponds to a modified sol-air temperature which is defined by:

  ', B', C' and D' are the reduced quadrupole coefficients (i.e. divided by the exponential of positive argument), fict ex h the fictive convective external coefficient, g in h the radiato-convective internal coefficient, abs VNIR  and his difference to the fictive convective coefficient, internal air temperature, i  the square root of time diffusion of the layer i,  the pulsation and j the imaginary variable. The temperature is obtained after an inverse Fourier transform of a next iteration. Periodically the radiative term k net IR  ~ is also updated as described in figure 2.



  related to all external facets (facets on S1 and S2 for the L-shape and facets on S1 for the T-shape) are then stored for future computations of the external temperature according to:

Fig. 8 .

 8 Fig. 8. Analysed scene during BATIR campaign (at left) and airborne IR image (at right)

Fig. 9 .

 9 Fig. 9. Comparison between temperature simulated by SUSHI and measured by thermocouple 4.3. Experimental validation of the full chain of simulation

Fig. 10 .

 10 Fig. 10. Comparison between radiance simulated by SOLENE/SUSHI/MOHICANS and measured by thermographyThe façade of the building called "LC4" is used for the validation of the full chain of programs. The figure10presents three images of this façade: (a) is the infrared image measured by the camera SC7000 (in the[8.4-8.7 µm] spectral band), (b) is the simulated image obtained at the outlet of the chain, (c) is the relative difference between the real and the simulated image and (d) is the vertical profile of radiance along the right side of the façade.These images were recorded at 05:00 PM on January 29th. Two days were used to take into account the thermal dynamic of the scene. The scene was triangulated with 9975 facets. Areas with known thermal bridges (and visible on the IR images) were meshed more densely. For convenience, grounds, stairs, roof, windows and the door contained few facets (and the steps of the stair wasn't modelled). At the moment, half-transparency materials aren't processed by the programs; this is why temperatures and radiances of windows aren't realistic. 2D thermal bridges induced by floors were taken into account by the simulation but not the thermal bridges induced by indoor stairs (on the left of the façade) and the 3D areas induced by the floors near the right corner. The image of the relative difference shows that the radiances of 1D-wall parts are correctly estimated with a relative difference less than 2.5%. Heat losses induced by the thermal bridges at floor junctions are relatively well retrieved by the simulation despite the size of the facets used here. Differences observed on the ground and the roof are explained the significant roughness of these surfaces (they were considered flat in the 3D model). The difference reaches 9% at the soffit (the underside of the roof): the radiance of this part is higher than revealed by the simulation because the internal temperature can be higher than the measured temperature and the material is unknown. From this part some warm air could flow from the attic against the façade, which explains the difference on the top of the 3 rd floor.

  Relative difference (%) (c)Radiance at sensor level (W.m -2 .sr -1 ) (b)Radiance at sensor level (W.m -2 .sr -1 ) (a)Radiance at sensor level (W.m -2 .sr -1 ) (d)

Table 2 .

 2 Description of wall

		Thickness and thermal properties
					)
	Layer	e (m)	k (W.m -1 .K -1 )	ρ.Cp (J.m -3 -1 ) .K	-1 .K b (J.s -1/2 .m -2	ξ (s 1/2 )
	Reinforced concrete	0.17	1.8	2.112e6	1950 184.1
	PSE	0.027	0.035 31.26e3	33.1	25.5
	Plaster (BA13)	0.013	1.0	1.44e6	1200 15.6

Acknowledgements

The authors thank PACA Region for its financial support, the laboratory CERMA for the provision of SOLENE, the laboratory IUSTI for the provision of the IR camera band II and the ONERA staff who participated to the measurement campaign BATIR.