J-J Roux

ARE 3D HEAT TRANSFER FORMULATIONS WITH SHORT TIME STEP AND SUN PATCH EVOLUTION NECESSARY FOR BUILDING SIMULATION?

A numerical model is developed to accurately simulate the transient thermal behaviour of rooms with sun-facing windows, with the use of a refined spatial and temporal discretization. For each node, the energy balance equations are developed based on a consideration of radiation, convection, air enthalpy and three-dimensional heat conduction. As buildings are exposed to rapid climatic variations (particularly incident solar radiation), we have added the different environmental conditions at short time-steps. The simulation considers the projection of solar radiation through a window onto interior walls, referred to as a sun patch. Therefore conduction transfer is treated in three dimensions. The indoor air temperature, the temperature of the cells in the walls and the surface temperatures are calculated at each time step using a variable-step Ordinary Differential Equation (ODE) solver. Results from this model are compared to well-known simulation tools using one-dimensional heat conduction without a sun patch.

INTRODUCTION

In France, the residential sector represents 42 % of energy consumption [START_REF] Ademe | Stratégie utilisation rationnelle de l'énergie[END_REF]. The French building code is guiding the construction sector to embrace low energy building practices. This regulation imposes limitations on maximal primary energy consumption, the use of the energy systems, and upper limits on air infiltration. In order to comply with these standards, low-energy buildings may incorporate high levels of insulation, heat recovery systems and passive solar building design techniques.

As low energy buildings are generally well insulated, overheating of rooms can be a serious issue. Therefore, it is important to assess the performance of current transient thermal models when adapted to low energy buildings, defined as those with particularly thick and insulated walls. Contemporary thermal models, like TrnSys (TRNSYS Mathematical Reference, 1996), Dymola (Dymola Version 2013) and Codyba (www.jnlog.com/codyba1.htm) are typically used to predict performance over the course of a year. As a result, certain features of heat transfer in building envelopes are simplified or neglected, for instance the location of a sun patch, or the treatment of heat conduction through a wall in one dimension only. Another drawback of TnrSys when performing short time step simulations of thick and highly insulated walls is its reliance on the Conduction Transfer Function method. Such problems have previously been reported by other authors (Flory-Celini, 2008;Delcroix et al., 2011;Savoyat et al., 2011). As a thermal model of a building envelope should take into account rapid environmental variations, we have decided to develop a new transient thermal model that provides a faithful description heat transfer through walls in 3D over a short timescale. This numerical three dimensional transient thermal model is introduced in the following section of this paper. In order to evaluate the importance of the novel features and the critical parameters of the model (control volume size, sampling rate, and sun patch simulation), the model was systematically compared to contemporary models under a range of simulation configurations.

MODEL DEVELOPMENT

The building envelope is modeled following a three dimensional approach. The model represents a single room with a large window. The thermal behavior to be included in the model is summarized as follows:

-Solar radiation exchange (short-wave and long wave) with the sun patch and the radiation distribution -Three-dimensional heat conduction -Enthalpy of a single air node -Convection between wall surfaces, the external and internal environment.

The walls of the model are simulated as a finite element mesh. Optimization of the meshing through the wall and its surfaces is described in a following section. Note that the model considers a single air node.

For each control volume, its temperature is given by the partial differential heat conduction equation in three dimensions, coupled to the sum of convective and radiative exchange :

(1) where

(2) and is the convective heat flow between the surface and its environment and and are the long wave and short wave radiation. Thermal conductivities in the three directions are designated = = [W/mK] and the volumetric heat capacity is denoted by [J/m 3 K]. The energy balance equation for the air in the room with temperature is:

(3)

Here is the exterior dry bulb temperature (°C), is the air flow (kg/s), the heat capacity of the air (J/kgK), the air density (kg/m 3 ), the volume of the room, the interior surface temperature and the convective transfer coefficient.

is the number of radiation balance equations corresponding to the number of surface mesh elements and the number of zones (here =1).

Short wave radiation

The total short wave radiation absorbed by the surfaces of a control volume is given by the vector the vector:

(4) where is the surfaces matrix and is the matrix of absorptivities of the control volumes of the internal wall for the short wave radiation. is the vector of radiation received by the mesh elements, obtained by resolving:

} (5)
Where is the reflectivity matrix. The view factors of the matrix are calculated following the Nusselt analog (unit hemisphere).

is the vector composed by primary radiation received by the control volumes. It results from the horizontal beam radiation and the diffuse radiation received by the mesh elements: = + , if the element i is in the sun patch, and

= if the element i is not included in the sun patch and are the direct and diffuse transmission coefficients of the glass and depend on the incidence angle of the beam, whereas Rb and Rd provide direct and diffuse radiation components on a titled surface with slope (Esveev andKudish, 2008)( Liu et al., 1969). The sun patch position has been calculated by a geometrical test: the boundary of the window is projected on an orthogonal plane to the beam. The control volumes of the walls are projected on the same plane, and thus those projected cells belonging to the projection of the window are identified. The short wave radiation received by the control volumes of external walls are determined depending on the position of the walls relative to the sun according to:

(7)
With Rr , and alb the ground albedo, which is the outside reflectivity of the ground.

Long wave radiation

Long wave radiation of the control volumes or the cells in the room of temperatures are given by: (

which is linearized to:

(9)

where is the unit matrix, the view factor matrix and emissivity matrix of the cells in the wall. The radiation coefficients were fixed to 5.8 W/m²K. Radiative heat transfer between external walls and the sky with effective temperature , and with the ground with temperature is given by:

(10) where are the radiative coefficients with the ground and the sky. These were also fixed to 5.8 W/m²K.

Convective heat transfer

Convective heat flows between the environment and the external surface of the walls are given by: (11) where W/m²K is the external convective coefficient. The convective heat flux with the internal environment is:

(12)
where the internal convective coefficient was assumed to be W/m²K

General resolution of the model

The input data are the building characteristics and the meteorological measurements. The numerical mesh of the room is generated with HEAT3 in cartesian coordinates [START_REF] Blomberg | HEAT CONDUCTION IN TWO AND THREE DIMENSIONS Computer Modelling of Building Physics Applications[END_REF]. Differential equations 1 and 2 are solved with a function of Matlab (ode23t). This function solves moderately stiff ordinary differential equations using the trapezoidal rule or Runge Kutta. This method is interesting because it solves the equations with an adapted and variable time-step. For fluctuant input data this method would automatically shorten the time step; similarly it would expand the time for periods with little fluctuation.

SIMULATIONS AND DISCUSSIONS

For all the simulations the same room dimensions were used: 3.01 m in width, 3.16 m in depth and 3.36 m in height. The walls, the ceiling and the roof are made of 10 cm thick concrete ( =2.7 W/mK) and the room has a 1 cm thick window located on the south-facing wall ( =1 W/mK)). Air flow in the room was neglected for all the simulations. Convective and radiative transfer coefficients are kept constant throughout the simulations.

The absorptivity and emissivity coefficients of the window and walls are listed in the table below.

Table 1: absorptivity and reflexivity of the wall and window

of the wall 1 of the window 0.05 of the wall 0.6 of the window 0.05 0.9 Ε 0.9

Minute-wise weather data were used from the meteorological station of Vaulx-en-Velin, France (latitude: 45°46'43''N and longitude : 4°55'21''E, height 170m) (http://idmp.entpe.fr/vaulx/mesfr.htm).

As a proof of principle test, the first model comprised the simple case of one control volume per internal surface and all radiation projected on the ground.

The model was subsequently transformed to more detailed configurations. As a first step, the importance of the weather data sampling was assessed. Then, the influence of cell size was analyzed and optimized for this room. Finally, simulations were performed for a time-varying sun patch.

First model: simplified hypotheses

For the first simulation the simplest case of one cell per surface was considered, neglecting reflections inside the room.

A visualization of this simple model, seen from the outside, is shown below generated in HEAT3 (Figure 1).

Figure 1: Shape and mesh of the room in HEAT3, the window appears in white

Each wall is described by only two cells. However, their connection to adjacent cells ensures that the heat transfer by conduction remains three dimensional.

For the same set of minute-wise weather data, the air temperature predicted by the model was compared to that obtained to Dymola (Dymola Version 2013) and TRNSYS (TRNSYS Mathematical Reference, 1996) in order to verify consistency for this simplified arrangement (Figure 2). The small discrepancies observed may arise from slightly different calculations of short and long wave, the one dimensional approach employed by TRNSYS and Dymola, and the solver algorithms. Differences between 0 and 1°C degrees are observed between our results and Dymola and differences between 0 and 2°C are observed between TrnSys and our model.

First model: sampling

Once we checked the general coherence of the model, we were interested in the sampling of weather data and its impact on the building simulation.

Weather data are often taken at the time step of one hour. This time step may be sufficiently accurate if the aim is to study the overall behavior of a building during a year. It is poorly adapted however, to transient simulations designed to study the detailed thermal behavior of the envelope. To overcome this, [START_REF] Escudero | Etude du comportement thermique des bâtiments. Proposition d'un modèle a deux dynamiques[END_REF] modeled the perturbations as mathematical functions and applied Shannon's sampling theorem. He concluded that to represent outside bulb temperature variation in time a 40 min time step is enough, and a 5 and 2 min time step are enough for the diffuse and direct radiation.

For the first more detailed simulation, we chose to keep the radiation data at the time step of one minute, in order to keep the short-time solar variations. So, the model needs to solve the differential equations ( 1) and ( 2) for each sample of weather data. For the second simulation, a set of hourly data was generated by integrating minute-wise data in the previous hour.

Simulation results using these aggregated data are shown superposed over those obtained with minutewise measured weather data in Figure 3. The same number of cells was used for the two simulations. In red, we can see simulations using one minute time step measured weather data as an input of the model and in green, results with hourly input data. Minute-wise data were retained for the following tests and simulations, in order to reveal any effects of these fluctuations on the numerical model. Note that this time step was not optimized. Indeed, a two or even three minute time step could be adequate.

First model: study of the spatial discretization

The importance of the cell size was first studied through the wall thickness. For a model, the greater the nodal discretization of the system, the finer the the calculated dynamics. Note that the time step of the thermal model and the narrowness of the spatial discretization of the described thermal system are interrelated.

When significant variations occur, for example fluctuations in solar radiation, the parts of the walls closest to the surfaces are the most affected.

Considering the perturbations studied, the mesh resolution was enhanced in the edge regions in order to provide more detailed information near the surface next to these gains. So, it seems obvious to refine the cells next to the surfaces in contact with the room's volume. Figure 4 shows the interior air temperature of the room for 2,4,8 and 13 cells and with a fine meshing near the surface.

Figure 4: Air temperature for different spatial discretization

With few control volumes, fluctuations in air temperature are smooth. As the control volumes are thicker, the first volumes near the room absorb more of the solar variations by conduction. The temperature differences, in green, vary from 0 to 2.3°C between a model with two and eight cells (Figure 5). These differences are largest when total incident solar radiation fluctuations are greatest. Indeed, we observe fewer differences for the third day, a day with less direct incident solar radiation.

For this wall thickness and with 13 cells, we observe very little difference with the model with eight cells.

In conclusion, 8 cells seem adequate and will be kept in the wall for the following tests as it seems enough to absorb the perturbations. More cells would cause longer calculation times. For the same hypotheses, we compared the previous result (with a finer mesh through the thickness of the wall) and a simulation with a finer discretization of the wall surfaces (Figure 6). The red scattered line is the air temperature variation in time for one cell per surface but eight cells perpendicular to the wall.

Figure 6: Influence of the spatial discretization of the walls' surfaces

In the figure, the blue line shows the results for a discretization of 144 cells per wall side and eight cells perpendicular to the wall. For this case, the three dimensional approach of the model has a greater influence on heat transfer. Thermal capacitance is more important as well as we observe that the response of air temperature to fluctuations in solar incident solar radiation. The air temperature differences between the two simulations are between 0 and 1.69 °C for the 4 th April.

Second model: sun patch integration Some of the gains which until now have been neglected are becoming more important in transient thermal modelling in low energy buildings [START_REF] Duforestel | Les outils de modélisation énergétique des bâtiments très basse consommation[END_REF]. It is the case of fast variation of incident solar radiation and their distribution or what we call sun patch.

The incoming incident solar radiation in a room is projected as a sun patch. In most of the modeling approaches, the incoming incident solar radiation is constantly projected on the floor. Other hypotheses are to project only 60% of the incident solar radiation on the floor.

However Wall (1997) shows that it is important to consider the sun patch distribution, especially for heating requirements of glazed spaces. Tittelein ( 2008) determines the sun position and projects the sun patch on the walls. He assigns to each wall a percentage of short wave radiation that varies in time. He shows that there is a difference of 8.5% for heating requirements between a model considering the sun patch and one without it. The next step now is to know whether assigning to each wall the same percentage of incident solar radiation is accurate enough. Indeed, a uniform distribution of solar heat flux on the walls seems adapted to one dimensional heat conduction; if we want to calculate accurately the impact of the perturbations it is important to locate the sun patch and therefore meshing the surfaces of the walls is necessary.

In order to locate the sun patch, the walls have been meshed in all the directions. In each direction we have 28 cells, 8 in each the wall's thickness and 12 on the surfaces. Internal reflections are integrated thanks to the view factors and absorptivity of the walls was fixed to 0.6 and reflexivity to 0.4. Finally, a three dimensional approach would be the best to do if we want to couple the sun patch movement and a transient thermal model. A better localization of the sun patch enables to locate overheating and the creation of plumes. These temperature differences will modify mean radiant temperature which is a comfort criterion.

The position of the sun patch depends on the position of the sun and the orientation of the building. On the figure below we show three moments in the day (morning, midday and afternoon) with the sun patch location (Figure 7 and8). We have represented the cells located in the inside of the room. In black, no sun patch touches the cells of the walls, in gray it is the window and in white we show the cells touched by the sun patch. Note that the gray cells represent the window and it is oriented south. We can see that for a simulation in April, the sun patch touches partially the floor, the wall oriented east and the one oriented west.

At 14h48 UTC on the 4 th April, we have the inside temperature distribution shown on Figure 9.At this moment of the day, the sun patch is mainly on the floor and on the wall oriented east.The sun patch is clearly located, because the cells receiving important short wave radiation reach 22°C. Reflections are observed on the adjacent cells. Temperature differences between one cell and another cell can reach 8.8°C. Besides, the impact of the solar sun patch inside the room has an important effect on heat conduction transfer inside the walls.For the same simulated day and time we looked at the crosssection of two walls by fixing a height (Figure 10). At this moment the air of the cell is at 18.54°C and the outside temperature is of 13.4°C. We can observe that the cells of the first figure (Figure 10), which are never in the sun patch, are quite homogenous in temperature. On the second figure we can see that at the left hand side we have high temperatures, which are located in the sun patch. There is clearly a temperature gradient, and the magnitude of the gradient rises towards the sun patch. Cell TB is never receiving direct solar radiation during this day but in the morning it will be hotter than the cell TA as it is located nearer to the sun patch at this moment. In the afternoon, when the sun patch is located on the opposite side, cell TA is warmer. Temperature differences for these two cells during one day are between 0 and 4.6°C. The room's air temperature (black scattered line) depends on all the surrounding cells' temperatures and therefore is located in between of the other curves. The heating of some areas of the wall due to the sun patch changes the transferred conduction fluxes to the outside. This is observed on the external surfaces where higher temperatures are observed where the nearby meshed volumes were touched by the sun patch (Figure 12) located on the inside of the eastern wall (Figure 9). We have compared the results from the first study done with a constant projection of the short wave radiation on the floor to this study with the sun patch model. In blue, we can see the room's air temperature when we project all the incident solar radiation on the floor with the same absorptivity and reflectivity than for the sun patch model. The black dotted line is the air temperature variation in time after integrating the sun patch. When integrating the sun patch we have the same incoming radiation through the window than with no sun patch, but the spatial distribution of the radiation touching the walls and control volumes are different. With the sun patch we have fewer cells touched by the direct radiation than when projecting the radiation on the floor. Therefore with the sun patch model, we tend to have higher temperatures for the cells touched by the sun patch but we have a smaller area than when supposing that the radiation touches the whole floor's surface. In average, with the sun patch model we have temperature differences between the surfaces and the air which are smaller than with the other model. That is the reason why, the air temperature has a smaller gradient when we add the sun patch (Figure 13). Over all, we obtain a less fluctuant temperature variation with the sun patch integration, than with the simple case. Differences between 0 -1.09°C are observed. We have started with a simplified model and with current taken hypotheses: few cells, all the short wave radiation on the floor and hourly weather data.

We have gradually transformed this model in a more complex model, by working on the sampling rate, the spatial discretization and the sun patch simulation.

In conclusion, taking hourly weather data can lead to smoothing too much the solar radiation fluctuations and it has an impact on the simulated results.

A minimum number of cells were necessary to model in three dimensions the heat transfers if we want to use short time steps. By adding the sun patch position moving in time, we have observed differences concerning the surfaces temperature distribution, the incident short wave length radiation calculation and therefore differences of the air of the room have been observed. In conclusion, if a detailed envelope model is required in the building area, using a 3D model and considering the rapid weather variations in time with the sun patch leads to important improvements, which have been quantified.

FUTURE WORK

These results have to be validated with experimental data, and therefore we are working on an in situ experimental set up.

Figure 2 :

 2 Figure 2: Comparison between TRNSYS -Dymolaour model

Figure 3 :

 3 Figure 3: Impact of hourly and minute time-step weather data on the air temperature simulation Short term variations in total radiation are clearly visible on the blue curve also shown in the figure.The weather fluctuations have an impact on the results simulated with the model for weather data with a time-step of one minute. Hourly data tend to smooth the fluctuations of the solar radiation and this has an impact on the building simulation results. Differences up to 1°C were observed on the air temperature calculation of the cell for the data sample.

Figure 5 :

 5 Figure 5: In black global incident solar radiation and in green the air temperature differences

Figure 8 :

 8 Figure 8: Sun patch location in the afternoon

Figure 7 :

 7 Figure 7: Sun patch location for the 4 th April: morning and midday

Figure 9 :

 9 Figure 9: Inside temperature distribution of the floor, window and the wall oriented east at 14h48 UTC

Figure 11 :

 11 Figure 11: Temperature of cell TA and TB and the temperature of the room's air

Figure 12 :

 12 Figure 12: Outside temperature distribution of the floor, window and the wall oriented east at 14h48 UTC on the 4 th April

Figure 13 :

 13 Figure 13: Air temperature of the cell for two simulations: with and without sun patch

CONCLUSION

ACKNOLEDGEMENTS

This work has been supported by the French National Agency for Research (ANR), in the framework of the project "SUPERBAT", ANR-10-HABISOL-004-06. I would also like to thank Leon Gaillard and Pierrick Haurant for their contributions.


Manuel Sánchez

Montero Responsable

Extremadura, Ceuta La O N G A C C E M En Andalucía

Melilla

Mercedes Guzmán Bejarano

A C C E M Sevilla

pre

El Observatorio Permanente de la Inmigración en Sevilla nació en Febrero de 2006 con el objetivo de crear una estructura de partenariado entre actores relacionados con el colectivo inmigrante tanto a nivel de intervención social como de planificación de políticas sociales.

Asimismo, a través de este dispositivo de observación cooperativa y de evaluación participativa se pretende dotar a los diferentes actores implicados de herramientas que permitan acercarse a la realidad de manera científica, analizarla y proponer estrategias de intervención en consonancia con las problemáticas y necesidades sociales en las que se ve inmerso el colectivo inmigrante.

La necesidad de articular este observatorio surge a propuesta de A.C.C.E.M., ante la conveniencia de que la toma de decisiones (a nivel de políticas sociales y económicas) en cualquier ámbito debe estar basada en un conocimiento profundo y dinámico de esta realidad de la población inmigrante. Este conocimiento permitirá confirmar y verificar la viabilidad de los planteamientos e iniciativas propuestas (planes, programas y proyectos de intervención social), así como su adecuación al público al que van dirigidos.

En sus inicios, el observatorio estaba compuesto por cuarenta y dos socios, y en la actualidad el número se ha duplicado llegando a ochenta y ocho entidades, en las cuales están representados los tres sectores (público, privado, y tercer sector compuesto por ONGs y sindicatos). Entre ellas se encuentran administraciones públicas a nivel regional y municipal, servicios públicos sociales, sanitarios y educativos, agrupaciones de empresarios, asociaciones de inmigrantes y pro-inmigrantes, sindicatos y federaciones de vecinos.

El territorio que abarca el dispositivo también se ha visto ampliado, ya que en un primer momento se centró en dos barrios de Sevilla, La Macarena y Macarena Norte, donde el porcentaje de población inmigrante es sensiblemente superior al del resto de la ciudad. Fue en 2007 cuando se dispuso que el área de observación se extendiera a toda la ciudad de Sevilla con vistas a tener una visión más global, pues en ningún caso los resultados de un estudio sobre esas dos zonas, se podrían extrapolar a la ciudad completa. A continuación se ofrecerán una serie de datos que permitan situar el contexto territorial al que hacemos referencia. Sevilla está situada al suroeste de la península Ibérica. Dentro del mapa político de España, pertenece a la comunidad autónoma de Andalucía de la que es capital. La ciudad de Sevilla cuenta con una población de 706.289 habitantes, de los cuales un 8,6 % son extranjeros. A continuación se muestran en una tabla la distribución por continentes de la población extranjera en Sevilla. Los principios sobre los cuales se asienta este dispositivo permanente de observación son: participación ciudadana, aproximación global al territorio, partenariado de los actores locales, acceso y uso de las nuevas tecnologías, accesibilidad a la información y calidad de la información.

Unión

Una vez finalizada la presentación de los objetivos, principios y componentes del partenariado así como la contextualización del territorio; pasaremos a tratar las acciones del Observatorio.

Teniendo en cuenta que con anterioridad a la creación del observatorio, no existía ninguna red que aglutinara a entidades de diversos sectores en torno a la realidad de la inmigración en el municipio de Sevilla, es posible comprender que en un primer momento el interés de los socios se centrara en establecer una estructura estable de partenariado (Para ser más exactos, hay que decir que antes del Observatorio habían surgido algunas iniciativas similares, pero ninguna de ellas tuvo éxito). A raíz de ese deseo de aprovechar al máximo las oportunidades que brinda el trabajo en red, dentro del Observatorio se conformaron a petición de los socios varios talleres (empleo, educación, salud y sensibilización) que continúan estando activos en la actualidad, y a los que se ha unido el taller de vivienda de reciente creación.

A continuación explicaremos los objetivos de cada uno de ellos y las acciones realizadas.

-El taller de empleo se formó con 35 socios. Pretende coordinar el trabajo entre las entidades para dar respuesta a las diferentes necesidades relacionadas con el empleo y la inmigración. Podemos destacar las siguientes acciones que se han llevado a cabo:

• Establecimiento de una tabla salarial para el sector del Servicio Doméstico que se utiliza en las entidades con bolsas de empleo. Este acuerdo ha tenido repercusión en otras provincias andaluzas, como es el caso de Córdoba, donde se ha firmado un convenio entre todas las entidades con bolsas de trabajo para personas inmigrantes en el que se comprometen a aplicar esas tablas salariales.

• Acceso a los cursos de formación del Servicio público de empleo por parte de los solicitantes de asilo en espera de resolución. Esto fue posible gracias a una instancia firmada por el taller de empleo y presentada en la Conserjería de Empleo.

• Después de obtener el listado oficial de profesiones "de difícil cobertura", se ha realizado un estudio de los diferentes convenios de trabajo con vistas a la elaboración de una tabla con las condiciones mínimas de contratación.

• Creación de un sistema de traspaso de ofertas de trabajo entre las entidades con bolsa de empleo, con el objetivo de que éstas no quedaran desiertas.

-El taller de educación se compone de 26 socios. La finalidad de este grupo de trabajo es reflexionar sobre las necesidades que se detectan en el ámbito de la educación y la inmigración. Los socios de esta mesa destacan el análisis y la puesta en común de los planes -de acogida y convivencia-de los centros escolares situados en zonas de Sevilla con mayor porcentaje de inmigrantes. También se han organizado Jornadas de puertas abiertas en los centros educativos con el objetivo de lograr la implicación y participación de las familias en el proceso educativo de sus hijos y más concretamente para que conozcan a los profesionales y los recursos de cada centro.

-El taller de salud está compuesto por 10 miembros. El objetivo principal es la coordinación de los diferentes organismos sanitarios (tanto del ámbito público como del tercer sector) que atienden a la población inmigrante para desarrollar actividades referentes a este ámbito. Se ha elaborado una guía de recursos de los diferentes servicios que las entidades prestan, para atender a la población inmigrante de una manera más ágil y completa.

-El taller de sensibilización está formado por 26 entidades. El objetivo principal de este grupo de trabajo es desarrollar actividades que propicien una adecuada relación entre las personas inmigrantes y la sociedad de acogida, dando a conocer la riqueza de sus valores, culturas, tradiciones, etc. Conforme a esto, se han organizado eventos que tenían como finalidad la celebración del día de la inmigración en los dos últimos años, 2006 y 2007. Conviene subrayar que, gracias al trabajo en red realizado entre los socios del Observatorio, se ha pasado de organizar cada entidad su propia celebración para este día, a llevar a cabo una Fiesta organizada conjuntamente entre todos los socios miembros del Observatorio.

-El taller de vivienda es el de más reciente creación. Los objetivos que pretende desarrollar son: crear una guía de recursos que permita facilitar información a ese respecto, y sensibilizar a los propietarios de viviendas que tengan intención de alquilarlas. Con vistas a conocer previamente las dificultades que encuentran cada una de las dos partes (la que desea encontrar un piso en alquiler y la que tiene intención de alquilar su vivienda), se ha acordado realizar un estudio en una zona de la ciudad donde se está empezando a establecer población inmigrante.

Como ya hemos dicho anteriormente, estos cinco talleres han nacido como fruto de la estructura de partenariado. Igualmente cabe destacar que el trabajo en red ha propiciado un mejor conocimiento y una mayor confianza entre las entidades que conforman el Observatorio. De ahí ha surgido la oportunidad de presentar proyectos en red entre varias entidades para solicitar subvenciones a las diferentes administraciones públicas.

Este año, en cambio, se han centrado los esfuerzos en la recogida y procesamiento de datos con vistas a realizar un diagnóstico de las necesidades de los inmigrantes atendidos en 19 entidades que forman parte del Observatorio. Ante todo, hay que subrayar los cursos de formación que han recibido los técnicos a cargo del Observatorio y que han sido cruciales para su capacitación. Estos cursos han sido impartidos por el equipo técnico de la Universidad de Besançon y han versado sucesivamente sobre Análisis Cuantitativo y Cualitativo, y Diseño Cartográfico utilizando diversas herramientas informáticas.

En los primeros meses del presente año, no estaba todavía establecida en la sede de A.C.C.E.M. Sevilla una dinámica de trabajo que permitiera una recogida masiva de datos. Pero, afortunadamente, se ha producido un cambio, y ahora la recogida está aumentando de un modo progresivamente acelerado. Una meta todavía sin alcanzar es la de implicar más activamente al resto de los socios en la tarea de recogida de información. Quizás el énfasis del pasado año en el conocimiento mutuo y el estrechamiento de lazos entre las entidades pertenecientes al Observatorio, distrajo la atención de la tarea básica de recogida de datos. El mantenimiento de esta actitud es un obstáculo con el que actualmente nos estamos enfrentando, y es de prever un cambio favorable muy cercano. En las reuniones plenarias el personal técnico a cargo del Observatorio está haciendo hincapié en la necesidad de reunir la mayor cantidad posible de datos: "La recogida de información es un paso sin el cual es imposible planificar, definir y evaluar cualquier política de desarrollo territorial sostenible".

Al lado de esa insistencia ante los responsables de las distintas entidades, se han llevado a cabo visitas a las diferentes sedes con vistas a formar al personal técnico y administrativo en el manejo de los formularios. En ocasiones también se ha instalado la herramienta informática de recogida de información (Pragma) en los ordenadores de esas sedes. Asimismo, con los socios pertenecientes al grupo diagnóstico se realiza un seguimiento sobre la recogida de datos, a la vez que se plantean alternativas si resulta necesario. En el caso de que una entidad no cuente con recursos humanos suficientes o que perciba como una sobrecarga inasumible esta nueva tarea, el resto de las entidades del grupo diagnóstico ponen a su disposición algún técnico que acuda determinados días a la entidad que lo ha solicitado. A pesar de esto, a día de hoy, los datos recogidos por otras entidades son todavía escasos y en su inmensa mayoría en formato papel. Como se dijo arriba, cabe esperar que estas circunstancias mejoren. Cuando ello se produzca, seguramente se podrá dar el importante paso adelante que supondría establecer un calendario de entrega de guiones recogidos.

A fecha de Septiembre del 2008 se han superado los 300 guiones recogidos y se espera reunir 500 hacia finales de año. Se han llevado a cabo dos controles de calidad de los datos recogidos. El primer control de calidad ha llevado a las personas implicadas a comprender la importancia de una recogida extremadamente cuidadosa.

Además, se está diseñando el plano de la ciudad de Sevilla dividido en distritos y secciones en el que se volcarán los datos territoriales proporcionados por las instituciones estadísticas oficiales utilizando programas de cartografía automática (Philcarto y Phildigit). Todo este material se utilizará durante la presentación de los resultados del diagnóstico que tendrá lugar a comienzos del 2009. Mediante esta contextualización podremos comparar los datos resultantes del diagnóstico con los del territorio y en último término verificar, o refutar en su caso, las hipótesis surgidas a raíz del diagnóstico. En definitiva, nos introduciremos de este modo en un proceso de aprendizaje sobre el territorio en cuestión.

Resto de África 2.401 América Central 2.415 América del Norte 886 América del Sur 20

	Europea	20.959
	Europa no Comunitaria	3.441
	Magreb	7.518
		.632
	Asia	3.020
	Oceanía	48
	Apátridas	19
	Total	61.339


Bertrand Putigny
email: bertrand.putigny@inria.fr

Brice Goglin
email: brice.goglin@inria.fr

Denis Barthou
email: denis.barthou@inria.fr

A Benchmark-based Performance Model for Memory-bound HPC Applications

Keywords: memory model, timing prediction, microbenchmarks, caches, multicore

The increasing computation capability of servers comes with a dramatic increase of their complexity through many cores, multiple levels of caches and NUMA architectures. Exploiting the computing power is increasingly harder and programmers need ways to understand the performance behavior.

We present an innovative approach for predicting the performance of memory-bound multi-threaded applications. It relies on micro-benchmarks and a compositional model, combining measures of microbenchmarks in order to model larger codes. Our memory model takes into account cache sizes and cache coherence protocols, having a large impact on performance of multi-threaded codes. Applying this model to real world HPC kernels shows that it can predict their performance with good accuracy, helping taking optimization decisions to increase application's performance.

I. Introduction

High performance computing requires proper software tuning to better exploit the hardware abilities. The increasing complexity of the hardware leads to a growing need to understand and to model its behavior so as to optimize applications for performance. While the gap between memory and processor performance keeps growing, complex memory designs are introduced in order to hide the memory latency. Modern multi-core processors feature deep memory hierarchies with multiple levels of caches, with one or more levels shared between cores.

Performance models are essential because they provide feed-back to programmers and tools, and give a way to debug, understand and predict the behavior and performance of applications. Coherence protocols, prefetchers, replacement policies are key automatic hardware mechanisms that improve peak performance, but they also make parallel architecture modeling much harder. Besides the detailed algorithms used by these mechanisms, the latencies involved or the size of the buffers implied are not precisely quantified by hardware manufacturers in general. Finally, even if the hardware resource organization may be known at runtime thanks to tools such as hwloc [START_REF] Broquedis | hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF], the way they are actually used depends on the application, namely how it allocates, reuses and shares buffers between the computing tasks.

These difficulties hinder precise performance modeling, in particular for memory performance. Many previous works have abstracted away the cache behavior by modeling essentially compulsory, capacity and conflict misses (the 3C). The resulting analytical models, through the computation of stack distances, capture data locality and offer a guidance for compiler optimizations (see [START_REF] Vera | Efficient and Accurate Analytical Modeling of Whole-Program Data Cache Behavior[END_REF], [START_REF] Cascaval | Estimating cache misses and locality using stack distances[END_REF], [START_REF] Andrade | Accurate prediction of the behavior of multithreaded applications in shared caches[END_REF] for instance). Most of these previous works target the cache behavior of single-threaded codes. Some recent works focus on how to model cache misses [START_REF] Andrade | Accurate prediction of the behavior of multithreaded applications in shared caches[END_REF] for multithreaded codes, or on how to model cache coherence traffic [START_REF] Lee | Thread Tailor: Dynamically Weaving Threads Together for Efficient, Adaptive Parallel Applications[END_REF]. However, they do not consider simultaneously the impact of cache misses, of coherence and contention coming from shared caches. A 5C model, accounting for compulsory, capacity, conflict misses, contention and coherence remains to be found in order to help study the impact of factors such as the data set size for each thread or the code scalability on multi-cores.

In this paper we present a novel performance model that allows detailed performance prediction for memorybound multi-threaded codes. This differs from the previous approaches by predicting the cumulated latency required to perform the data accesses of a multi-threaded code. The model resorts to micro-benchmarks in order to take into account the effects of the hierarchy of any number of caches, compulsory misses, capacity misses (to some extent), coherence traffic and contention. Micro-benchmarks offer the advantage of making the approach easy to calibrate on a new platform and able to take into account complex mechanisms. In addition these benchmarks can be used as hardware test-beds, e.g. to choose between several architectures which one best suits the needs. This paper is organized as follows: first, Section II presents the scope and an overview of our work. Our model is later detailed in Section III as a combination of hardware and software models. Then Section IV presents the benchmarks used to build the hardware models. Finally Section V details the use of our model to predict real world code performance.

II. Scope and Model's Overview

Our model takes as input the source code to analyze. The code can be structured with any number of loops and statements, and we assume the data structures accessed are only scalars and arrays. Parallelism is assumed to be expressed in OpenMP parallel loops. The iteration count of the loops have to be known at the time of the analysis, and the array indexes have to be affine functions of surrounding loop counters and of thread ids.

Predicting the time to access memory usually requires to build a full theoretical performance model of the memory hierarchy. This work is however difficult due to the complexity of all the hardware mechanisms involved in modern architectures. Instead we choose to build a memory model upon benchmarks that are used to capture the hardware behavior of common memory access patterns while hiding the impact and complexity of these mechanisms inside the benchmark output.

The main difficulty with this approach is to find a set of benchmarks that characterizes hardware precisely, and to keep this set as small as possible. Then we try to rebuild the application memory access pattern by combining the outputs of these basic benchmarks. We found that in cache coherent architectures, the state of the target cache lines has a large impact on performance. Concurrent accesses to shared data buffers lead to cache-line bounces between cores due to the need to maintain coherence between the existing data copies in their caches. Thus we build a set of benchmarks that gives us insight about the read and write latency to cache lines for every of the state of the MESI protocol [START_REF] Papamarcos | A low-overhead coherence solution for multiprocessors with private cache memories[END_REF]. Indeed most cache coherent processors use protocols that are based upon MESI.

In order to predict the time needed to access memory for a given application, we decompose it into a memory access pattern. This pattern tells us the amount of memory access and how it is accessed (i.e. reading or writing). If we suppose that we know the full state of memory (i.e. the cached addresses and their locations), we can i) reconstruct the state of memory after every access, and ii) read each access duration from the output of our benchmarks. By taking memory events one after another we can track the state of data in caches and construct a formula that will predict the time needed to access memory for the whole application.

III. Program and Memory Models

In Section III-A we present how programs are represented in order to be able to apply the prediction. Then we detail the view of memory used in our work in Section III-B. Finally Section III-C describes how the model is used for predicting execution times.

A. Program Model

Input codes are OpenMP parallel codes, using only parallel for loop constructs. Our model represents programs by only considering their memory access patterns since our work focuses on memory-bound codes. OpenMP codes are transformed into a simplified code, where statements are memory accesses, each statement accessing a memory region.

Read/write array accesses in OpenMP parallel for loops are transformed into read or write statements accessing all the elements read or written by a thread. The parallel for loop is then transformed into a parallel OpenMP region. For instance, the following DAXPY computation with n threads, assuming N is a multiple of n:

double X[N] ,Y[N ] ; #pragma omp p a r a l l e l f o r f o r ( i =0; i <N; i ++) Y[ i ] = a * X[ i ] + Y[ i ] ;
is modeled as the following code: In this representation, the sequence of reads and writes, where initially two reads alternates with one write, has been replaced by two streams of reads followed by a stream of writes. In order to represent the capacity of modern architectures to load and store multiple elements per cycle, we consider an execution model where read and write statements can be dispatched to parallel units. For instance, on an architecture able to perform one read and one write at a time, the two read accesses can only be performed in sequence, while the write access can be performed concurrently to the reads. This corresponds to the asymptotic optimal use of the hardware units, that can be reached through code optimization.

Similarly, element-wise read/write accesses in sequential loops are subsumed into array region accesses, using read and write statements. For instance, the code on the left hand side is modeled as the code on the right hand side.

f o r ( i =0; i <N; i ++) X[ i ] = k ; -> write ( [ 0 : N-1 : 8 ] ,X ) ; k = Y[ i ] ;
read ( [ 0 : N-1 : 8 ] ,Y ) ;

In the original code, the first statement depends on the second statement of the previous iteration. This loopcarried dependence is not represented in our model. In the execution model proposed above, both read and write array region statements will be executed simultaneously.

That represents the fact that on modern hardware, an optimized version of the original code is still able to execute one load (X[i+1]) and one store (Y[i]) per cycle.

In the previous code there is only one data chunk for X, while in the parallel DAXPY code, the same region corresponds to N/T data chunks. In our memory model (presented in the following section), we keep track of the state of each data chunk in the coherency protocol. To do so, the data chunks accessed by a program have to be considered as atomic pieces of data, in particular, two different chunks cannot share array elements, as for array X and Y in the previous examples. This can be obtained by decomposing read and write statements into statements accessing smaller data chunks, so that the chunks represent a partition of the array elements. For instance, in order to manipulate the same data chunks in the previous code and in the parallel DAXPY, the two read and write statements in the previous code are decomposed into N/T statements:

f o r ( t =0; t<T ; t++) write ( [ t * N/T : ( t +1) * N/T-1 : 8 ] ,X ) ; f o r ( t =0; t<T ; t++) read ( [ t * N/T : ( t +1) * N/T-1 : 8 ] ,Y ) ;
where T denotes the number of threads used by the program. This assumes that the number of threads is constant for the execution of the program. Splitting array regions can lead in the extreme case to one region for one array element. We assume this operation does not create a number of region dependent on the size of the array.

This representation only handles exactly array index regions that can be represented by triplets. For OpenMP, this is limited to parallel loops with static scheduling, MASTER sections (only thread 0) or SINGLE (can be modeled with thread 0). For other cases, or for triangular sequential loops, the chunks are then over-approximated.

B. Memory Model

The objective of the memory model is to predict the time required to perform the access statements, taking into account the impact of the coherency protocol. To this end, we associate a state in the coherency protocol to chunks and tabulate the time to perform an access statement to a chunk of any given size, with any number of threads and in any state. As the impact of the cache hierarchy is taken into account by a latency function, the memory hierarchy is entirely modeled as one level of

Cache (∞) Core #0 Cache (∞) Core #1 Cache (∞) Core #2 Cache (∞) Core #3 Cache (∞) Core #4 Cache (∞) Core #5
Figure 1: The memory hierarchy is modeled as one level coherent, private and infinite capacity caches. coherent, private and infinite capacity caches, with one cache per core, as depicted in Figure 1. We will assume in the following that threads are pinned to cores. A thread has therefore a corresponding cache. We define a latency function, giving the time to access a chunk of data as a function depending: on the size of the chunk and on its stride, on the state of this data in the caches, and on the number of threads that access data simultaneously.

The state of a chunk consists in one of the MESI states and on the list of threads that have the chunk in their own cache. Adapting this technique to variants of the MESI protocol, such as MESIF or MOESI, is straightforward. These chunk states can be: The transitions between these states correspond to the actions performed on the chunks. A chunk can be either read by a set of threads t ∈ T (action denoted L T for load) or written by one thread t (action denoted S t for store). Figure 2 defines the transition function δ for these states, according to the type of access and the ids of the threads accessing the chunk. For instance, consider a chunk in state Exclusive, in the cache of thread 0. This state is denoted (E, {0}). Writing this chunk with thread 1 changes its state to δ((E, {0}), S 1 ) = (M, {1}). Each chunk maintains its own state and each state keeps track of the threads ids that have a valid copy of the chunk in their cache. At the beginning of the program, the data chunks are in the state (I, ∅). Executing read and write statement in our program model changes chunk states and the latencies for these transitions are defined by micro-benchmarks.

M
M, {t} E, {t} S, T I, ∅ S t Su, u = t, t ← u L {t} L T , T = {t} S t Su , u = t, t ← u L T , T = { t } L {t} S t S t , t ∈ T Lω , ω ⊂ T Lω , ω ⊂ T , T ← ω ∪ T S t L {t} L T , |T | > 1

C. Time Prediction

The prediction of the memory access time of an OpenMP parallel code relies on two steps: The time measurement of all transitions in Automaton 2, for all data sizes and thread numbers, and the definition of a semantics connecting the read and write statements with the transitions within this automaton.

All transitions of the Automaton 2 are associated to timing functions obtained through micro-benchmarks. These transitions depend on the type of access (read/write), and of the number of threads performing this same transition. Indeed, some read transitions L T can be taken by threads in a subset T of threads. The time to execute a transition depends additionally of the initial state of the chunk, its size and stride, and of all transitions taken by other threads concurrently (sharing bandwidth to memory and last level caches). In order to limit the number of microbenchmarks, the timing function considered for a transition will only depend on the initial state of the chunk, its size and stride, the total number of threads, assuming they are all performing the same transition, on a different chunk with same characteristics. These benchmarks are thus taking correctly into account capacity misses and limited bandwidth when threads access simultaneously different chunks. However, the latency of a parallel access to the same data chunk is not considered and will be approximated by a parallel access to multiple chunks of the same size. Considering a data chunk c read by thread t in a parallel region with T threads, the time to perform this read will be denoted L(c.state, c.stride, c.size, T ). The description of how the benchmarks are built is detailed in the following section.

The timing functions provided by the benchmarks are then used to build a timing function for the program. The states are computed after each read and write statement. To aggregate these timing functions, we consider separately sequential phases in the program from parallel phases (within OpenMP parallel loops). Besides, two considerations are taken into account:

• At least one read and one write memory access can be executed in parallel on modern architectures. For each phase, the total duration for a sequence of memory access statements is computed as the maximum of the aggregated time to perform the reads, and of the aggregated time to perform the writes. • A sequence of read/write chunk statements can represent element-wise accesses that are interleaved. We consider therefore the size of the data accessed during a read/write chunk statement as the total size of all chunks accessed during this phase.

12 θ = max t=1..T (θ L (t), θ R (t))
of accesses, the total working set size for each thread is first computed (line 1). Read/write accesses are assumed to be on data of this size (parameter s(t) in L,S lines 7, 9), and the time is scaled down to the actual size of the chunk. States are then updated according to the type of access.

Read (and write) statements can be executed in parallel. As all chunks may not be in the same state, these memory accesses may not take the same time for all chunks. Zhang et al. [START_REF] Zhang | Does cache sharing on modern CMP matter to the performance of contemporary multithreaded programs?[END_REF] have shown that the sharing pattern of threads on multi-threaded application are often very regular: in this case all chunks accessed within a read or write in a parallel phase are in the same state. In order to ensure prediction even when different chunks are in different states, the latency for all chunk accesses is simply assumed to be the longest one among all accesses (max in line 12).

IV. Benchmarking Framework

This section presents the benchmarking framework we developed [START_REF] Putigny | Mbench: memory benchmarking framework for multicores[END_REF]. In the following, memory latencies are presented as bandwidths. These metrics are equivalent since the data size accessed is known but bandwidth makes performance comparison easier on graphical plots.

A. Benchmark Design

We build a set of elementary benchmarks to measure the performance of each memory level (i.e. cache levels and main memory). In order to achieve near peak memory performance, these codes are in assembly and perform read or write access to a buffer of a given size. These building blocks are used to time each of the 15 transitions of automaton 2. To do so, these benchmarks first set a buffer into the initial state of the automaton, then access the buffer (by a read or write) and record the bandwidth achieved. Moreover, another set of benchmarks evaluates the possible slowdown due to contention when the transitions are taken in parallel by several sets of threads. Our approach differs significantly from other tools such as the Stream benchmark [START_REF] Mccalpin | Memory Bandwidth and Machine Balance in Current High Performance Computers[END_REF] because it takes coherence and possible contention into account.

The performance of benchmarks involving either the state (S, T ) or a transition L T depends on the number of caches holding the same chunk of memory or the number of threads accessing this data, hence it depends on their location in the machine. For instance, to execute the transition from state S to state M with a store requires to use at least 2 threads on a 8-core socket, hence there are 28 possible configurations. The same benchmark with 3 threads in parallel theoretically requires 56 tests. In order to avoid the combinatorial explosion of running all these configurations for every benchmark, only one configuration is tested by benchmarks for a given number of threads. This leads to n different runs total instead of to 2 n-1 if we ran all combinations for all amounts of threads. Fortunately the performance of a benchmark with a given number of threads does not depend on their location within the socket since the inner-socket architecture of our target platform is flat.

B. Benchmark Outputs

All results in the remaining sections of the paper are obtained on a dual-socket 8-core 2 GHz Xeon E5-2650 host based on the Intel Sandy Bridge architecture.

Figure 3 presents the timing for the transitions corresponding to the same action, a load L {u} when u is not among the threads of the initial state (either M, E or S). This transition hence corresponds to a miss, with possible coherency messages. Note that in this case, counting the number of cache hits and misses, as commonly done [START_REF] Pesterev | Locating Cache Performance Bottlenecks Using Data Profiling[END_REF], is not sufficient to predict performance. Indeed, for these three configurations, the numbers of hits and misses are exactly the same, while the effective latencies are very different. Latency highly depends here on the initial coherency state of the data. There is almost a factor 2 difference between a load to modified data (6-8 GB/s for up to 64KB of data) w.r.t. a load to shared data (14-18 GB/s for the same range), when this data fits in L1. The reason is that making a read request for data in the state M in another L1 cache requires first to invalidate it in the L1 (and possibly L2). When data is large enough, it moves during the write from L1 to the L3. The read then shows L3 performance (around 18GB/s for sizes larger than 8MB).

•••• • • • • •••• • • • • •• •• • • • • • • •• • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • •• ••• • • • ••• •• • • • •• • • • • • • • • ••• 6
To assess the effects of contention occurring when the same caches are accessed concurrently simultaneously, we generate benchmarks where all threads run the same code, on different data. Threads do not share data, only bandwidth and the capacity of the last level cache. For a given transition, we analyze the performance degradation defined as the ratio between the bandwidth measured for one benchmark over the bandwidth measured when n similar benchmarks are running in parallel. A ratio of 1 means that there is no contention, each thread can use the same bandwidth as it would if it was running alone on the machine. A ratio greater than 1 is the factor dividing the available bandwidth per thread. This ratio does not necessarily represent contention within caches or on the memory bus. We actually observed no cache contention on the Intel Sandy Bridge microarchitecture used for our tests, while the AMD Bulldozer micro-architecture shows some. The performance loss can also result from the limited capacity of the shared L3, observed when threads saturate it. They cause some parallel capacity misses, which looks like cache contention on the benchmark outputs.

Figure 4 shows the bandwidth degradation ratio when a store hit occurs. As shown on Figure 4 there is almost no contention on private caches: Every independent cache can deliver the same bandwidth when it is accessed alone or when all private caches are accessed at the same time. This is particularly interesting because no contention appears even when coherence traffic is involved. However, as one would expect, contention seems to appear in shared resources such as L3 and memory. As explained above, the L3 cache contention is actually caused by parallel capacity misses caused by multiple threads sharing the overall L3 size.

••••• • • • ••••• • • • • •••• • • • ••••• • • • ••••• • • • ••••• • • • ••••• • • • ••••• • • • • • • • • • • • • ••• • • • • •• • • • • • • • • • • • • • • •• ••• 0 2 4
More interestingly the ratio depends on the state of cache lines accessed: accessing modified cache lines (in L3) exhibits less performance than accessing clean lines. Also, we can see on Figure 4 that performance is more degraded by writing to exclusive cache lines than by writing to shared ones. We could not explain this behavior and it justifies further the idea to hide hardware complexity by using benchmarks: they capture such puzzling hardware behaviors more than abstracted analytical models, and we just use their outputs in our model.

V. Experiments

In this section we present several memory-bound applications that we modeled in order to predict their performance. Compute-bound applications are not considered because memory accesses are overlapped with computation and thus do not need to be optimized much.

One possible use of our model is to select the best working set size to achieve best performance. Another would be to select the minimal number of threads to use for a given computation in order to reach some level of performance. In order to illustrate these two approaches in the next sections, we will present comparison between our predictions and real applications. We only selected a few graphs in order to show interesting or unexpected results.

A. A Simple Code: Dot Product

The dot product computation is a simple code, consisting in the load of 2 different chunks. For our experiment, we assume all chunks are written beforehand by thread 0, therefore all chunks are in state (M, {0}). We choose to initialize vectors this way because it shows first, that the initialization phase can be critical for further performance, second because that is what many users would do in the first place. The performance prediction is compared to measurements of the MKL library dotproduct.

For both arrays read by dot product and for all of their chunks, the transition taken in Automaton 2 is a L T from state (M, {0}). When there is more than one thread, this leads to state (S, T ), while there is no state change with only thread 0. Since there is no chunk store in the dot product, the timing function is the sum of the timing functions for these transitions. Figure 5 shows that our model is able to predict the behavior of dotproduct function calls from the MKL library.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 10
The figures show that when only using one thread the dot product computation is faster when data set fits in the L1 cache while it is faster in L3 when using several threads. Again, this comes from coherence overhead: with a single thread, there is no coherence to maintain, thus we get better performance with faster memory. However with more threads, coherence gets involved due to the initialization of the vectors and performance is degraded. But when data only fits in the last level of cache, which is shared between all cores, then no cache coherence is required anymore, and code achieves better performance.

However, if we are careful with vector initialization (i.e. computing threads initialize the chunks they read), the dotproduct kernel can exhibit super-linear speedup as show on Figure 6 This super-linear speedup is due to Yet, with very large chunks, the dotproduct kernel can have a poor speedup even with a careful data initialization as shown on Figure 7. It increases up to 4 with 6 threads and then stagnates. The predicted speed-up is close to the measured one. As seen in Section IV-B, parallel capacity misses appear when all threads use the same shared cache, leading to poor scalability. Even with 8 threads running this kernel, each of them still manipulates 8 MB which does not fit in local caches. This is an example showing that our model handles parallel capacity misses: when multiple threads compete for the memory blocks of a shared cache, our performance model is able to take into account the capacity misses that occur. Moreover, while this is not observed here, memory contention could be predicted likewise.

B. Comparing Iterative Methods

Iterative Krylov methods are crucial algorithmic patterns for many numerical simulations, in particular for solving Partial Differential Equations. They are all built from the same building blocks (BLAS-1 and BLAS-2 functions) but exhibit different behaviors in term of convergence (application dependent) and time per iteration (architecture dependent). We compare in the following the performance prediction and measurement of one iteration of different Krylov methods, solving the equation Ax = b where x is an unknown vector, b a constant vector and A a matrix. The methods considered are Conjugate Gradient Normal Residual (CGNR), Conjugate Residual Normal Error (CRNE) and Biconjugate Gradient Stabilized (BiCGSTAB) [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. The code we use is adapted from a Lattice QCD simulation code [START_REF] Barthou | Automated code generation for lattice quantum chromodynamics and beyond[END_REF] and uses OpenMP C code. The matrix A is a 8000 sparse matrix of complex in double precision, consisting in small blocks of dense matrices on the diagonal and subdiagonals. As we are focusing only on the time per iteration, the iterative methods are arbitrarily stopped after 200 iterations.

All these methods use different access patterns, different working sets per thread, and have possible coherency issues (with dot product in particular). The scalability of the time per iteration for BiCGSTAB method is shown in Figure 8. The performance increases up to a maximum The timings of different methods are compared, for different thread numbers. Figure 9 shows that the time prediction is accurate enough to predict how performance changes when parallelism increases, and how different methods compare.

VI. Related Work

Several methods are commonly used to optimize software by observing and predicting performance. One is to simulate the full hardware, for instance with cycle accurate simulators [START_REF] Hughes | Rsim: Simulating Shared-Memory Multiprocessors with ILP Processors[END_REF], [START_REF] Covington | The Efficient Simulation of Parallel Computer Systems[END_REF]. Such predictions are very precise and permit collection of large amount of performance metrics. However they are time consuming and require a deep knowledge of the hardware in order to implement all architecture features, including prefetchers or cache replacement policies, with enough precision to provide cycle accurate simulation. Developing simulation software is a long process for each newly supported platform,

Normalized time 0.0 0.5 1.0 1.5 2.0 B iC G S T A B 1 B iC G S T A B 2 B iC G S T A B 4 B iC G S T A B 8 C G N R 1 C G N R 2 C G N R 4 C G N R 8 C R N E 1 C R N E 2 C R N E 4 C R N E 8 C R 1 C R 2 C R 4 C R 8
Real algorithm Model and it highly depends on the hardly-available hardware documentation. However, such simulation is able to explain precisely what happens for a given code fragment.

Our approach hides this complexity in the benchmarks and tries to remain portable by using a memory model that matches most widely-available modern processors and coherence protocols.

Profiling can also be used in order to record performance metrics for HPC applications tuning. Profiling has the same drawbacks as simulation since it slows down application performance. Tools such as Valgrind or Cachegrind [START_REF] Nethercote | Valgrind: A Program Supervision Framework[END_REF] can present an overhead up to 100 times the normal program execution time. Our approach is only time consuming when running the benchmarks once for each platform. Again, these approaches are complementary to the technique we propose, since profiling and simulation are able to provide more details concerning the hardware mechanisms triggered during a code execution.

Hardware counters are tools for performance tuning. More elaborated tools have been developed in order to ease the use of hardware counters [START_REF] Mucci | PAPI: A Portable Interface to Hardware Performance Counters[END_REF], [START_REF] Jarp | Perfmon2: A leap forward in Performance Monitoring[END_REF]. The advantage of counters compared to simulators is that it is lightweight: There is no overhead aside from the library initialization. However counters are not enough to optimize software. Indeed once a bottleneck of the application is found (let say, too many TLB misses), one needs a way to link the information back to source code in order to tackle the problem. Also, as discussed in Section IV-B, the overhead of misses significantly varies with cache states.

Hackenberg et al. compared cache coherence of real world CPUs in [START_REF] Hackenberg | Comparing Cache Architectures and Coherency Protocols on x86-64 Multicore SMP Systems[END_REF]. They show that cache coherency and cache data states are to be taken into account when modeling memory hierarchy. Williams et al. proposed an ingenious work in modeling both memory and computation in order to predict best achievable performance of a given code depending on its arithmetic intensity [START_REF] Williams | Roofline: An Insightful Visual Performance Model for Multicore Architectures[END_REF]. Ilic et al. extended the model to support caches and data reuse [START_REF] Ilic | Cache-aware Roofline model: Upgrading the loft[END_REF]. Compute-bound applications are handled while we are not able to predict computation performance. However, our model is able to predict in a better way applications with heavy coherence traffic. This also allows us to point out that bad performance of some applications can come from a huge overhead due to cache coherence. The references confirm the relevance of our approach for modeling memory access performance.

Our approach differs from all existing ones as the memory model is based on benchmarks. Benchmarks, especially the ones focusing on memory, have been developed in order to understand memory or application performance [START_REF] Treibig | Performance patterns and hardware metrics on modern multicore processors: Best practices for performance engineering[END_REF], [START_REF] Mccalpin | STREAM: Sustainable Memory Bandwidth in High Performance Computers[END_REF]. They are a great way to understand architecture behavior, however they can not be directly used to optimize software. Once our model is built for a given architecture, we are able to predict both software scalability and achievable memory bandwidth. By understanding the memory model or predicted scalability, one can see if performance is limited by memory contention or because of a cache coherence unfriendly memory access pattern. Also, our contribution is based on a model, there is no need to run the targeted application to predict its performance.

VII. Conclusion and Future Work

Code simulation and performance prediction become critical for performance analysis and software tuning. This paper presented an innovative model that predicts the performance of memory-bound OpenMP applications by composing the output of micro-benchmarks based on the state of data buffers in hardware caches. The model successfully predicts the execution time for several commonlyused application patterns, in particular those based on BLAS, and the effects of coherency and contention on performance. Additionally, our model incorporates the effects of compulsory and, to some extent, capacity misses. The systematic use of micro-benchmarks makes the model accurate without paying the cost to model complex hardware mechanisms such as prefetchers and cache coherency protocol implementations.

For future work, we plan to add support for multiple processor sockets. Loss of performance resulting from inter-socket interactions are ignored in our current model. While benchmarks are able to measure threads binded on different sockets, we need to plug them into the model, avoiding the combinatorial explosion coming from the exploration of the different bindings. Besides, we are thinking of adding automatic ways to detect coherence issues, for instance by defining some metrics based on the model and looking at applications feedback with hardware counters.

Figure 2 :

 2 Figure 2: Automaton for tracking chunk's states and for selecting the benchmarks that represent each memory access step.

1 : 4 forall the t = 1 ..T do 5 c 6 if Si is read then 7 θ 9 else 10 θ

 141567910 Now, the aggregation time for reads (or writes) in sequential phases of the program simply corresponds to the sum of the individual time for read/write statements. Algorithm 1 shows how the time is computed and states are updated. L and S are the functions obtained from the micro-benchmarks. To take into account the interleaving Algorithm Compute the time θ of read/write statements in a parallel region. Input: S1 . . . Sn // Sequence of statements in a parallel region Input: L, S // Functions defining transition timings 1 forall the t = 1..T do // Compute working set of thread t 2 s(t) ← n i=1 chunk(Si, t).size 3 for i = 1..n do ← chunk(Si, t) // Chunk accessed by thread t, Si L (t)+ = L(c.state, c.stride, s(t), T ) * c.size/s(t) 8 c.state ← δ(c.state, Lt) R (t)+ = S(c.state, c.stride, s(t), T ) * c.size/s(t) 11 c.state ← δ(c.state, St)

Figure 3 :

 3 Figure 3: Load Miss benchmarks.

Figure 4 :

 4 Figure 4: Store Hit bandwidth degradation for 8 threads.

Figure 5 :

 5 Figure 5: Dotproduct performance for 1 and 4 threads.

Figure 7 :

 7 Figure 7: Dotproduct strong scalability on 32 MB vectors.Yet, with very large chunks, the dotproduct kernel can have a poor speedup even with a careful data initialization as shown on Figure7. It increases up to 4 with 6 threads and then stagnates. The predicted speed-up is close to the measured one. As seen in Section IV-B, parallel capacity misses appear when all threads use the same shared cache, leading to poor scalability. Even with 8 threads running this kernel, each of them still manipulates 8 MB which does not fit in local caches. This is an example showing that our model handles parallel capacity misses: when multiple threads compete for the memory blocks of a shared cache, our performance model is able to take into account the capacity misses that occur. Moreover, while this is not observed here, memory contention could be predicted likewise.

Figure 8 :

 8 Figure 8: Speedup of the Bi-Conjugate gradient, Stabilized (BiCGSTAB).of 5 threads according to our model (6 according to experimental results).

Figure 9 :

 9 Figure 9: Performance for one iteration of different iterative methods. Performance figures are normalized, 1 corresponds to the time for the best sequential method. Two bars are given for Method n, corresponding to the normalized time when using n threads: in black, measured time, in red the predicted time.

  {t} This data chunk is in state Modified in the cache of the thread t. E {t} This chunk is only in the cache of the t th thread, in Exclusive state.

	S Ω	The chunk is in Shared state for all threads in Ω.
	I	This chunk is not present in any cache. At the
		beginning of the program, all chunks are in the
		Invalid state.

Therefore, for any array X and any mapping function f associated to X, the state of a chunk X[f (t, N, T )], for any t, N, T , is updated according to the accesses to it. This state will be denoted X[f (t, N, T )].state.


A Rodler

; J Virgone

J-J Roux

J L Hubert

IMPACT OF THE SUN PATCH ON HEATING AND COOLING POWER EVALUATION: APPLIED TO A LOW ENERGY CELL

Keywords: sun patch, fast climatic variations, heating and cooling power, insulated cell

. In this paper we present a 3D envelope model taking into account the minute-wise sun patch evolution. Simulation results are analysed for a low energy cell.

A numerical model has been developed in order to simulate the transient thermal behaviour, with a refined spatial (3D) and temporal (down to one minute time step) discretization of the single room. For each node of the grid, the energy conservation equations are developed. They traduce balance between short-wave and long-wave irradiations, convection, air enthalpy and three-dimensional heat conduction. The particularities of the program are that it projects the sun patch on the inner walls, the conduction is treated in three dimensions and climatic minute-wise variations are taken into account.

As main results, the surface temperature evolution, the air temperature evolution and the heating or cooling power necessary to maintain an inner air set-point temperature are calculated at each time-step. Heating or cooling power is compared to the power calculated with no sun patch incorporation (solar loads only on the floor). Conclusions are made on the importance of the integration of the sun patch and its impact on the observed results.

INTRODUCTION

Electrical energy savings has become an important issue. Among the consumption sources, the building's sector represents 40% of the consumed energy in France. The heating power represents 2/3 of the total power consumption of the buildings sector [START_REF]Stratégie utilisation rationnelle de l'énergie, chapitre II : « Les bâtiments[END_REF] so that, in January 2009 for example, heating power represented a third of the total maximal power demand [START_REF]La pointe d'électricité en France[END_REF] (Figure 1). Therefore, a particular attention is given to the development of low energy buildings, which are strongly insulated and require less energy demand.

Very often, the overall power charged by the supplier is done according to simplified models, which do not evaluate the power peaks in time and do not integrate fast climatic variations in time [START_REF] Weinmann | Prévision et justification des consommations d'électricité pour 3 catégories de bâtiment[END_REF]. On Figure 2, we can observe power discrepancies: the power charged by the supplier compared to the power measured does not overlap. The power charged by the supplier does not evaluate variations in time of the power consumption.

In dynamic simulations, the predominant method used to estimate the heat load is done by using the heat balance equation of the zone [START_REF]Multizone Building modeling with Type56 and TRNBuild[END_REF][START_REF] Energyplus | Energy Plus Engineering Reference[END_REF][START_REF]CODYBA[END_REF]. The heat is supplied to the air node by convection and/or radiation. Then, the heating load is estimated in order to reach the set-point temperature. When the temperature of the zone is higher than its set-point, the value of the load is set to zero. This method is coupled to the transient behaviour of the zone. In TrnSys for example, in order to describe a building zone, the walls have to be defined. They are considered as "black boxes" and the "thermal history" of the walls is defined using the transfer function method [START_REF] Mitalas | Calculation of Heat Conduction Transfer Functions for Multi-Layer Slabs[END_REF]. Using TrnSys for simulations of high thermal inertia walls, when they reach a certain value of thickness, can lead to important problems [START_REF] Savoyat | Integration of thick wall in TRNSYS simulation[END_REF][START_REF] Flory-Celini | Modélisation et positionnement de solutions bioclimatiques dans le bâtiment résidentiel existant[END_REF][START_REF] Delcroix | Conduction Transfer Functions in TRNSYS multizone building model: Current implementation, limitations and possible improvements[END_REF]. Also, heat conduction transfer through a wall is evaluated in one dimension, only. The aim of the actual work is to better evaluate the power peaks for low energy buildings. In order to do this, it is important to better model the envelop so that it takes into account more accurately the weather variations in time. In order to consider weather variations in time we use one minute time-step weather data and the sun patch is located for each minute. In order to locate the sun patch, we have discretized the walls in control volumes and used the 3D heat conduction equation in order to model more accurately the heat exchanges through the walls.

Then, the heating or cooling load is evaluated for given set-point temperatures.

METHOD

The building envelop is modeled following a three dimensional approach. The model represents a single room with a window.

The walls of the model are discretized. However, only a single air node is taken, considering that the air stratification has been neglected in this study in regard to the surface temperature distribution. The partial differential heat conduction equation in three dimensions for the temperature of the control volumes is applied to all the volumes of the cells. The heat conduction equation has been coupled to the convective and radiation exchanges.

= ( ) + -( ) - + ( ) + -( ) - + ( ) + -( ) - + Ф where (1)
Ф = Ф + Ф L + Ф C (2)
ФCONV is the convective heat flow between the surfaces and its environment and ФLW and ФSW are the long wave and short wave radiations. Thermal conductivities in the three directions are designated as x = y = z [W/mK] and the volumetric heat capacity is denoted by C [J/m 3 K]. A particular attention is given to the short wave radiations absorbed by the surfaces of a control volume, as it will depend on the sun patch location. The short wave radiations absorbed by the surfaces of a control volume come together in the vector {ФSWI} and are calculated by:

{Ф I } = [ ][ ]{ } (3)
where [S] is the surfaces matrix and [aSWI] is the matrix of absorptivity of the internal wall for the short wave radiations. {ESW} is the vector of radiations received by the meshes obtained resolving:

[ ]{ } = [ ]{ °} + [ ][][ ]{ } (4)
where [] is the reflectivity matrix, [FF] the view factors of the matrix which are calculated following the Nusselt analog. {ESW °} is the vector composed by primary radiations received by the control volumes. It results from the horizontal beam radiation Gb and the diffuse radiation Gd received by the meshes:

, ° = + , if the element i is in the sun patch , ° =
, if the element i is not included in the sun patch (5) τb and τd are the direct and diffuse transmission coefficients of the glass and depend on the incidence angle of the beam, whereas Rb and Rd allow to calculate direct and diffuse radiation on titled surfaces, with p the slope of the surface. The sun patch position has been calculated by a geometrical test: the boundary of the window is projected on an orthogonal plane to the beam. The control volumes of the walls are projected on the same plane, and thus those projected cells belonging to the projection of the window are identified. The energy balance equation for the air in the cavity with temperature is:

∁ = ∑ ∁ - = + ∑ ℎ - = (6)
where is the exterior dry bulb temperature (°C), is the air flow (kg/s), ∁ the heat capacity of the air (J/kgK), the air density (kg/m 3 ), Vc the volume of the cavity, the interior surface temperature and ℎ the convective transfer coefficient.

is the number of radiation balance equations corresponding to the number of surface mesh elements and the number of zones (here =1).

When wanting to maintain a set point temperature of 20°C for example, we maintain the set point temperature to 20°C and heat power is instantaneously injected on the air node. Differential equations 1 and 2 are solved with a function of Matlab (ode23t). This function solves moderately stiff ordinary differential equations using the trapezoidal rule or Runge Kutta. This method solves the equations with an adapted and variable time-step, depending on the fluctuations of the data.

CASE STUDIES AND HYPOTHESES

We have modelled a single cell, similar to an existing in situ experimental set up. Five walls among the six are strongly insulated (U = 0.099 W/m²K) and therefore heat losses through these walls are very small. Slightly higher losses take place through the south oriented wall, which has a window (Uwall = 0.277 W/m²K and Uwindow = 5.95 W/m²K). This room has a small thermal inertia and is highly insulated. The dimensions of the room are: 2.98 m in width, 2.89 m in depth and 2.83 m in height and the window's area measures 1.3×1.3m². The weather data is existing one minute time-step data from Vaulx-en-Velin of 2011 [START_REF]International Daylight Measurement Programme[END_REF]. The six walls and the window describing the cell contain 11000 control volumes, with a refined meshing near to the surfaces of the walls.

The first case study evaluates the cooling power necessary to maintain an inner set point temperature of 25°C, in July. Again, two models are studied: the first projects the incoming radiation on the floor (M1) and the second locates the sun patch throughout the day (M2).

The second case study consists in evaluating the power necessary to maintain an inner set point temperature of 20°C, in January. Two models are used: M1 and M2.

RESULTS OF DYNAMIC SIMULATIONS

For the first study, we evaluated cooling power for two days: 11 th and 12 th July 2011. Inside temperature was fixed to 25°C and cooling power was calculated for M1 and M2. For both models the incoming radiant flux is the same but it has been projected differently on the cells. On the 12 th July (15h UTC) we observe the inside temperature distribution for the two models M1 and M2. Figure 3a/ shows that the sun patch touches some of the cells of wall 3 and 2. In Figure 3b/ the incoming radiation has been projected on all the cells of the floor. We can observe that the floor cells' are warmer than the other cells of the room. This last model follows the traditional taken hypotheses by most of the simulation tools.

Figure 3: a/ with sun patch (M2) b/ incoming radiation projected on the floor (M1)

Despite the 25°C set point temperature, we notice differences on the temperature field (Figure 3) between both models. Maximal variations of 0.4°C are observed between the cells of the sun patch and the floor when comparing both models. This can explain the discrepancies between the two models when evaluating the power.

Cooling power is between 0 and 50W and gets more important when the temperature and the incident radiation is rises (Figure 4). Cooling power of M1 and M2 are compared. Discrepancies are observed, represented by the red area. Differences are larger when important solar fluctuations occur, during the day (red area). During this period the two models will project spatially the incoming radiation differently. With the sun patch model (M2) we have less cooling demand. For these two models and days of July, we have a relative difference of 83.9 %.

For the second case, we modelled the heating power necessary to maintain 20°C (Figure 5). Powers between 30 and 45 W are calculated. We observe discrepancies between these two models only during the day time, as for the other configuration. With the sun patch model, more heating is required during the day. A relative difference of 3.6 % is obtained between these two models.

CONCLUSIONS

A 3D transient envelope model has been developed. In this model, we discretized the walls in small control volumes, with a refined meshing towards the surfaces. This spatial discretization enables to know where the sun patch is located at each time step. The weather data was sampled at a minute time-step, in order to take into account the fast climatic variations in time. We applied this model to a low energy cell case, strongly insulated. A general similar trend between the results with the sun patch model and the results of S. Dautin [START_REF] Dautin | Réduction de modèles thermiques de bâtiments: amélioration des techniques par modélisation des sollicitations météorologiques[END_REF] was found, knowing that we have a 3D model and in her work only 1D exchange were modeled. In the observed results of this paper, it seems important to consider the sun patch position for days with clear sky conditions. Higher discrepancies are observed between the two models for the cooling power demand than for the heating power demand. Taking too simplified hypotheses, like done with model M1, for clear sky conditions, could lead to important errors. These results were found for an important insulated cell but future studies for different configurations of cells will be realized in order to draw broaden conclusions on the impact of the sun patch on power demand. These results have to be validated with experimental data, and therefore we are working on an in situ experimental set up.

Figure 1 :

 1 Figure 1: Heating power among the overall power consumption

Figure 4 :Figure 5 :

 45 Figure 4: Cooling demand for M1 and M2