Nejmeddine Bahri
email: nejmeddine.bahri@gmail.com

Abdessamad El Ansari

Mohamed Maazouz

Ali Ahaitouf

Nouri Masmoudi

HEVC Video encoder implementation On Texas Instruments Platforms

Keywords: HEVC encoder, DSP, ARM processor, RTOS, Cross-Compilation I

this paper presents a high efficiency video encoder implementation (HEVC) on two different Texas Instruments (TI) platforms: the BeagleBoard-xM based on ARM processor and the TMS320C6678 DSP. The new features of these processors such as multicore architecture, high frequency processor and low power consumption motivate researchers to develop an embedded HEVC video encoder which could be exploited in several multimedia applications such as High Definition (HD)TV, smart cameras, HD video conference, Ultra HD video surveillance systems. Different operating systems (OS) with different compilers are tested in order to obtain an optimized HEVC encoder implementation. Experimental results show that HEVC encoder DSP-based solution using SYS BIOS real time OS allows saving up to 62% of encoding time compared to Linux-c6x OS with the same DSP and about 71% of encoding time compared to BeagleBoard-based solution without inducing any performance degradation in terms of video quality or bitrate.

INTRODUCTION

Nowadays, the fast evolution of digital cameras technology makes users more consumers of UHD images and videos applications. Facing the huge amount of data transmission, the big memory storing requirement and the transmission bandwidth constraints, video encoding with high compression performance and low computational cost is absolutely required.

The High Efficiency Video Coding HEVC [START_REF] Bross | JCTVC-L1003_v30[END_REF] is the new generation of video coding standard. It is developed to overcome the huge amount of UHD video data. The HEVC encoder ensures a high encoding efficiency with respect to the already commercialized H264/AVC encoder. It saves up to 50% of bitrate with the same visual quality. This efficiency is accompanied by a tremendous computational complexity due to many new features included in this standard.

Develop an embedded solution for this encoder will be an interesting work to be incorporated in several multimedia devices. In fact, smart cameras, mobile phones, smart TV, MP4 players and surveillance systems usually include embedded processors to run video applications.

Currently, the majority of embedded implementations focus on the HEVC decoder which is less complicated and easier for implementation compared to the encoder one [START_REF] Chi | Parallel HEVC Decoding on Multi-and Many-core Architectures Journal of Signal Processing Systems[END_REF] [3] [START_REF] Pescador | Complexity analysis of an HEVC decoder based on a digital signal processor[END_REF] and [START_REF] Pescador | A DSP HEVC decoder implementation based on OpenHEVC[END_REF]. Regarding the HEVC encoder, there are two implementation aspects: some works focus on software optimizations of different encoder modules such as intra prediction [START_REF] Na | Edge-based fast mode decision algorithm for intra prediction in HEVC[END_REF], motion estimation [START_REF] Belghith | A new fast motion estimation algorithm using fast mode decision for highefficiency video coding standard Journal of Real-Time Image Processing[END_REF]. Others are interested in hardware implementations on different platforms such as Field Programmable Gate Array (FPGA) platform and x86 processors with Single Instruction Multiple Data (SIMD) architecture by using openMP library to accelerate encoding speed [START_REF] Miyazawa | Real-time hardware implementation of HEVC video encoder for 1080p HD video[END_REF] to [START_REF] Chen | Efficient SIMD optimization of HEVC encoder over X86 processors[END_REF].

In this paper, hardware implementations of HEVC encoder on two Texas Instruments embedded systems running at the same processor frequency are presented. BeagleBoard-xM platform based on ARM processor and the TMS320C6678 multicore DSP platform are chosen to run the HEVC encoder in order to evaluate and compare their performances. As preliminary work, our implementation is limited to a single core implementation. The motivation of this work is to perform an embedded encoder implementation that could be integrated on several devices such smart cameras, smart TV, UHD surveillance systems etc.

The remainder of this paper is organized as follows: section 2 gives a brief description of the target platforms features. Section 3 presents the HEVC encoder implementation strategies on the BeagleBoard-xM and DSP platforms using Linux operating system. Optimization steps by exploiting a real time operating system are detailed in section 4. Section 5 concludes the paper and presents some perspectives.

II. PLATFORMS DESCRIPTION

A. BeagleBoard-xM Platform

BeagleBoard-xM [START_REF]BeagleBoard-xM Rev C system Reference Manual[END_REF] is a low cost ARM Cortex A8 board. It is based on the DM3730 processor manufactured by Texas Instruments running at 1 GHz with extra memory @ 512MB of low-power DDR RAM as shown in Figure 1.

The BeagleBoard-xM is designed specifically to address the Open Source Community. It is characterized by low power consumption and integrates a large amount of On-Chip memory. It supports high-capacity microSD slot and 4-GByte microSD card which makes developers able to boot new operating systems as Angstrom Linux, Ubuntu, and Android without any dependencies on fixed non-volatile memory like NAND flash. Several interfaces are linked to this board such as 10/100 Ethernet for Internet Protocol networks and four high speed USB 2.0 ports used for USB peripherals. It includes a camera port to directly import videos from a connected camera. All these features make the BeagleBoard-xM very interesting solution to run several embedded multimedia applications like image and video processing applications, Robotics etc.

B. TMS320C6678 DSP Platform

DSP processors represent an attractive solution for embedded systems implementations and high performance applications thanks to software flexibility, low power consumption, competitive price tag and time to market. TMS320C6678 DSP [START_REF]TMS320C6678 Multicore Fixed and Floating-Point Digital Signal Processor Data Manual[END_REF] belongs to the latest generation of multicore DSPs made by Texas Instrument. It includes eight cores running each @ 1GHz which allows delivering 64000 MIPS performance. It is based on Very Long Instruction Word (VLIW) architecture and uses Single Instruction Multiple Data (SIMD) set instruction. Each C66x core integrates a large amount of on-chip memory. In addition to 32KB of L1 program and data cache, it includes 512KB of internal memory per core that can be configured as mapped RAM or cache. The platform also integrates 4MB of Multicore Shared Memory and 512MB of DDR3-1333 external memory. This device also includes 64MB of NAND Flash memory to run Linux-c6x Kernel. This platform comes with the TI's Multicore Software Development Kit (MCSDK) for SYS/BIOS Real Time Operating System (RTOS).

A. HEVC encoder implementation on BeagleBoard-xM

BeagleBoard-xM platform is based on DM3730 processor. This device supports high-level operating systems such as Windows CE, Linux, ubuntu and QNX etc. This platform comes with a 4GB microSD card that supports an Angstrom Linux validation image. To perform a functional HEVC video encoding application on the BeagleBoard-xM platform, we are processed as follows:

• A cross compile toolchain for ARM BeagleBoard-xM architecture is installed on the host machine in order to compile and build the HEVC encoder software.

• An Angstrom Linux kernel is installed in the SD card of the BeagleBoard-xM.

• The HEVC encoder is then compiled with the ARM cross compiler. An executable cross compiled binary output is generated. It is then copied with the input videos and the encoder configuration file ".cfg" into the SD card.

• The BeagleBoard-xM, which is connected with the host machine via the serial port RS232, boots on Linux from this SD card.

• minicom or picocom tool, which are a text-based serial port communications programs could be used to ensure communication between the host and the Beagleboard platform. Thus, execution of the binary executable can be started.

B. HEVC encoder implementation on TMS320C6678 DSP platform

As a preliminary work, only a mono-core HEVC encoder implementation is presented. For future work, a multicore implementation will be performed.

The HEVC software cannot be directly executed on the DSP using C6000 compiler with SYS/BIOS, the Real time operating system (RTOS) for DSP processor. The HEVC software includes several functions and data structures that are not compatible with DSP development frameworks. Thus, several changes should be performed in order to obtain a DSP functional version as Pescador's implementation [START_REF] Pescador | A DSP HEVC decoder implementation based on OpenHEVC[END_REF].

To avoid these modifications and easily obtain HEVC DSPbased solution, increasing the level of abstraction by exploiting a Linux gcc compiler instead of the standard C6000 SYS/BIOS compiler could be an efficient solution. In fact, TMS320C6678 includes a 64MB NAND Flash memory that can support a Linux operating system. So, the key idea is exploiting Linux instead of SYS/BIOS and a µclinux-c6x cross compiler instead of the standard C6000 compiler to compile and build the HEVC software.

The implementation steps of the HEVC encoder on the TMS320CC6678 DSP using a Linux OS are summarized as follows:

• Linux-c6x kernel, which is an embedded Linux dedicated for C6000 DSP family, is installed and loaded into the Nand Flash memory [START_REF]Linux-c6x-2.0-Beta2 release[END_REF].

• The HEVC encoder software is build using µclinux-c6x cross compiler.

• To ensure that the DSP exploits its entire memory in processing and loading Linux files system, the HEVC binary executable, the input video sequences and the encoder configuration file ".cfg" are kept on the host machine in a shared folder. The DSP will access to this shared folder with the Network File System (NFS) protocol [START_REF] Postel | File Transfer Protocol[END_REF].

• The EVM switches are configured according to the NAND flash boot mode.

• The DSP and the host machine are interconnected via both the Ethernet link using a static IP addresses and the serial RS232 interface to ensure between them.

• Another folder is created in the DSP file system. It will be used as the mount point. NFS command is then used to access to the shared folder on the host machine where the HEVC binary image and the video sequences have been saved. The DSP can now access to the HEVC executable binary code and the application could be executed without modifying any function or data structure or take into account data allocations into the DSP memory.

C. Experimental results with Linux Operating Systems

To evaluate our implementations performances, the HEVC encoder reference software (HM 12.1) [START_REF] Mccann | High Efficiency Video Coding (HEVC) Test Model 12 (HM 12) Encoder Description[END_REF] is executed on the both platforms that have the same processor frequency. In fact, BeagleBoard-xM ARM processor and the C6678 DSP core are both running at 1GHz. Experimental conditions are similar to common test conditions [START_REF] Bossen | Common test conditions and software reference configuration[END_REF] in HEVC standard. The largest coding unit (LCU) is fixed to 64x64 with a maximum depth level equal to 4. The entropy coder is CABAC (Context Adaptive Binary Arithmetic Coding). Different video sequences recommended by the JCT-VC [START_REF] Antonio | HEVC: A Review, Trends and Challenges[END_REF] are used for tests. The number of encoded frame for each video sequence is 100. All the simulations are evaluated with four QPs (Quantification Parameter) 22, 27, 32 and 37. As preliminary tests, only Class D (416x240) video resolution and "ALL INTRA" configuration are used to evaluate the performance of these implementations.

Table I presents the average encoding time (seconds) for one frame using Linux Operating System on BeagleBoard-xM and DSP TMS320C6678 platforms. Regarding quality degradation or bitrate increase, these two criteria have not changed since we have not made any modifications in the HEVC software. Only encoding time is taken into account to evaluate our implementations performance. Time saving is also computed as follows:

(%) = - ( 1
)
Experiments confirm that the HEVC video encoder represents a high computational complexity. An important execution time is required to encode one frame even for low video resolution. Experimental results show that our single core DSP implementation can save about 20% in average of encoding time compared to BeagleBoard-xM. This time saving may return to the C6678 external memory bandwidth where most of required encoding data are allocated compared to that of BeagleBoard-xM. In fact, C6678 includes a 64-bits DDR-3 external memory interface (EMIF) running at 1333 MHz whereas DDR's BeagleBoard is a 32 bits running at 166 MHz. This drastically affects the encoding time since the HEVC encoder includes an important number of Load/Store data instructions. Moreover, the C6678 includes a large amount of on-Chip memory (512KB of internal memory per core and 4MB of Multicore Shared Memory) characterized by high bandwidth which may reduce memory latency and speed up the encoding run-time. Based on these results and in order to speed up the encoding time, we have decided to keep working with the DSP platform and moving to use the SYS BIOS real time operating system with the C6000 compiler and the Integrated Development Environment (IDE) Code Composer Studio v5 [START_REF]Code Composer Studio (CCS) Integrated Development Environment (IDE)[END_REF].

IV. HEVC ENCODER IMPLEMENTATION USING REAL TIME OPERATING SYSTEM

A. HEVC implementation details on DSP platform

To obtain a DSP-based functional version using C6000 compiler, some modifications have been done in the C++ HEVC reference software in order to be compatible with DSP frameworks. The most important changes in the original code were:

• Use the string library instead of the memory library.

• Redefine the function "find" included in the class

TComIterator to adapt the input parameters.

• Include math library and use abs intrinsic function that calculates the absolute value of a variable.

• Redefine some C++ classes to adapt some variables unsupported by the DSP C++ compiler.

• Define the type "bool" in the encoder in order to be recognized by the C6000 compiler.

• Development of the function strdup, that is unavailable in the DSP development environment.

• Replace the functions used to measure the performance with PC timers by the appropriate functions that use DSP timers.

• Configure the compiler output format to ELF (Executable and Linkable Format) instead of COFF (Common Object File Format), needed to support some features of the C++ code.

• Check the compiler option "Support C ++ run-time type information -rtti" used to determine the type of a variable during the execution of the program.

B. Experimental results with SYS BIOS Real Time Operating System

By applying the changes indicated above, we have succeeded to obtain a DSP-based functional version with SYS BIOS operating system and C6000 compiler.

To compare between the two operating systems: Linux-C6x and SYS BIOS, the same test conditions and the same video sequences are used. Table II Experimental results show that the encoding time is improved by about 27% when using SYS BIOS and C6000 compiler compared to the implementation version using Linux-C6x and µclinux-c6x compiler. Hence, the average Class D encoding time is decreased from 17.6 seconds with Linux OS to 12.77 seconds with SYS BIOS. This enhancement returns to that C6000 compiler used with SYS BIOS outperforms the uClinux-C6x compiler used with Linux-C6x OS in terms of generating a well optimized assembly code that takes into account the internal architecture of the TMS320C6678 DSP platform.

C. System optimizations of HEVC implementation with SYS BIOS Real Time Operating System

Code Composer studio (CCS) environment gives the possibility to apply various optimizations either from build options project or from SYS BIOS configuration file based on the XDCtools [START_REF]XDCtools User's guide[END_REF]. This last one is a software product installed with CCS that contains all of the tools necessary to create, test, deploy, install, and use RTSC components [START_REF]Real-Time Software Components (RTSC)[END_REF] (Real Time Software component called packages) of SYS BIOS and to configure the application.

1) Configuration of program sections Map

During the first implementation of HEVC encoder on the C6678 DSP with SYS BIOS, all program memory sections have been allocated in the external DDR3 memory to avoid memory allocation problems since its size (512 MByte) is sufficient to support all data allocations of HEVC encoder. To optimize our implementation, we have decided to use the local internal L2RAM memory (512 KByte of size) and the shared internal MSMCSRAM memory (4 MByte of size) to allocate some program sections in order to get a faster data access compared to the external DDR3 memory.

Therefore, as shown in Table III, the ".text" section which presents the segment that contains the executable program code is configured in the shared MSMCSRAM memory (4 MByte) as its size is about 1 MByte. The "heap" segment which concerns dynamic allocation, is kept in the external DDR3 memory to support large memory allocations. Other sections, such as ".stack" which includes local variables, ".bss" for global variables, ".cinit" for initial values are allocated in the internal L2RAM memory since their sizes do not exceed 250 KByte.

2) Cache memory activation

In addition to program sections Map configuration, SYS BIOS gives the possibility to configure the target platform via the RTSC tool. Hence, to optimize our implementation, we have configured a part of L2RAM memory as a cache to speed up the data processing. In the first implementation, the whole L2 memory was configured as data memory. Since the L2RAM size is 512 KByte and knowing that the memory sections that have been allocated in don't exceed the 256 KByte, in this case the L2 cache size has been set to 256 KB. The L1D memory (Data memory Level 1) and L1P (Program memory Level 1) are always configured as 32Kbyte of cache even for the first implementation with SYS BIOS.

3) Project Build options optimizations

In order to generate a well optimized code level parallelism and software pipeline, the CCS environment provides the possibility to set the project build options that the C6000 compiler must take into account when compiling the code and generating the executable file. Consequently, the build options of our HEVC encoder project are fixed as follows:

• Build configuration mode: Release

• Optimization level (--opt_level,O): 3

• Optimize for code size (--opt_for_space,-ms) : 3

• Optimize for speed (--opt_for_speed, -mf): 5

• Inline function called once --single_inline: checked

D. Experimental results with system optimizations using the SYS BIOS real time operating system

The HEVC encoder is executed on the C6678 DSP after applying all the proposed system optimizations. The same test conditions are used in this implementation. Table IV shows a comparison between the experimental results of the nonoptimized version and the optimized one with SYS BIOS real time operating system. Experimental results show that the applied optimizations have ensured up to 48% of encoding time saving compared to the non-optimized version. In fact, the encoding time of a Class D image is decreased from 12.7 seconds to 6.52 seconds in average.

As result, our optimized implementation of the HEVC encoder on the C6678 DSP using SYS BIOS real time operating system allows saving up to 71% of encoding time compared to the implementation on the BeagleBoard-xM. In fact, the encoding time is decreased from 22.27 seconds on BeagleBoard-xM to 6.52 seconds on the C6678 DSP. Moreover, the use of SYS BIOS and C6000 compiler instead of Linux-c6x has saved up to 62% of encoding time which is decreased from 17.60 seconds with Linux-c6x to 6.52 seconds with SYS BIOS. Finally, we can consider this work as an important step to develop an embedded HEVC encoder despite the high complexity of this standard. This work will be the starting point for future researches. We have succeeded to obtain a functional HEVC encoder embedded on two different Texas Instruments processors without inducing any encoding degradation in terms of video quality or bitrate.

In order to improve the encoding time and reduce the computational complexity of the HEVC encoder, several optimizations could be proposed such as proposing fast coding algorithms and more optimized implementations exploiting the fully platform architecture. It is also possible to use pragmas instructions (MUST_ITERATE, UNROLL, INLINE etc) to optimize loops and functions. Moreover, a low-level assembly language programming of some functions such as SAD (Sum of Absolute Difference), SSE (Sum of Square Error), HAD (Transform HADAMARD) could also significantly improve the encoding speed as these functions are several times called for each Coding Unit (CU). Finally, exploiting the multicore architecture of the C6678 DSP and moving to a parallel implementation could be an interesting solution to accelerate the encoding process.

V. CONCLUSION

This paper proposed a HEVC video encoder implementation on two different platforms: BeagleBoard-xM based on ARM processor and the TMS320C6678 DSP, running both at the same processor frequency 1GHz. Different operating systems with different compilers have been tested in order to obtain an optimized implementation of this encoder. We have compared between the use of Linux operating systems and SYS BIOS real time operating system in terms of encoding time. Experimental results have shown that HEVC encoder DSP-based solution using SYS BIOS allows saving up to 62% of encoding time compared to the same DSP solution when using Linux-c6x operating system. Moreover, the DSP-based solution with SYS BIOS has saved about 71% of encoding time compared to BeagleBoard-based solution. Our implementations have not induced any performance degradation in terms of video quality or bitrate. Future work will be focused on the parallel implementation of the HEVC encoder on the multicore TMS3206678 DSP including eight C66x cores. This may give the possibility to reach real-time encoding when combined with architectural, software, and algorithmic optimizations.

Figure 1 .

 1 Figure 1. BeagleBoard-xM Platform

Figure 2 .

 2 Figure 2. TMS320C6678 DSP Platform III. HEVC ENCODER IMPLEMENTATION USING LINUX OPERATING SYSTEMS

TABLE I

 I

	.	AVERAGE ENCODING TIME FOR CLASS D (416X240) VIDEO
			SEQUENCES
	Class D (416x240) video sequences	QP	Average encoding time for one frame on BeagleBoard (s)	Average encoding time for one frame on C6678 DSP (s)	Time saving (%)
		22	23,07	18,23	20,98
	Basketball	27	20,24	16,09	20,50
	Pass	32	17,75	14,57	17,92
		37	15,90	13,35	16,04
		22	29,59	22,72	23,22
	BQSquare	27 32	26,15 22,72	20,05 17,97	23,33 20,91
		37	20,10	16,12	19,80
		22	28,77	22,25	22,66
	Blowing	27	24,79	19,18	22,63
	Bubbles	32	21,02	16,71	20,50
		37	17,88	14,62	18,23
		22	26,58	20,48	22,95
	Race	27	23,54	18,25	22,47
	Horses	32	20,28	16,71	17,60
		37	17,88	14,33	19,85
	Average time		22,27	17,60	20,60

  presents the average encoding time for class D (416x240) video sequences with different QP values.

TABLE II .

 II HEVC ENCODING TIME ON DSP PLATFORM USING LINUX-C6X AND SYS BIOS FOR CLASS D VIDEO SEQUENCES

	Class D (416x240) video sequences	QP	Average encoding time for one frame on C6678 using Linux-c6x (s)	Average encoding time for one frame on C6678 using SYS BIOS (s)	Time saving (%)
		22	18,23	13,70	24,85
	Basketball	27	16,09	11,83	26,48
	Pass	32	14,57	10,39	28,69
		37	13,35	9,32	30,19
		22	22,72	16,74	26,32
	BQSquare	27 32	20,05 17,97	14,61 12,71	27,13 29,27
		37	16,12	11,32	29,78
		22	22,25	16,53	25,71
	Blowing	27	19,18	14,03	26,85
	Bubbles	32	16,71	11,90	28,79
		37	14,62	10,27	29,75
		22	20,48	15,45	24,56
	Race	27	18,25	13,55	25,75
	Horses	32	16,71	11,77	29,56
		37	14,33	10,26	28,40
	Average time		17,60	12,77	27.43

TABLE III

 III

	.	PROGRAM SECTIONS MAP CONFIGURATION
	Program.sectMap ["heap"]	= "DDR3";
	Program.sectMap [".text"]	= "MSMCSRAM";
	Program.sectMap [".far"]	= "L2SRAM";
	Program.sectMap [".cinit"]	= "L2SRAM";
	Program.sectMap [".rodata"]	= "L2SRAM";
	Program.sectMap [".switch"]	= "L2SRAM";
	Program.sectMap [".args"]	= "L2SRAM";
	Program.sectMap [".bss"]	= "L2SRAM";
	Program.sectMap [".neardata"]	= "L2SRAM";
	Program.sectMap [".stack"]	= "L2SRAM";
	Program.sectMap [".data"]	= "L2SRAM";
	Program.sectMap [".cio"]	= "L2SRAM";
	Program.sectMap [".vecs"]	= "L2SRAM";

TABLE IV

 IV

	.	HEVC ENCODING PERFORMANCE ON C6678 DSP USING SYS BIOS WITH SYSTEM OPTIMIZATIONS
			Encoding time for	Encoding time for
	Class D (416x240) video sequences	QP	one frame on C6678 DSP using SYS BIOS without system	one frame on C6678 DSP using SYS BIOS with system	Time saving (%)
			optimizations (s)	optimizations (s)
		22	13,70	6,97	49,12
	Basketball	27	11,83	6,07	48,69
	Pass	32	10,39	5,35	48,51
		37	9,32	4,81	48,39
		22	16,74	8,46	49,46
	BQSquare	27 32	14,61 12,71	7,41 6,51	49,28 48,78
		37	11,32	5,80	48,76
		22	16,53	8,26	50,03
	Blowing	27	14,03	7,12	49,25
	Bubbles	32	11,90	6,10	48,74
		37	10,27	5,30	48,39
		22	15,45	7,86	49,13
	Race	27	13,55	6,95	48,71
	Horses	32	11,77	6,06	48,51
		37	10,26	5,30	48,34
	Average time		12,77	6,52	48,88


Emmanuel Frénod

L Émel

Lmam

Alexandre Mouton

Two-scale numerical simulation of the weakly compressible 1D isentropic Euler equations

Results

In this section, we claim the mathematical results leading to the construction of our twoscale numerical method.

After recalling existence and properties of solutions of (1.1), we give a first theorem giving the homogenized model and convergence properties. Then, we set out a numerical approximation for this homogenized model based on Roe's finite volume method which

Introduction

This paper enters a work programme concerning the development of two-scale numerical methods to solve PDEs with oscillatory singular perturbations linked to physical phenomena. Recently, these methods have been tested on simple problems. For instance, in Ailliot, Frénod and Monbet [START_REF] Ailliot | Long term object drift in the ocean with tide and wind[END_REF], such a method is used to manage the tide oscillation for long term drift forecast of objects in the ocean; in Frénod, Salvarani and Sonnendrücker [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale PIC-method[END_REF], such a method is used to simulate a charged particle beam in a periodic focusing channel.

The ultimate goal of this work programme is to propose efficient numerical methods to simulate plasmas submitted to strong magnetic field. Of course, simulations of magnetic confinement fusion are part of this ultimate goal.

Before going further towards this ultimate goal, several questions concerning the behaviour of the concerned methods in front of non-linearities and non-smooth solutions need to be investigated. The field of weakly compressible 1D isentropic Euler equations offers a relatively confortable framework in order to tackle these questions.

Indeed, the mathematical context established in Grenier [START_REF] Grenier | Oscillatory pertubation of the Navier-Stokes equations[END_REF], Klainerman and Majda [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of the compressible fluids[END_REF][START_REF] Klainerman | Compressible and imcompressible fluids[END_REF], Fortenbach, Frénod, Munz and Sonnendrücker [START_REF] Fortenbach | Multiple scale consideration for sound generation in low Mach number flow[END_REF], Frénod, Raviart and Sonnendrücker [START_REF] Frénod | Two-scale expansion of a singularly perturbed convection equation[END_REF], Frénod and Sonnendrücker [START_REF] Frénod | Long time behavior of the Vlasov equation with a strong external magnetic field[END_REF][START_REF] Frénod | The finite Larmor radius approximation[END_REF], Majda [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF], Métivier and Schochet [START_REF] Métivier | The incompressible limit of the non-isotropic Euler equations[END_REF], Munz [START_REF] Munz | Computational fluid dynamics and aeroacoustics for low Mach number flow[END_REF], Schochet [START_REF] Schochet | Asymptotic for symmetric hyperbolic systems with a large parameter[END_REF][START_REF] Schochet | Fast singular limit of hyperbolic pdes[END_REF][START_REF] Schochet | Symmetric hyperbolic systems with a large parameter[END_REF][START_REF] Schochet | The compressible Euler equations in a bounded domain existence of solutions and the incompressible limit[END_REF], offers a nice marked out way. This is the first motivation of the present paper.

A second motivation originates from the fact that robust numerical methods such as finite volumes set out in Godunov [START_REF] Godunov | A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations[END_REF], LeVeque [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] and Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] have been developped in order to solve compressible or incompressible Euler equations. Nevertheless, when they are used to simulate the flow in an experiment with conditions inducing a small Mach number, the CPU time cost of these methods is too high for operational applications because of the very small time step required in order to capture high speed travelling waves that appear in this case.

The precise aim of this paper is to develop a two-scale numerical method to solve

∂ t u ǫ + 1 2 ∂ x (u ǫ ) 2 + ∂ x q ǫ (ρ ǫ ) + 1 ǫ ∂ x ρ ǫ = 0 , ∂ t ρ ǫ + ∂ x (ρ ǫ u ǫ ) + 1 ǫ ∂ x u ǫ = 0 , u ǫ |t = 0 = u 0 , ρ ǫ |t = 0 = ρ 0 , (1.1)
where (x, t) ∈ T 1 × [0, T ). In this model, ǫ is the Mach number, u ǫ = u ǫ (x, t) is the dimensionless speed of the fluid, ρ ǫ = ρ ǫ (x, t) is linked to the dimensionless density ρǫ = ρǫ (x, t) by the relation

ρ ǫ = ρǫ -1 ǫ , (1.2)
and

q ǫ (ρ ǫ ) = γ -2 2 (ρ ǫ ) 2 + ǫ q ǫ 1 (ρ ǫ ) , (1.3)
where γ is the adiabatic coefficient of the fluid and q ǫ 1 (ρ ǫ ) is regular. This form of the isentropic Euler equations is set out in Grenier [START_REF] Grenier | Oscillatory pertubation of the Navier-Stokes equations[END_REF] (p. 494) and Métivier and Schochet [START_REF] Métivier | The incompressible limit of the non-isotropic Euler equations[END_REF]. Equation (1.1) is obtained from the classical form of the 1D isentropic Euler equations

∂ t u ǫ + u ǫ ∂ x u ǫ + 1 ǫ 2 ρǫ ∂ x p(ρ ǫ ) = 0 , ∂ t ρǫ + ∂ x ρǫ u ǫ = 0 , u ǫ (x, 0) = u 0 (x) , ρǫ (x, 0) = 1 + ǫ ρ 0 (x) , (1.4)
where usually, pressure function is given by p(ρ ǫ ) = (ρ ǫ ) γ γ . Using the relation (1.2) between ρ ǫ and ρǫ , we obtain

1 ǫ 2 ρǫ ∂ x p(ρ ǫ ) = 1 ǫ 2 (1 + ǫ ρ ǫ ) ∂ x (1 + ǫ ρ ǫ ) γ γ = (1 + ǫ ρ ǫ ) γ-2 ǫ ∂ x ρ ǫ . (1.5)
Making and expansion of (1 + ǫ ρ ǫ ) γ-2 , we obtain

1 ǫ 2 ρǫ ∂ x p(ρ ǫ ) = 1 ǫ (∂ x ρ ǫ ) 1 + (γ -2) ǫ ρ ǫ + ǫ 2 q ǫ 1 (ρ ǫ ) , (1.6)
where q ǫ 1 (ρ ǫ ) is regular, and, by introducing this result in the first equation of (1.4), we deduce the first equation of (1.1). In order to obtain the second equation of (1.1), we introduce the relation (1.2) in the second equation of (1.4).

In order to achieve this precise aim, inspired by two-scale convergence theory developped in Nguetseng [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] and Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], we establish a homogenized model describing the mean behaviour of (u ǫ , ρ ǫ ). This homogenized model neither contains nor generates high frequency oscillations but only their mean action. We also establish a way to reconstruct the oscillating solution (u ǫ , ρ ǫ ) from obtained mean behaviour. The homogenized model and the reconstruction procedure involve an additional variable which is a copy of the time variable and which allows the transfert of informations between the mean behaviour and the high frequency oscillations.

is the two-scale numerical method that allows to solve (1.1). Finally, we give a theorem giving some properties about the convergence of the considered numerical scheme.

It is well known (see Grenier [START_REF] Grenier | Oscillatory pertubation of the Navier-Stokes equations[END_REF], Klainerman and Majda [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of the compressible fluids[END_REF], Majda [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF], Metivier and Schochet [START_REF] Métivier | The incompressible limit of the non-isotropic Euler equations[END_REF], or Schochet [START_REF] Schochet | Asymptotic for symmetric hyperbolic systems with a large parameter[END_REF][START_REF] Schochet | Symmetric hyperbolic systems with a large parameter[END_REF][START_REF] Schochet | The compressible Euler equations in a bounded domain existence of solutions and the incompressible limit[END_REF]) that, with initial data satisfying

u 0 , ρ 0 ∈ H s (T 1 ) with s > 3 2 , (2.1)
there exists an existence time T > 0, independent of ǫ, such that the system (1.1) admits a unique solution (u

ǫ , ρ ǫ ) in C 0, T ; H s (T 1 ) 2 ∩ C 1 0, T ; H s-1 (T 1 ) 2 for any ǫ > 0. Fur- thermore, the sequence (u ǫ , ρ ǫ ) ǫ > 0 is bounded in C 0, T ; H s (T 1 ) 2 ∩ C 1 0, T ; H s-1 (T 1 )
2 for the usual norm.

Having this existence result at hand, we can claim the following theorem giving the homogenized model.

Theorem 1. If we fix ǫ 0 > 0 and s ≥ 2, there exists a constant K > 0 such that, for any ǫ ∈ ]0, ǫ 0 ], we have

T 0 u ǫ (•, t) -U •, t ǫ , t 2 L 2 (T 1 ) dt 1 2 ≤ Kǫ , T 0 ρ ǫ (•, t) -R •, t ǫ , t 2 L 2 (T 1 ) dt 1 2 ≤ Kǫ , (2.2)
where U and R are defined by

U (x, τ, t) = F (x -τ, t) + B(x + τ, t) + u 2π , R(x, τ, t) = F (x -τ, t) -B(x + τ, t) + ρ 2π , (2.3
)

with u = T 1 u 0 (x) dx, ρ = T 1 ρ 0 (x) dx, and F, B ∈ C 0, T ; H s (T 1 ) satisfying T 1 F (x, t) dx = T 1 B(x, t) dx = 0 ∀ t ∈ [0, T ) , (2.4)
and

∂ t F + ∂ x γ + 1 4 F 2 + 2u + (γ -1)ρ 4π F = 0 , ∂ t B + ∂ x γ + 1 4 B 2 + 2u -(γ -1)ρ 4π B = 0 , F | t = 0 = 1 2 u 0 + ρ 0 - u + ρ 2π , B | t = 0 = 1 2 u 0 -ρ 0 - u -ρ 2π .
(2.5)

A way to interpret this theorem is that if we can compute some good approximations F h and B h of the solutions F and B of (2.5), we are able to reconstruct some functions U h and R h in the same way as we reconstruct U and R in (2.3). Then the obtained functions (x, t) → U h (x, t ǫ , t) and (x, t) → R h (x, t ǫ , t) are good approximations of (x, t) → u ǫ (x, t) and (x, t) → ρ ǫ (x, t) respectively.

Based on this idea, we build our two-scale numerical method: firstly, we consider a uniform mesh on T 1 × [0, T ] with space step h and time step k, and we use the notations x i = ih and t n = nk. Then we build F h and B h such that

F h (x, t) = F n i and B h (x, t) = B n i ∀ (x, t) ∈ [x i-1/2 , x i+1/2 [×[t n , t n+1 [ , (2.6)
where F n i and B n i are computed with Roe's finite volume approximation of (2.5):

F n+1 i = F n i - k 2h γ + 1 4 (F n i+1 ) 2 -(F n i-1 ) 2 + 2u + (γ -1)ρ 4π (F n i+1 -F n i-1 ) - γ + 1 4 (F n i+1 + F n i ) + 2u + (γ -1)ρ 4π (F n i+1 -F n i ) + γ + 1 4 (F n i + F n i-1 ) + 2u + (γ -1)ρ 4π (F n i -F n i-1 ) , (2.7)
B n+1 i = B n i - k 2h γ + 1 4 (B n i+1 ) 2 -(B n i-1 ) 2 + 2u -(γ -1)ρ 4π (B n i+1 -B n i-1 ) - γ + 1 4 (B n i+1 + B n i ) + 2u -(γ -1)ρ 4π (B n i+1 -B n i ) + γ + 1 4 (B n i + B n i-1 ) + 2u -(γ -1)ρ 4π (B n i -B n i-1
) .

(2.8)

Finally, we compute U h and R h by

U h (x, τ, t) = F h (x -τ, t) + B h (x + τ, t) + u 2π , R h (x, τ, t) = F h (x -τ, t) -B h (x + τ, t) + ρ 2π , (2.9)
and we have the following convergence result.

Theorem 2. If s ≥ 2, then the approximations U h (•, τ, •) and R h (•, τ, •) converge to U (•, τ, •) and R(•, τ, •) in L 1 [0, T ) × T 1 for any τ ∈ T 1 . Furthermore, if s ≥ 3,
u ǫ -U (•, • ǫ , •) ⇀ 0 weakly - * , ρ ǫ -R(•, • ǫ , •) ⇀ 0 weakly - * , (2.10)
in L ∞ 0, T ; H s (T 1 ) when ǫ → 0, with U and R solution of (2.3)-(2.5).

3 Construction of the two-scale numerical method

Homogenized model

We present here the construction of the model (2.5). We recall that we start from the weakly compressible 1D isentropic Euler equations (1.1)- (1.3). Introducing the functions f ǫ , b ǫ defined by

f ǫ (x, t) = 1 2 u ǫ x + t ǫ , t + ρ ǫ x + t ǫ , t - u + ρ 2π , b ǫ (x, t) = 1 2 u ǫ x + t ǫ , t -ρ ǫ x + t ǫ , t - u -ρ 2π , (3.1)
we can rewrite the system (1.1):

∂ t f ǫ (x, t) + ǫ 2 ∂ x q ǫ 1 f ǫ (x, t) -b ǫ x + 2t ǫ , t + ρ 2π + ∂ x γ + 1 4 f ǫ (x, t) 2 + γ -3 4 b ǫ x + 2t ǫ , t 2 + 2u + (γ -1)ρ 4π f ǫ (x, t) + (3 -γ)ρ 4π b ǫ x + 2t ǫ , t + 3 -γ 2 f ǫ (x, t)b ǫ x + 2t ǫ , t = 0 , (3.2)
∂ t b ǫ (x, t) + ǫ 2 ∂ x q ǫ 1 f ǫ x - 2t ǫ , t -b ǫ (x, t) + ρ 2π + ∂ x γ + 1 4 b ǫ (x, t) 2 + γ -3 4 f ǫ x - 2t ǫ , t 2 + 2u -(γ -1)ρ 4π b ǫ (x, t) - (3 -γ)ρ 4π f ǫ x - 2t ǫ , t + 3 -γ 2 f ǫ x - 2t ǫ , t b ǫ (x, t) = 0 , (3.3)
equipped with

f ǫ | t = 0 = f 0 = 1 2 u 0 + ρ 0 - u + ρ 2π , (3.4)
b ǫ | t = 0 = b 0 = 1 2 u 0 -ρ 0 - u -ρ 2π . (3.5)
Since the sequences (f ǫ ) ǫ > 0 and (b ǫ ) ǫ > 0 are bounded in L ∞ 0, T ; H s (T 1 ) , there exist two functions F and B in L ∞ 0, T ; H s (T 1 ) such that for some subsequences always denoted (f ǫ ) ǫ > 0 and (b ǫ ) ǫ > 0 , we have

f ǫ ⇀ F weakly-* in L ∞ 0, T ; H s (T 1 ) , (3.6) b ǫ ⇀ B weakly-* in L ∞ 0, T ; H s (T 1 ) , (3.7)
for ǫ → 0.

A first property of F and B is that their averages are equal to 0: to show this, we integrate (3.2)-(3.3) with respect to x to obtain

∂ t T 1 f ǫ (x, t) dx = ∂ t T 1 b ǫ (x, t) dx = 0 .
Then we deduce that, for all t ∈ [0, T ), we have

T 1 f ǫ dx = T 1 f 0 dx = 0 , T 1 b ǫ dx = T 1 b 0 dx = 0 . (3.8)
These results, combined with the fact that f ǫ and b ǫ weakly-* converge to F and B respectively, lead to the results (2.4).

Furthermore, usual compactness results (for example Lions [START_REF] Lions | Quelques méthodes de résolution de problèmes aux limites non linéaires[END_REF]) yield that the functional space

U = g ∈ L ∞ 0, T ; H s (T 1 ) : ∂ t g ∈ L ∞ 0, T ; H s-1 (T 1 ) , (3.9)
provided with the usual product norm, is compactly embedded in L ∞ 0, T ; H s-1 (T 1 ) . As a consequence, since (f ǫ ) ǫ > 0 and (b ǫ ) ǫ > 0 are bounded in U, we have

f ǫ ǫ → 0 -→ F strongly in L ∞ 0, T ; H s-1 (T 1 ) , (3.10)
b ǫ ǫ → 0 -→ B strongly in L ∞ 0, T ; H s-1 (T 1 ) .

(3.11)

Since f ǫ and b ǫ are continuous in t, we can remark that F, B ∈ C 0, T ; H s-1 (T 1 ) , so the convergence result (3.10)- (3.11) is also true in C 0, T ; H s-1 (T 1 ) . Having now this convergence result at hand, we can look for the constraint equations on F and B. For that, we multiply the equations (3.2)-(3.3) by a regular function ϕ with compact support on T 1 × [0, T ) and we integrate with respect to x and t on T 1 × [0, T ). We obtain

- T 0 T 1 f ǫ (x, t) ∂ t ϕ(x, t) dx dt - T 1 f 0 (x) ϕ(x, 0) dx - T 0 T 1 γ + 1 4 f ǫ (x, t) 2 ∂ x ϕ(x, t) dx dt - T 0 T 1 2u + (γ -1)ρ 4π f ǫ (x, t) ∂ x ϕ(x, t) dx dt - T 0 T 1 γ -3 4 b ǫ x + 2t ǫ , t 2 ∂ x ϕ(x, t) dx dt - T 0 T 1 (3 -γ)ρ 4π b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt - T 0 T 1 3 -γ 2 f ǫ (x, t) b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt - ǫ 2 T 0 T 1 q ǫ 1 f ǫ (x, t) -b ǫ x + 2t ǫ , t + ρ 2π ∂ x ϕ(x, t) dx dt = 0 , (3.12)
and

- T 0 T 1 b ǫ (x, t) ∂ t ϕ(x, t) dx dt - T 1 b 0 (x) ϕ(x, 0) dx - T 0 T 1 γ + 1 4 b ǫ (x, t) 2 ∂ x ϕ(x, t) dx dt - T 0 T 1 2u -(γ -1)ρ 4π b ǫ (x, t) ∂ x ϕ(x, t) dx dt - T 0 T 1 γ -3 4 f ǫ x - 2t ǫ , t 2 ∂ x ϕ(x, t) dx dt + T 0 T 1 (3 -γ)ρ 4π f ǫ x - 2t ǫ , t ∂ x ϕ(x, t) dx dt - T 0 T 1 3 -γ 2 f ǫ x - 2t ǫ , t b ǫ (x, t) ∂ x ϕ(x, t) dx dt - ǫ 2 T 0 T 1 q ǫ 1 f ǫ x - 2t ǫ , t -b ǫ (x, t) + ρ 2π ∂ x ϕ(x, t) dx dt = 0 . (3.13)
Because of (3.10) and (3.11), passing to the limit in the four first terms of (3.12) and (3.13) is straightforward. Because of the factor ǫ 2 in front of the last term of (3.12) and (3.13), we deduce that these terms converge to 0. In order to find the limit of

T 0 T 1 b ǫ x+ 2t ǫ , t ∂ x ϕ(x, t) dx dt we firstly make a change of variables T 0 T 1 b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = T 0 T 1 b ǫ x, t) ∂ x ϕ x - 2t ǫ , t dx dt ,
then, we write

T 0 T 1 b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = T 0 T 1 b ǫ x, t) -B(x, t) ∂ x ϕ x - 2t ǫ , t dx dt + T 0 T 1 B x, t) ∂ x ϕ x - 2t ǫ , t dx dt .
Then, using the Cauchy-Schwarz inequality, we show that

lim ǫ → 0 T 0 T 1 b ǫ x, t) -B(x, t) ∂ x ϕ x - 2t ǫ , t dx dt = 0 . (3.14) If we define ϕ ǫ (x, t) = φ(x, t ǫ , t) where φ(x, τ, t) = ϕ(x -2τ, t), the sequence (ϕ ǫ ) ǫ > 0 is bounded in L ∞ 0, T ; H s (T 1 )). Then, there exists a function Φ in L ∞ 0, T ; L 2 (T 1 ; H s (T 1 )) such that lim ǫ → 0 T 0 T 1 ψ x, t ǫ , t ϕ ǫ (x, t) dx dt = T 0 T 1 T 1 ψ(x, τ, t) Φ(x, τ, t) dτ dx dt (3.15)
for any regular function ψ defined on T 1 × T 1 × [0, T ). As described in Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], Φ is called the two-scale limit of ϕ ǫ . Furthermore, following Marusic-Paloka and Piatnitski [START_REF] Marusic-Paloka | Homogenization of nonlinear convection-diffusion equation with rapidly oscillating coefficients and strong convection[END_REF] and Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], since φ is regular and 2π-periodic in τ , it is an easy game to prove that

Φ = φ in L ∞ 0, T ; L 2 (T 1 ; H s (T 1 )) . As a consequence, we have Φ(x, τ, t) = ϕ(x -2τ, t) and ϕ ǫ ⇀ T 1 ϕ(x -2τ, t) dτ weakly-* in L ∞ 0, T ; H s (T 1 ) . (3.16)
Using these results, we obtain lim

ǫ → 0 T 0 T 1 B(x, t)∂ x ϕ x - 2t ǫ , t dx dt = - T 0 T 1 T 1 ∂ x B(x, t) ϕ(x -2τ, t) dτ dx dt = - 1 2 T 0 T 1 B(x, t) T 1 ∂ τ ϕ(x -2τ, t) dτ dx dt = 0 , (3.17)
and finally

lim ǫ → 0 T 0 T 1 b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = 0 . (3.18)
In the same way, we obtain

lim ǫ → 0 T 0 T 1 f ǫ x - 2t ǫ , t ∂ x ϕ(x, t) dx dt = 0 , (3.19)
lim ǫ → 0 T 0 T 1 b ǫ x + 2t ǫ , t 2 ∂ x ϕ(x, t) dx dt = 0 , (3.20)
lim ǫ → 0 T 0 T 1 f ǫ x - 2t ǫ , t 2 ∂ x ϕ(x, t) dx dt = 0 . (3.21)
To find the limit of

T 0 T 1 f ǫ (x, t) b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt, we notice that T 0 T 1 f ǫ (x, t) b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = T 0 T 1 f ǫ (x, t) -F (x, t) B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt + T 0 T 1 F (x, t) b ǫ x + 2t ǫ , t -B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt + T 0 T 1 f ǫ (x, t) -F (x, t) b ǫ x + 2t ǫ , t -B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt + T 0 T 1 F (x, t) B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt .
With Cauchy-Schwarz and Hölder's inequalities, we easily show that

T 0 T 1 f ǫ (x, t) -F (x, t) B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt → 0 , (3.22) T 0 T 1 F (x, t) b ǫ x + 2t ǫ , t -B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt → 0 , (3.23) T 0 T 1 f ǫ (x, t) -F (x, t) b ǫ x + 2t ǫ , t -B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt → 0 , (3.24)
when ǫ → 0. If we define the functions B ǫ (x, t) = B x + 2t ǫ , t , we remark that the sequence (B ǫ ) ǫ > 0 is bounded in L ∞ 0, T ; H s-1 (T 1 ) and then admits a two

-scale limit B in L ∞ 0, T ; L ∞ (T 1 ; H s-1 (T 1 )) satisfying B(x, τ, t) = B(x + 2τ, t) ∀ (x, τ, t) ∈ T 1 × T 1 × [0, T ) , (3.25)
B ǫ ⇀ T 1 B(x + 2τ, t) dτ = 0 weakly-* in L ∞ 0, T ; H s (T 1 ) . (3.26)
This means that, for any function ψ ∈ C 0, T ; H s-1 (T 1 ) , we have

lim ǫ → 0 T 0 T 1 B ǫ (x, t) ψ(x, t) dx dt = 0 . (3.27) Setting then ψ = F ∂ x ϕ in (3.27), we obtain lim ǫ → 0 T 0 T 1 F (x, t) B x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = 0 . (3.28) Hence lim ǫ → 0 T 0 T 1 f ǫ (x, t) b ǫ x + 2t ǫ , t ∂ x ϕ(x, t) dx dt = 0 . (3.29)
The same method also gives

lim ǫ → 0 T 0 T 1 b ǫ (x, t) f ǫ x - 2t ǫ , t ∂ x ϕ(x, t) dx dt = 0 . (3.30)
With these results, we can pass to the limit in (3.12)-(3.13) and get

- T 0 T 1 F (x, t) ∂ t ϕ(x, t) dx dt - T 1 f 0 (x) ϕ(x, 0) dx - T 0 T 1 γ + 1 4 F (x, t) 2 + 2u + (γ -1)ρ 4π F (x, t) ∂ x ϕ(x, t) dx dt = 0 , (3.31)
and

- T 0 T 1 B(x, t) ∂ t ϕ(x, t) dx dt - T 1 b 0 (x) ϕ(x, 0) dx - T 0 T 1 γ + 1 4 B(x, t) 2 + 2u -(γ -1)ρ 4π B(x, t) ∂ x ϕ(x, t) dx dt = 0 . (3.32)
Remembering the definition (3.4)-(3.5) of f 0 and b 0 , we recognize here the weak formulation of (2.5).

Under the hypothesis of Theorem 1 about the initial data in (2.5), we deduce that the solution (F, B) belongs to C 0, T ; H s (T 1 )

2 and is unique in this space (we can apply the Theorem 3.6.1 of Serre [START_REF] Serre | Systems of conservations laws 1[END_REF] for example). Finally, from uniqueness, we have all the convergence results above for the whole sequences (f ǫ ) ǫ > 0 , (b ǫ ) ǫ > 0 and not only for some subsequences.

Convergence of the homogenized model: proof of Theorem 1

In this section, we will prove inequalities (2.2). Before going further, we remark that convergence (3.6)-(3.7) is already proved and that, thanks to (2.3) and (3.1), these inequalities are equivalent to the theorem below.

Theorem 4. For ǫ 0 > 0 and s ≥ 2 fixed, there exists a constant K > 0 such that

T 0 f ǫ (•, t) -F (•, t) 2 L 2 (T 1 ) dt 1 2 ≤ Kǫ , (3.33) T 0 b ǫ (•, t) -B(•, t) 2 L 2 (T 1 ) dt 1 2
≤ Kǫ , (3.34)

for any ǫ ∈ ]0, ǫ 0 ].
Proof of theorem 4: to prove this theorem, in the same spirit of Frénod, Raviart and Sonnendrücker [START_REF] Frénod | Two-scale expansion of a singularly perturbed convection equation[END_REF], we introduce the functions γ ǫ and δ ǫ defined by

γ ǫ (x, t) = 1 ǫ f ǫ (x, t) -F (x, t) -W x, t ǫ , t , δ ǫ (x, t) = 1 ǫ B(x, t) -b ǫ (x, t) -V x, t ǫ , t , (3.35)
with W and V defined on

T 1 × T 1 × [0, T ) by W (x, τ, t) = - γ -3 4 τ 0 ∂ x B(x + 2θ, t) 2 -2F (x, t)B(x + 2θ, t) - ρ π B(x + 2θ, t) dθ , V (x, τ, t) = γ -3 4 τ 0 ∂ x F (x -2θ, t) 2 -2F (x -2θ, t)B(x, t) + ρ π F (x -2θ, t) dθ .
(3.36)

Concerning sequences (γ ǫ ) ǫ ∈ ]0,ǫ 0 ] and (δ ǫ ) ǫ ∈ ]0,ǫ 0 ] , we have the property below.

Lemma 5. Under the hypothesis of theorem 4, there exists a contant M > 0, independent of ǫ ∈ ]0, ǫ 0 ] and t ∈ [0, T ), such that

∂ t γ ǫ (•, t) 2 L 2 (T 1 ) + δ ǫ (•, t) 2 L 2 (T 1 ) ≤ M 3 2 + √ 2 γ ǫ (•, t) 2 L 2 (T 1 ) + δ ǫ (•, t) 2 L 2 (T 1 ) + √ 2M , (3.37)
for any ǫ ∈ ]0, ǫ 0 ] and t ∈ [0, T ).

Having this inequality at hand, we can apply Gronwall's lemma to find a constant L > 0, independent of t and ǫ, such that

γ ǫ (•, t) 2 L 2 (T 1 ) + δ ǫ (•, t) 2 L 2 (T 1 ) ≤ L , ∀ (t, ǫ) ∈ [0, T )×]0, ǫ 0 ] . (3.38)
Finally, we define the constant K by

K = √ T √ L + max sup t ∈ [0,T ) τ ∈ T 1 W (•, τ, t) H s (T 1 ) , sup t ∈ [0,T ) τ ∈ T 1 V (•, τ, t) H s (T 1 ) , (3.39)
and we obtain the inequalities (3.33)-(3.34), giving the theorem. Now, to finish the proof of Theorem 4, we have to prove Lemma 5. Firstly, we notice that ∂ t γ ǫ and ∂ t δ ǫ read

∂ t γ ǫ (x, t) = -α f ǫ (x, t) + F (x, t) + ζB x + 2t ǫ , t -β + ∂ x γ ǫ (x, t) + -α∂ x f ǫ (x, t) + F (x, t) + ζ∂ x B x + 2t ǫ , t γ ǫ (x, t) + ζ 2 B x + 2t ǫ , t + b ǫ x + 2t ǫ , t -2f ǫ (x, t) - ρ π ∂ x δ ǫ x + 2t ǫ , t + ζ 2 ∂ x B x + 2t ǫ , t + b ǫ x + 2t ǫ , t -2f ǫ (x, t) δ ǫ x + 2t ǫ , t + Γ ǫ (x, t) , (3.40)
∂ t δ ǫ (x, t) = -α B(x, t) + b ǫ (x, t) + ζF x - 2t ǫ , t -β -∂ x δ ǫ (x, t) + -α∂ x B(x, t) + b ǫ (x, t) + ζ∂ x F x - 2t ǫ , t δ ǫ (x, t) + ζ 2 f ǫ x - 2t ǫ , t + F (x - 2t ǫ , t -2b ǫ (x, t) + ρ π ∂ x γ ǫ x - 2t ǫ , t + ζ 2 ∂ x f ǫ x - 2t ǫ , t + F (x - 2t ǫ , t -2b ǫ (x, t) γ ǫ x - 2t ǫ , t + ∆ ǫ (x, t) , (3.41)
with α, β + , β -, ζ, Γ ǫ (x, t) and ∆ ǫ (x, t) defined by

α = γ + 1 4 , β + = 2u + (γ -1)ρ 4π , β -= 2u -(γ -1)ρ 4π , ζ = γ -3 2 , (3.42) Γ ǫ (x, t) = -α∂ x W x, t ǫ , t f ǫ (x, t) + F (x, t) -β + ∂ x W x, t ǫ , t -ζ∂ x f ǫ (x, t)V x + 2t ǫ , t ǫ , t -B x + 2t ǫ , t W x, t ǫ , t + ζ 2 ∂ x V x + 2t ǫ , t ǫ , t B x + 2t ǫ , t + b ǫ x + 2t ǫ , t - ζρ 2π ∂ x V x + 2t ǫ , t ǫ , t - 1 2 ∂ x q ǫ 1 ρ ǫ x + t ǫ , t -∂ t W x, t ǫ , t (3.43)
and

∆ ǫ (x, t) = -α∂ x V x, t ǫ , t B(x, t) + b ǫ (x, t) -β -∂ x V x, t ǫ , t -ζ∂ x W x - 2t ǫ , t ǫ , t b ǫ (x, t) -F x - 2t ǫ , t)V x, t ǫ , t + ζ 2 ∂ x W x - 2t ǫ , t ǫ , t f ǫ x - 2t ǫ , t + F x - 2t ǫ , t + ζρ 2π ∂ x W x - 2t ǫ , t ǫ , t + 1 2 ∂ x q ǫ 1 ρ ǫ x - t ǫ , t -∂ t V x, t ǫ , t .
(3.44)

Since the sequences (f ǫ ) ǫ ∈ ]0,ǫ 0 ] and (b ǫ ) ǫ ∈ ]0,ǫ 0 ] are bounded in C 0, T ; H s (T 1 ) with s ≥ 2, F and B are in C 0, T ; H s (T 1 ) , and q ǫ 1 is regular, we can find some constants C i > 0 (i = 1, . . . , 8), independent of ǫ and t, such that 2

T 1 -α f ǫ (x, t) + F (x, t) -β + + ζB x + 2t ǫ , t γ ǫ (x, t) ∂ x γ ǫ (x, t) dx ≤ C 1 γ ǫ (•, t) 2 L 2 (T 1 ) , (3.45) 2 T 1 -α b ǫ (x, t) + B(x, t) -β -+ ζF x - 2t ǫ , t δ ǫ (x, t) ∂ x δ ǫ (x, t) dx ≤ C 2 δ ǫ (•, t) 2 L 2 (T 1 ) , (3.46)
2

T 1 -α∂ x f ǫ (x, t) + F (x, t) + ζ∂ x B x + 2t ǫ , t γ ǫ (x, t) 2 dx ≤ C 3 γ ǫ (•, t) 2 L 2 (T 1 ) , (3.47) 2 T 1 -α∂ x b ǫ (x, t) + B(x, t) + ζ∂ x F x - 2t ǫ , t δ ǫ (x, t) 2 dx ≤ C 4 δ ǫ (•, t) 2 L 2 (T 1 ) , (3.48)
2

T 1 ζ 2 ∂ x B x + 2t ǫ , t + b ǫ x + 2t ǫ , t -2f ǫ (x, t) δ ǫ x + 2t ǫ , t γ ǫ (x, t) dx ≤ C 5 γ ǫ (•, t) L 2 (T 1 ) δ ǫ (•, t) L 2 (T 1 ) , (3.49) 2 T 1 ζ 2 ∂ x F x - 2t ǫ , t + f ǫ x - 2t ǫ , t -2b ǫ (x, t) γ ǫ x - 2t ǫ , t δ ǫ (x, t) dx ≤ C 6 γ ǫ (•, t) L 2 (T 1 ) δ ǫ (•, t) L 2 (T 1 ) , (3.50) 2 T 1 Γ ǫ (x, t) γ ǫ (x, t) dx ≤ C 7 γ ǫ (•, t) L 2 (T 1 ) , (3.51) 2 T 1 ∆ ǫ (x, t) δ ǫ (x, t) dx ≤ C 8 δ ǫ (•, t) L 2 (T 1 ) . (3.52)
On the other hand, making in the first integral the change of variable x → x + τ and in the second one x → xτ , rearranging the terms and using the definitions of γ ǫ and δ ǫ , we have

ζ T 1 b ǫ (x + 2τ, t) + B(x + 2τ, t) -2f ǫ (x, t) - ρ π ∂ x δ ǫ (x + 2τ, t) γ ǫ (x, t) dx + ζ T 1 f ǫ (x -2τ, t) + F (x -2τ, t) -2b ǫ (x, t) + ρ π ∂ x γ ǫ (x -2τ, t) δ ǫ (x, t) dx = ζ T 1 b ǫ (x + τ, t) -B(x + τ, t) γ ǫ (x -τ, t) ∂ x δ ǫ (x + τ, t) dx + 2ζ T 1 F (x -τ, t) -f ǫ (x -τ, t) γ ǫ (x -τ, t) ∂ x δ ǫ (x + τ, t) dx -2ζ T 1 F (x -τ, t) -B(x + τ, t) + ρ 2π γ ǫ (x -τ, t) ∂ x δ ǫ (x + τ, t) dx + ζ T 1 f ǫ (x -τ, t) -F (x -τ, t) δ ǫ (x + τ, t) ∂ x γ ǫ (x -τ, t) dx + 2ζ T 1 B(x + τ, t) -b ǫ (x + τ, t) δ ǫ (x + τ, t) ∂ x γ ǫ (x -τ, t) dx + 2ζ T 1 F (x -τ, t) -B(x + τ, t) + ρ 2π δ ǫ (x + τ, t) ∂ x γ ǫ (x -τ, t) dx = ζ T 1 ǫ ∂ x γ ǫ (x -τ, t) -∂ x δ ǫ (x + τ, t) γ ǫ (x -τ, t) δ ǫ (x + τ, t) dx + ζ T 1 ǫ∂ x γ ǫ (x -τ, t) W (x -τ, τ, t) + 2V (x + τ, τ, t) δ ǫ (x + τ, t) dx -ζ T 1 ǫ∂ x δ ǫ (x + τ, t) V (x + τ, τ, t) + 2W (x -τ, τ, t) γ ǫ (x -τ, t) dx -2ζ T 1 ǫ∂ x δ ǫ (x + τ, t) γ ǫ (x -τ, t) 2 dx + 2ζ T 1 ǫ∂ x γ ǫ (x -τ, t) δ ǫ (x + τ, t) 2 dx + 2ζ T 1 R(x, τ, t) δ ǫ (x + τ, t)∂ x γ ǫ (x -τ, t) -γ ǫ (x -τ, t)∂ x δ ǫ (x + τ, t) dx .
The sequences (ǫγ ǫ ) ǫ ∈ ]0,ǫ 0 ] and (ǫδ ǫ ) ǫ ∈ ]0,ǫ 0 ] are bounded in C 0, T ; H s (T 1 ) , so we can find some constants C i > 0 (i = 9, . . . , 13), independent of ǫ and t, such that

ζ T 1 ǫ ∂ x γ ǫ (x -τ, t) -∂ x δ ǫ (x + τ, t) γ ǫ (x -τ, t) δ ǫ (x + τ, t) dx ≤ C 9 γ ǫ (•, t) L 2 (T 1 ) γ ǫ (•, t) L 2 (T 1 ) , (3.53) ζ T 1 ǫ∂ x γ ǫ (x -τ, t) W (x -τ, τ, t) + 2V (x + τ, τ, t) δ ǫ (x + τ, t) dx ≤ C 10 δ ǫ (•, t) L 2 (T 1 ) , (3.54) ζ T 1 ǫ∂ x δ ǫ (x + τ, t) V (x + τ, τ, t) + 2W (x -τ, τ, t) γ ǫ (x -τ, t) dx ≤ C 11 γ ǫ (•, t) L 2 (T 1 ) , (3.55)
2ζ

T 1 ǫ∂ x δ ǫ (x + τ, t) γ ǫ (x -τ, t) 2 dx ≤ C 12 γ ǫ (•, t) 2 L 2 (T 1 ) , (3.56) 2ζ T 1 ǫ∂ x γ ǫ (x -τ, t) δ ǫ (x + τ, t) 2 dx ≤ C 13 δ ǫ (•, t) 2 L 2 (T 1 ) . (3.57)
For the last integral, we simply compute:

2ζ

T 1 R(x, τ, t) δ ǫ (x + τ, t)∂ x γ ǫ (x -τ, t) -γ ǫ (x -τ, t)∂ x δ ǫ (x + τ, t) dx = 2ζ T 1 ∂ τ R(x, τ, t) δ ǫ (x + τ, t) γ ǫ (x -τ, t) dx -2ζ∂ τ T 1 R(x, τ, t) δ ǫ (x + τ, t) γ ǫ (x -τ, t) dx = 2ζ T 1 ∂ τ R(x, τ, t) δ ǫ (x + τ, t) γ ǫ (x -τ, t) dx = -2ζ T 1 δ ǫ (x + τ, t) γ ǫ (x -τ, t) ∂ x U (x, τ, t) dx ≤ C 14 γ ǫ (•, t) L 2 (T 1 ) δ ǫ (•, t) L 2 (T 1 ) (3
∂ t γ ǫ (•, t) 2 L 2 (T 1 ) + δ ǫ (•, t) 2 L 2 (T 1 ) ≤ M γ ǫ (•, t) 2 L 2 (T 1 ) + δ ǫ (•, t) 2 L 2 (T 1 ) + γ ǫ (•, t) L 2 (T 1 ) + δ ǫ (•, t) L 2 (T 1 ) + γ ǫ (•, t) L 2 (T 1 ) δ ǫ (•, t) L 2 (T 1 ) , (3.59) with M = max(C 1 + C 3 + C 12 , C 2 + C 4 + C 13 , C 5 + C 6 + C 9 + C 14 , C 7 + C 11 , C 8 + C 10 ).
Inequality (3.37) is directly obtained from (3.59), ending the proof of lemma 5 and then the proof of theorem 4.

Properties of the finite volume scheme: proof of the Theorem 2

Thanks to (2.3) and (2.9), convergence of

U h (•, τ, •), R h (•, τ, t) to U (•, τ, •), R(•, τ, t) for any τ ∈ T 1 is equivalent to the convergence of (F h , B h ) to (F, B).
Furthermore, the equations satisfied by F and B are of the form

∂ t q + ∂ x f (q) = 0 , q(x, 0) = q 0 (x) , T 1 q(x, t) dx = 0 , ∀ t , (3.60)
with f (q) = αq 2 + βq, α and β = β ± defined in (3.42). In the same way, the numerical method (2.6)-(2.7)-(2.8) is of the form

Q h (x, t) = Q n i ∀ [x i-1/2 , x i+1/2 [×[t n , t n+1 [ , Q n+1 i = Q n i - k h F(Q n i+1 , Q n i ) -F(Q n i , Q n i-1 ) , F(Q i , Q i-1 ) = 1 2 f (Q i ) + f (Q i-1 ) - 1 2 α(Q i + Q i-1 ) + β (Q i -Q i-1 ) , (3.61)
with a space step h = 2π N x + 1 and a time step k.

We notice that the CFL condition at the n-th time step is given by

k h max 0 ≤ i ≤ Nx α(Q n i + Q n i-1 ) + β = ν ≤ 1 . (3.62)
With these notations, proving Theorem 2 is equivalent to proving the theorem below.

Theorem 6. If s ≥ 2, the approximation Q h converges to the solution q of (3.60) in L 1 [0, T ) × T 1 norm. Furthermore, if s ≥ 3, the local truncation error of the numerical scheme (3.61) is first order accurate.

Proof of theorem 5: the total variation of a function q ∈ L 1 [0, T ) × T 1 is defined by

T V T (q) = lim sup η → 0 1 η T 0 T 1 q(x + η, t) -q(x, t) + q(x, t + η) -q(x, t) dx dt . (3.63)
For an approximation Q h computed with the numerical method (3.61), we have

T V T (Q h ) = T /k n = 0 k T V (Q n ) + Q n+1 -Q n 1 , (3.64)
with T V (Q) and Q 1 defined by

T V (Q) = Nx i = 0 |Q i -Q i-1 | and Q 1 = h Nx i = 0 |Q i | . (3.65)
Introducing the sets L and K defined by

L = q ∈ L 1 [0, T ) × T 1 : T 1 q dx = 0 ∀ t , K = q ∈ L : T V T (q) ≤ R , (3.66)
where R is a constant depending on q 0 , it is well known that K is a compact subset of L.

Since the numerical flux involved in (3.61) is continuous and satisfies F(q, q) = f (q) ∀ q , (3.67) the considered scheme is consistent with the conservation law (3.60). Hence, applying LeVeque [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], proving that Q h converges to the solution q of (3.60) in L 1 [0, T ) × T 1 norm for h → 0 reduces to prove that the scheme is TV-stable. In other words, we need to prove that for h ∈ [0, h 0 ], the approximation Q h lies in some fixed set K where R only depends on the initial data q 0 , the final time T and the function f .

The TV-stability is the consequence of two lemmas. Lemma 7. The numerical method (3.61) is TVD (Total Variation Diminishing), i.e.

T V (Q n+1 ) ≤ T V (Q n ) ∀ n . (3.68)
Proof of the lemma 7:

we develop Q n+1 i+1 -Q n+1 i : Q n+1 i+1 -Q n+1 i = Q n i+1 -Q n i - k 2h f (Q n i+2 ) -f (Q n i ) -f (Q n i+1 ) + f (Q n i-1 ) -α(Q n i+2 + Q n i+1 ) + β (Q n i+2 -Q n i+1 ) + 2 α(Q n i+1 + Q n i ) + β (Q n i+1 -Q n i ) -α(Q n i + Q n i-1 ) + β (Q n i -Q n i-1 ) = - k 2h 1 -sg α(Q n i+2 + Q n i+1 ) + β α(Q n i+2 + Q n i+1 ) + β (Q n i+2 -Q n i+1 ) + 1 - k h α(Q n i+1 + Q n i ) + β (Q n i+1 -Q n i ) + k 2h 1 + sg α(Q n i + Q n i-1 ) + β α(Q n i + Q n i-1 ) + β (Q n i -Q n i-1 ) ,
where sg stands for the usual sign function. Using the CFL condition (3.62), we write:

Q n+1 i+1 -Q n+1 i ≤ k 2h 1 -sg α(Q n i+2 + Q n i+1 ) + β α(Q n i+2 + Q n i+1 ) + β |Q n i+2 -Q n i+1 | + 1 - k h α(Q n i+1 + Q n i ) + β |Q n i+1 -Q n i | + k 2h 1 + sg α(Q n i + Q n i-1 ) + β α(Q n i + Q n i-1 ) + β |Q n i -Q n i-1 | .
We deduce then

T V (Q n+1 ) = Nx i = 0 Q n+1 i+1 -Q n+1 i ≤ Nx i = 0 Q n i+1 -Q n i = T V (Q n ) (3.69)
giving the lemma.

As a consequence of this lemma, we deduce that there exists a constant M 1 > 0 which only depends on the initial data q 0 and such that

T V (Q n ) ≤ M 1 ∀ n .
(3.70) Lemma 8. For any n, we have

Q n+1 -Q n 1 ≤ 2M 1 h.
Proof of lemma 8: we develop

Q n+1 i -Q n i : Q n+1 i -Q n i = - k 2h f (Q n i+1 ) -f (Q n i ) -α(Q n i+1 + Q n i ) + β (Q n i+1 -Q n i ) -f (Q n i ) + f (Q n i-1 ) + α(Q n i + Q n i-1 ) + β (Q n i -Q n i-1 ) = - k 2h 1 -sg α(Q n i+1 + Q n i ) + β α(Q n i+1 + Q n i ) + β (Q n i+1 -Q n i ) + k 2h 1 -sg α(Q n i + Q n i-1 ) + β α(Q n i + Q n i-1 ) + β (Q n i -Q n i-1 ) .
Then, using the CFL condition, we obtain

Q n+1 -Q n 1 = h Nx i = 0 Q n+1 i -Q n i ≤ h Nx i = 0 k 2h α(Q n i+1 + Q n i ) + β |Q n i+1 -Q n i | + h Nx i = 0 k 2h α(Q n i + Q n i-1 ) + β |Q n i -Q n i-1 | ≤ 2h Nx i = 0 |Q n i -Q n i-1 |
and we complete the proof by applying lemma 7.

End of the proof of theorem 5: combining all these results, we finally have the inequality

T V T (Q h ) ≤ 1 + 2h k T M 1 . (3.71)
Remarking that h k is bounded by a constant M 2 which only depends on the initial data q 0 , we finally obtain that there exists a constant C > 0 which only depends on q 0 and T , and such that

T V T (Q h ) ≤ C (3.72)
for any space step h. We conclude that the numerical method (3.61) is TV-stable, yielding the convergence of Q h to q. Now, we have to prove that, assuming that s ≥ 3, the local truncation error is first order accurate. With the notation above, the expression of the local truncation error reads

e(x i , t n ) = q n i -q n+1 i k - 1 2h f (q n i+1 ) -f (q n i-1 ) -α(q n i+1 + q n i ) + β (q n i+1 -q n i ) + α(q n i + q n i-1 ) + β (q n i -q n i-1 ) , (3.73)
where q n i = q(x i , t n ). According to the signs of α(q n i + q n i-1 ) + β and α(q n i+1 + q n i ) + β, we distinguish four different cases.

• First case: we assume that α(q n i+1 + q n i ) + β ≥ 0 and α(q n i + q n i-1 ) + β ≥ 0. Then, the local truncation error becomes

e(x i , t n ) = q n i -q n+1 i k - f (q n i ) -f (q n i-1 ) h . (3.74) Since s ≥ 3, q ∈ C 2 [0, T ) × T 1 , so we can write e(x i , t n ) = -∂ t q(x i , t n ) -∂ x f q(x i , t n ) + O(k) + O(h) = O(h + k) . (3.75)
• Second case: we assume that α(q n i+1 + q n i ) + β ≤ 0 and α(q n i + q n i-1 ) + β ≤ 0. Then, the local truncation error becomes

e(x i , t n ) = q n i -q n+1 i k - f (q n i+1 ) -f (q n i ) h . (3.76) Since s ≥ 3, q ∈ C 2 [0, T ) × T 1 , so we can write e(x i , t n ) = -∂ t q(x i , t n ) -∂ x f q(x i , t n ) + O(k) + O(h) = O(h + k) . (3.77)
• Third case: we assume that α(q n i+1 + q n i ) + β ≥ 0 and α(q n i + q n i-1 ) + β ≤ 0. Then, the function λ defined by

λ(x) = α q(x, t n ) + q(x -h, t n ) + β (3.78) admits a zero denoted x * in [x i , x i+1 ]. Since s ≥ 3, λ is of class C 2 and verifies λ(x i ) = λ(x * ) + O(x i -x * ) = O(h) , λ(x i+1 ) = λ(x * ) + O(x i+1 -x * ) = O(h) . (3.79)
In the same way, we have

f ′ (q n i ) = 2αq n i + β = λ(x i+1 ) + O(h) = O(h) . (3.80)
The local truncation error is of the form

e(x i , t n ) = q n i -q n+1 i k = -∂ t q(x i , t n ) + O(k) = f ′ q(x i , t n ) ∂ x q(x i , t n ) + O(k) .
(3.81) Then, using (3.80), we obtain e(x i , t n ) = O(h + k) .

(3.82)

• Fourth case: we assume that α(q n i+1 + q n i ) + β ≤ 0 and α(q n i + q n i-1 ) + β ≥ 0. We proceed as we did in the precedent case, i.e. we deduce that the function λ defined by (3.78) admits a zero x * ∈ [x i , x i+1 ], and then the result (3.80) is true. Hence,

e(x i , t n ) = q n i -q n+1 i k - f (q n i+1 ) -f (q n i-1 ) h = -∂ t q(x i , t n ) -2∂ x f q(x i , t n ) + O(h + k) = -f ′ (q n i )∂ x q(x i , t n ) + O(h + k) = O(h + k) . (3.83)
We conclude that, in any case, the local truncation error is O(h + k), i.e. first order accurate.

Numerical results

The first goal of this section is to numerically show that u

ǫ (x, t) -U h (x, t ǫ , t) = O(ǫ) and ρ ǫ (x, t)-R h (x, t ǫ , t) = O(ǫ).
Secondly, we will briefly analyze the gain in term of CPU time our two-scale numerical method brings when compared with a classical method consisting in solving directly (1.1). Lastly, by using our method, we will explore the simulation of experiments which were not accessible before because of too small Mach number.

Convergence in ǫ

We consider here the initial data u 0 and ρ 0 defined by

u 0 (x) = 1 + cos(x) 2 and ρ 0 (x) = 1 + sin(x) 2 , (4.1)
and we take a uniform mesh with N x + 1 points denoted x 0 , . . . , x Nx . As illustrated in figures 1 and 2, the functions U h (x, t ǫ , t) and R h (x, t ǫ , t) are very close to u ǫ and ρ ǫ respectively. The numerical experiment showed in these figures is made with ǫ = 0.05, γ = 1, and N x = 1023, T = 2.5. In order to quantify this good accuracy of U h (x, t ǫ , t) with u ǫ and R h (x, t ǫ , t) with ρ ǫ , we compute the errors u ǫ (x, t) -U h (x, t ǫ , t) and ρ ǫ (x, t) -R h (x, t ǫ , t) in L p [0, T ) × T 1 norm (p = 1, 2, ∞), for several values of ǫ ranging in {0.01, 0.03, 0.05, 0.07, 0.1}, and we obtain the figures and the array below:

3: Error u ǫ (x, t) -U h (x, t ǫ , t) in L 1 , L
K ∞ ǫ with K ∞ ≈ 1.3792216 in L ∞ [0, T ) × T 1 norm. Concerning the error ρ ǫ (x, t) -R h (x, t ǫ , t)
, we observe that it decreases when ǫ → 0 as

K ′ 1 ǫ with K ′ 1 ≈ 3.5842834 in L 1 [0, T ) × T 1 norm, as K ′ 2 ǫ with K ′ 2 ≈ 1.1780144 in L 2 [0, T ) × T 1 norm,

and and as

K ′ ∞ ǫ with K ′ ∞ ≈ 1.690496 in L ∞ [0, T ) × T 1 norm.

CPU time cost

One of the motivations of the present paper is the high CPU time cost of classical methods when used with small Mach number, as we explained it in the introduction. For example, if we apply Roe's finite volume scheme on the model (1.1), we must verify the CFL condition

k h max i , n ûn i ± √ P n i ǫ ≤ 1 , (4.2)
where h is the space step, k is the time step, ûn i and P n i are defined by

ûn i = u n i 1 + ǫρ n i + u n i-1 1 + ǫρ n i-1 1 + ǫρ n i + 1 + ǫρ n i-1 , P n i = (1+ǫρ n i ) γ -(1+ǫρ n i-1 ) γ γ ǫ(ρ n i -ρ n i-1 ) if ρ n i = ρ n i-1 , (1 + ǫρ n i ) γ-1 else . (4.3)
Hence, in order to garantee the stability of the scheme, the smaller ǫ is, the smaller k must be (and the higher the CPU time cost will be). One advantage of our two-scale numerical method is that its CFL condition (3.62) does not depend on ǫ, so k does not have to diminish with ǫ.

Value

Roe's method Two-scale numerical method on (2. The array above gives the comparison in terms of CPU time cost between Roe's method applied on the non-homogenized model (1.1), and our two-scale numerical method. These results have been obtained with the computer characteristics below:

• processor: Pentium c M 715,

• memory: 512 Mo DDR-RAM,

• operating system: SuSe c 9.1 Pro,

• compiler: gcc 3.3.3-41.

In particular, we distinguish the time cost for the computation of F h and B h on one hand, and the time cost for the reconstruction of U h (x, t ǫ , t) and R h (x, t ǫ , t) on the other hand.

Numerical experiments with small Mach number

In this paragraph, we present some numerical results obtained with our two-scale numerical method on simulations of experiments inducing a very small Mach number. These results are obtained with the initial data (4.1), N x + 1 = 1024 points in x, and γ = 1. We can notice that these results have been obtained after a few minutes of computation just as described in the table below:

Value of ǫ Computation of (F h , B h ) Computation of (U h , R h ) 10 -4 3 m 45 s 71" 42 s 84" 10 -5 3 m 45 s 95" 42 s 87" 10 -6 3 m 45 s 75" 42 s 98" Furthermore, if we continue the simulation to final time T ≈ 3.2 for ǫ ranging in {10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 , 10 -6 }, we remark that discontinuities appear at time t ≈ Table 3: CPU time costs.

3.11 independently of ǫ and that, for ǫ ranging in {0.01, 0.03, 0.05, 0.07, 0.1}, these discontinuities also appear in the solution (u ǫ , ρ ǫ ) of (1.1) at the same time.

Conclusion

In this paper, we have developped a two-scale numerical method for the weakly compressible 1D Euler equations. By using two-scale convergence tools from Nguetseng [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] and Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], we have developped a new model independent of the Mach number ǫ. Then, using a finite volume scheme on this model, we have obtained numerical approximations that converge to the solutions of the weakly compressible 1D Euler equation when we refine the mesh on T 1 × [0, T ] and when ǫ → 0.

Furthermore, we have proved that this two-scale numerical method is first order accurate in ǫ. One the other hand, numerical results have confirmed this accuracy in ǫ and also have proved that, even with a very small Mach number, the two-scale numerical method allows us to simulate experiments with a reasonable CPU time cost.

Motived by the behaviour of the two-scale method, we are now ready to develop such a method on a model which describes a plasma submitted to a strong magnetic field such as the Vlasov-Poisson model established in Frénod and Sonnendrücker [START_REF] Frénod | The finite Larmor radius approximation[END_REF]. Furthermore, we only studied the behaviour of the two-scale numerical method in front of smooth initial data: it can be interesting to study the behaviour of the method on non-smooth initial data. Concerning the periodicity, we have built our numerical method in a periodic framework for mathematical reasons. However, it must be possible to extend this method to almost periodic framework in order to study problems like the noise generated by the blades of a turbine or a fan. If we go further, it can be interesting to improve the method in order to use it for studying the noise generated by multidimensional structures like cooling systems for computers or air conditioning, even if the links between time scales and space scales are not clear in this framework.

Figure 2 :

 2 Figure 2: Comparison of ρ ǫ (•, t) and R h (•, t ǫ , t) at times t = 0 and t = 2.5.

  with K = 2.8893757 K*epsilon with K = 1.0358362 K*epsilon with K = 1.3792216

Figure

  Figure 3: Error u ǫ (x, t) -U h (x, t ǫ , t) in L 1 , L 2 and L ∞ norms.

Table 1 :

 1 2 and L ∞ norms. Figure 4: Error ρǫ (x, t) -R h (x, t ǫ , t) in L 1 , L 2 and L ∞ norms. Errors u ǫ (x, t) -U h (x, t ǫ , t) and ρ ǫ (x, t) -R h (x, t ǫ , t).The results above show an error u ǫ (x, t) -U h (x, t ǫ , t) in L 1 norm decreasing when ǫ → 0 as K 1 ǫ with K 1 ≈ 2.8893757, as K 2 ǫ with K 2 ≈ 1.0358362 in L 2 [0, T ) × T 1 norm, and and as

	0.40
	0.36					Ο
	0.32
	0.28
					Ο
	0.24
	0.20
				Ο		×
	0.16
	0.12				×	+
	0.08		Ο	×	+
			×	+
	0.04	Ο ×	+
		+
	0
	0	Ο + ×	0.02 error in L^1 norm error in L^2 norm error in Max norm	0.04	0.06	0.08 K*epsilon with K = 3.5842834 0.10 K*epsilon with K = 1.1780144 K*epsilon with K = 1.690496	0.12

Table 2 :

 2 Comparison in terms of CPU time cost between Roe's method on (1.1) and the two-scale numerical method on (2.3)-(2.5) with the initial data (4.1), N x = 1023, γ = 1.

	3)-(2.5)

  Figure8: U h (x, t ǫ , t) (left) and R h (x, t ǫ , t) (right) with ǫ = 10 -6 at times t = 0, t = 1, t = 2.

			Function x -> U(x,t/epsilon,t) with t = 2 and epsilon = 0.0001 Function x -> U(x,t/epsilon,t) with t = 0 and epsilon = 0.000001					Function x -> R(x,t/epsilon,t) with t = 2 and epsilon = 0.0001 Function x -> R(x,t/epsilon,t) with t = 0 and epsilon = 0.000001
	3.0	3.0															3.0	3.0
	2.7	2.7															2.7	2.7
	2.4	2.4															2.4	2.4
	2.1	2.1															2.1	2.1
	1.8	1.8															1.8	1.8
	1.5	1.5															1.5	1.5
	1.2	1.2															1.2	1.2
	0.9	0.9															0.9	0.9
	0.6	0.6															0.6	0.6
	0.3	0.3															0.3	0.3
	0.0	0.0															0.0	0.0
	0.000	0.000 0.628	0.628 1.255	1.255 1.883	1.883 2.511	2.511 3.139	3.139 3.766	3.766 4.394	4.394 5.022	5.022 5.649	5.649 6.277	6.277 0.000	0.000 0.628	0.628 1.255	1.255 1.883	1.883 2.511	2.511 3.139	3.139 3.766	3.766 4.394	4.394 5.022	5.022 5.649	5.649 6.277	6.277
		U(x,t/epsilon,t) U(x,t/epsilon,t)												R(x,t/epsilon,t) R(x,t/epsilon,t)
	Figure 6: U h (x, t ǫ , t) (left) and R h (x, t ǫ , t) (right) with ǫ = 10 -4 at time t = 2.
			Function x -> U(x,t/epsilon,t) with t = 1 and epsilon = 0.000001							Function x -> R(x,t/epsilon,t) with t = 1 and epsilon = 0.000001
	3.0																3.0
	2.7				Function x -> U(x,t/epsilon,t) with t = 0 and epsilon = 0.00001			2.7				Function x -> R(x,t/epsilon,t) with t = 0 and epsilon = 0.00001
	2.4	3.0															2.4	3.0
	2.1	2.7															2.1	2.7
	1.8	2.4															1.8	2.4
	1.5	2.1															1.5	2.1
	1.2	1.8															1.2	1.8
	0.9	1.5															0.9	1.5
	0.6	1.2															0.6	1.2
	0.3	0.9															0.3	0.9
	0.000 0.0	0.6 0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277	0.000 0.0	0.628 0.6	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
		0.3 U(x,t/epsilon,t)														0.3 R(x,t/epsilon,t)
		0.0																0.0
		0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
				U(x,t/epsilon,t)														R(x,t/epsilon,t)
			Function x -> U(x,t/epsilon,t) with t = 2 and epsilon = 0.000001							Function x -> R(x,t/epsilon,t) with t = 2 and epsilon = 0.000001
	3.0				Function x -> U(x,t/epsilon,t) with t = 1 and epsilon = 0.00001			3.0				Function x -> R(x,t/epsilon,t) with t = 1 and epsilon = 0.00001
	2.7	3.0	Function x -> U(x,t/epsilon,t) with t = 0 and epsilon = 0.0001					2.7	3.0		Function x -> R(x,t/epsilon,t) with t = 0 and epsilon = 0.0001
	3.0 2.4	2.7															3.0 2.4	2.7
	2.7 2.1	2.4															2.7 2.1	2.4
	2.4 1.8	2.1															2.4 1.8	2.1
	2.1 1.5	1.8															2.1 1.5	1.8
	1.8 1.2	1.5															1.8 1.2	1.5
	1.5 0.9	1.2															1.5 0.9	1.2
	1.2 0.6	0.9															1.2 0.6	0.9
	0.9 0.3	0.6															0.9 0.3	0.6
	0.6 0.0	0.3															0.6 0.0	0.3
	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
	0.3	0.000 0.0 U(x,t/epsilon,t) 0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277 0.3	0.000 0.0 R(x,t/epsilon,t) 0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
				U(x,t/epsilon,t)														R(x,t/epsilon,t)
	0.0																0.0
	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
		U(x,t/epsilon,t)														R(x,t/epsilon,t)
					Function x -> U(x,t/epsilon,t) with t = 2 and epsilon = 0.00001							Function x -> R(x,t/epsilon,t) with t = 2 and epsilon = 0.00001
		3.0																3.0
			Function x -> U(x,t/epsilon,t) with t = 1 and epsilon = 0.0001								Function x -> R(x,t/epsilon,t) with t = 1 and epsilon = 0.0001
	3.0	2.7															3.0	2.7
	2.7	2.4															2.7	2.4
	2.4	2.1															2.4	2.1
	2.1	1.8															2.1	1.8
	1.8	1.5															1.8	1.5
	1.5	1.2															1.5	1.2
	1.2	0.9															1.2	0.9
	0.9	0.6															0.9	0.6
	0.6	0.3															0.6	0.3
	0.3	0.000 0.0	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277 0.3	0.000 0.0	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
				U(x,t/epsilon,t)														R(x,t/epsilon,t)
	0.0																0.0
	0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277		0.000	0.628	1.255	1.883	2.511	3.139	3.766	4.394	5.022	5.649	6.277
		U(x,t/epsilon,t)														R(x,t/epsilon,t)

Figure

5

: U h (x, t ǫ , t) (left) and R h (x, t ǫ , t) (right) with ǫ = 10 -4 at times t = 0, t = 1.

Figure

7

: U h (x, t ǫ , t) (left) and R h (x, t ǫ , t) (right) with ǫ = 10 -5 at times t = 0, t = 1, t = 2.


Georgios Bouloukakis

Ioannis Moscholios

Nikolaos Georgantas

Valérie Issarny

Simulation-based Queueing Models for Performance Analysis of IoT Applications

Keywords: Queueing Models, Mobile Connectivity, Internet of Things, QoS Analysis

To facilitate the development of Internet of Things (IoT) applications, numerous middleware protocols and APIs have been introduced. Such applications built atop reliable or unreliable protocols and they expose different characteristics. Additionally, with regard to the application context (e.g., emergency response operations), several Quality of Service (QoS) requirements must be satisfied. To study QoS in IoT applications, the provision of a generic performance analysis methodology is required. Queueing network models offer a simple modeling environment, which can be used to represent IoT interactions by combining multiple queueing model types for building queueing networks. The resulting networks can be used for performance analysis through analytical or simulation models. In this paper, we present several types of queueing models that represent different QoS settings of IoT interactions, such as intermittent mobile connectivity, message drop probabilities, message availability/validity and resource constrained devices. Using MobileJINQS, we simulate our models demonstrating the significant effect on response times and message success rates when varying QoS settings.

I. INTRODUCTION

The Internet of Things (IoT) promises the integration of the physical world into computer-based systems. IoT devices, featuring sensing capabilities, are deployed in a variety of application domains, such as smart buildings, community spaces, intelligent transportation to name a few. By exploiting their information, there are new opportunities to improve peoples' safety and quality of life. For instance in [START_REF]Community Seismic Network[END_REF], [START_REF] Cochran | A novel strongmotion seismic network for community participation in earthquake monitoring[END_REF], IoT devices are used by volunteers to monitor homes/offices for possible indications of seismic activity. Such an application generates critical information and is expected to function correctly, timely and reliably. Hence, identifying a general methodology for the performance modeling and analysis of IoT applications is of prime importance.

IoT devices can be mobile, low-powered, inexpensive and with regard to the environmental context, the resulting applications expose different characteristics. In particular, several

The work is supported by the research associate team ACHOR (inria.fr/en/associate-team/achor).

Quality of Service (QoS) constraints can be identified: i) intermittent connectivity of mobile devices; ii) limited memory of inexpensive devices; iii) limited availability of data in emergency response scenarios, etc. Furthermore, IoT applications built atop reliable/unreliable protocols and Application Programming Interfaces (APIs) which introduce additional QoS constraints [START_REF] Caro | Comparison of two lightweight protocols for smartphone-based sensing[END_REF], [START_REF] Lee | Correlation analysis of mqtt loss and delay according to qos level[END_REF]. Hence, these constraints raise important challenges for the supporting of IoT performance analysis.

Existing efforts [START_REF] Mehmeti | Performance analysis of "on-thespot" mobile data offloading[END_REF]- [START_REF] Wu | Tradeoff analysis for mobile cloud offloading based on an additive energy-performance metric[END_REF] concerning the design and evaluation of mobile systems under specific constraints (e.g., intermittent availability, limited resources, etc) and for specific use cases rely on the field of Queueing Theory [START_REF] Gross | Fundamentals of queueing theory[END_REF]. Key IoT protocols have been evaluated with regard to performance metrics such as response times and delivery success rates [START_REF] Durkop | Performance evaluation of m2m protocols over cellular networks in a lab environment[END_REF], [START_REF] Fysarakis | Which iot protocol? comparing standardized approaches over a common m2m application[END_REF]. However, such efforts are protocol-specific and thus, the research community has analyzed the performance of well known interactions paradigms [START_REF] Aldred | On the notion of coupling in communication middleware[END_REF], [START_REF] Kattepur | Analysis of timing constraints in heterogeneous middleware interactions[END_REF] that may abstract different IoT protocols. The publish/subscribe paradigm has been analyzed using formal models [START_REF] He | Formal analysis of publish-subscribe systems by probabilistic timed automata[END_REF] and evaluated using models such as Queueing Petri Nets [START_REF] Kounev | A methodology for performance modeling of distributed event-based systems[END_REF] or Queueing networks [START_REF] Lazowska | Quantitative system performance: computer system analysis using queueing network models[END_REF].

Queueing networks have been extensively applied to represent and analyze communication and computer systems. A queueing network is a network of connected service centers which provides analytical or simulation solutions for performance measures (e.g., response time). Combining separate service centers in order to form a queueing network, enable us to model and represent an IoT interaction. Towards this, in [START_REF] Bouloukakis | Timeliness evaluation of intermittent mobile connectivity over pub/sub systems[END_REF]- [START_REF] Bouloukakis | Enabling emergent mobile systems in the iot: from middleware-layer communication interoperability to associated qos analysis[END_REF] we have identified a methodology where middleware protocol nodes (clients, servers, brokers, etc) are represented as queues and the exchanged messages as jobs served. To model QoS settings such as the intermittent connectivity [START_REF] Bajaj | Toward enabling convenient urban transit through mobile crowdsensing[END_REF] of mobile Things, we have investigated a separated queueing center via the so called ON/OFF queueing model in [START_REF] Bouloukakis | Timeliness evaluation of intermittent mobile connectivity over pub/sub systems[END_REF], [START_REF] Bouloukakis | Performance modeling of the middleware overlay infrastructure of mobile things[END_REF]. Then, we demonstrated the applicability of our approach by modeling the performance of publish/subscribe [START_REF] Bouloukakis | Timeliness evaluation of intermittent mobile connectivity over pub/sub systems[END_REF], [START_REF] Bouloukakis | Leveraging cdr datasets for context-rich performance modeling of largescale mobile pub/sub systems[END_REF] and data streaming protocols [START_REF] Bouloukakis | Performance modeling of the middleware overlay infrastructure of mobile things[END_REF]. In [START_REF] Bouloukakis | Enabling emergent mobile systems in the iot: from middleware-layer communication interoperability to associated qos analysis[END_REF] we applied our approach to analyze the performance of IoT devices with heterogeneous QoS settings. In this paper, we provide additional queueing models that represent the aforementioned QoS constraints of IoT applications. The key contributions of this work are: Fig. 1: M/M/1 and ON/OFF queues.

1) The description of a generic ON/OFF model extended with message join probabilities and message losses.

2) The introduction of additional features for queueing models representing resource constrained devices and message availability.

3) The comparison of different queueing models in terms of response time and message success rate through simulation-based experiments. The rest of this paper is organized as follows: Section II, provides a description of our queueing models that represent a variety of QoS settings found in IoT applications. In Section III, we compare the performance of queueing models when applying message join probabilities, message availabilities and finite buffer capacities. We finally complement this paper with a brief conclusion and future work in Section IV.

II. QUEUEING MODELS

In this section, we define the individual queueing models that are used as part of the simulation-based queueing networks of our methodology.

A. M/M/1 model

This queue models uninterrupted serving (transmission, reception or processing) of messages as part of an end-to-end IoT interaction. It corresponds to the most common M/M/1 queue (see Fig. 1a), featuring Poisson arrivals and exponential service times.

An M/M/1 queue (q m/m/1 ) is defined by the tuple:

q m/m/1 = (λ, µ). ( 1
)
where λ is the input rate of messages to the queue and µ is the service rate for the processing of messages. Let D be the service demand for the processing delay of each message (i.e., D = 1/µ). Based on standard solutions for the M/M/1 queue [START_REF] Lazowska | Quantitative system performance: computer system analysis using queueing network models[END_REF], the time that a message remains in the system (corresponding to queueing time + service time; we also call it mean response time) is given by:

∆ m/m/1 = D 1 -λD . (2)

B. ON/OFF model

To deal with the mobile peer's connections and disconnections we introduce the Intermittent (ON/OFF) queue, which is depicted in Fig. 1b. Messages arrive according to a Poisson process with rate λ > 0, and are placed in a queue waiting to be "served" (waiting area in Fig. 1b). Messages are served with rate µ >0, which is exponentially distributed. We assume that the server is subject to an on-off procedure. That said, it remains in the ON-state for an exponentially distributed time with parameter θ ON , during which it serves messages (if any). Let T ON be the time the server is ON -i.e., T ON = 1/θ ON . Upon the expiration of this time the server enters the OFF-state during which it stops working (stops serving relevant messages) for an exponentially distributed time period with rate θ OFF . Let T OFF be the time the server is OFF -i.e., T OFF = 1/θ OFF . Accordingly, an ON/OFF queue q on/off is defined by the tuple:

q on/off = (λ, λ on/off , µ, T ON , T OFF ). ( 3
)
where λ is the input rate of messages to the queue, λ on/off is the output rate of messages, and µ is the service rate for the processing of messages (if any) during T ON . The output process λ on/off is intermittent, because no messages exit the queue during T OFF intervals. Without loss of generality, we make the following assumption: if T ON expires and there is a message currently being served, the server interrupts its processing and will continue in the next T ON period. Let ∆ on/of f be the the mean response time (the time that a message remains in the system) for the q on/off queue. In our recent works [START_REF] Bouloukakis | Timeliness evaluation of intermittent mobile connectivity over pub/sub systems[END_REF], [START_REF] Bouloukakis | Performance modeling of the middleware overlay infrastructure of mobile things[END_REF], we have investigated an analytical solution for estimating ∆ on/of f by considering average periods of connections and disconnections (T ON , T OFF ). This is given by:

∆ on/of f = E(n) on/of f λ . (4)
where E(n) on/of f is the average number of messages in the system (server + queue). Its formula is defined in [START_REF] Bouloukakis | Performance modeling of the middleware overlay infrastructure of mobile things[END_REF].

C. Probabilistic ON/OFF model

To reduce queueing delays in the ON/OFF model (see subsection II-B), we introduce the probabilistic ON/OFF model which is depicted in Fig. 2a. In such a model, messages arrive in the system according to a Poisson process with rate λ. If the server is online (ON) then messages are placed in a queue waiting to be served. Otherwise, if the server is offline (OFF), messages decide to join the queue with probability ζ or leave the system with probability 1-ζ. A probabilistic ON/OFF queue q prob on/off is defined by the tuple:

q prob on/off = (λ, λ on/off , µ, ζ, T ON , T OFF ). ( 5
)
where λ is the input rate of messages, ζ is the probability of joining the queue, λ on/off is the output rate of messages and µ is the service rate for the transmission of messages (if any) during T ON . It is worth noting that if ζ = 1, then all messages join the system during OFF periods. This is equivalent with the ON/OFF model presented in the subsection II-B. Hence, the probabilistic ON/OFF model generalizes the ON/OFF model. ∆ on/of f prob is the the mean response time (the time that a message remains in the system) for the q prob on/off queue. We are currently working on identifying an analytical model for estimating ∆ on/of f prob . Using the probabilistic ON/OFF model, system designers will be able to appropriate tune/improve delays of components that may disconnect for an average period T OFF and transmit messages during T ON over a limited bandwidth. Reducing the arrival rate of messages during T OFF , decreases queueing delays and prevents network flooding during T ON . Such a component may be a resource-constrained (in terms of energy and memory) IoT sensor that transits messages over limited bandwidth and may be on sleep-mode periodically (ON/OFF). This probabilistic behavior (i.e., the probability of transmitting messages during T OFF ) must be defined from the corresponding application-layer of the IoT device.

D. ON/OFF loss model

The ON/OFF models defined in the above subsections do not assume message drops (or losses) from the systemany message is either served or buffered and then served. The probabilistic ON/OFF model assumes messages losses, however such messages do not enter the queue since they decide to join the queue with probability ζ.

To model message losses we introduce the ON/OFF loss queue (see Fig. 2b). Similarly to the above ON/OFF queues, the arrival rate (λ) of messages to the system is Poisson. If the server is ON then messages are placed in a queue waiting to be "served" with service rate µ. Otherwise, (i.e., server is OFF) all messages that join the queue (ζ = 1), are "served" with service rate µ loss and afterwards they exit the system (i.e., we have lost messages). We assume that lost messages are "served" with service rate µ loss > µ.

The latter assumption is introduced due to the following reason: we have identified the ON/OFF loss queue in order to model the performance of an IoT mobile device that employs a middleware protocol. Such a protocol introduces message losses since it builds atop UDP and it does not setup any logical session between the sender and the receiver IoT devices. Messages are sent without any guarantee to the other endpoint. As already defined, the service rate µ can be parameterized based on the network transmission delay of the corresponding end-to-end interaction (sender to receiver). On the other hand, to parameterize the µ loss we must identify at which point of the network messages are lost. For instance, assuming that the ON/OFF loss queue models a wireless receiver and its intermittent connectivity, then µ loss can be parameterized through the corresponding network transmission delay. This delay can be defined based on the interaction between the sender and the access point that the receiver connects and disconnects. An ON/OFF loss queue q loss on/off is defined by the tuple:

q loss on/off = (λ, λ out , λ loss , µ, µ loss , T ON , T OFF ). ( 6
)
where λ is the input rate of messages, λ out and µ are the output and processing rates of messages during the ON periods, λ loss and µ loss are the lost and processing rates during the OFF periods. The λ out /λ loss rates are both exponential, since arrived messages split based on the T ON /T OFF intervals. ∆ on/of f loss is the the mean response time for the q loss on/off queue. We are currently working on identifying an analytical model for estimating ∆ on/of f loss .

E. Additional Features

Up to now, we have defined queueing models having buffers with infinite capacity and arriving messages with infinite lifetime (i.e., a time-to-live period which defines the message's availability inside the queueing network). This certainly affects the response time but also the rate of messages successfully served over the total number of arriving messages. However, modeling IoT interactions with such characteristics may not be realistic. For instance, upon a long disconnection period (e.g., 30 mins) of an IoT sensor, the produced data/messages may exceed the sensor's buffer capacity and/or some of the oldest data may become obsolete for the receiving application/user. Accordingly, in this subsection we introduce the corresponding features that take into account the above constraints. These features can be applied to the queueing models defined in the previous subsections.

1) Queueing network with lifetime messages: As already pointed out, a queueing network is a network of connected queues which can be used to model the performance of a system. For instance, in the context of our work, we model IoT interactions by creating queueing networks which consist of queueing models presented in the previous subsections. As depicted in Fig. 3a, messages arrive with a rate λ in order to be processed in the first queue of the queueing network. An arriving message carries a lifetime period attributed to it upon its creation, which represents the message validity inside the queueing network. Hence, a message may enter the queueing network and as soon as its lifetime elapses, the message leaves the network and is considered as expired.

To consider a message as expired, we take into account the time the message spends in both queue and server at each queue. Assuming that a queueing network consists of a single M/M/1 queue, a message reneges if its service does not begin by a certain lifetime period (which includes its expected service demand as well). Based on the queueing theory literature, such a model is studied as an M/M/1 queue with reneging or impatient customers [START_REF] Montazer-Haghighi | On a multiserver markovian queueing system with balking and reneging[END_REF]- [START_REF] Yue | Optimal performance analysis of an M/M/1/N queue system with balking, reneging and server vacation[END_REF]. In this work, we provide the simulation of the above M/M/1 model through our simulator. Furthermore, we enrich the ON/OFF models to support reneging or impatient customers. Hence, system designers are able to create queueing networks that may include different queueing models (described in the previous subsection) enriched with lifetime periods.

2) Queues with finite capacity: This is a well known queueing model feature, where a specific buffer size is applied to the queue that ensures having max K messages in the system (queue + server). This prevents from storing too many messages for too long in devices with limited hardware capacity (memory, hard disk). In particular, as depicted in Fig. 3b, messages arrive in the queue with λ in . Before a message enters the queue the following condition is checked: new queue size + message in service > K. If the condition is true, the message is dropped. Otherwise, the message enters the queue to be processed.

Based on the literature, an M/M/1 queue with finite capacity is notated as M/M/1/K [START_REF] Gross | Fundamentals of queueing theory[END_REF]. In our models, we represent M/M/1 and the ON/OFF queues presented in the previous subsections with finite capacity by adding the system size (K) to the corresponding definition (see tuples [START_REF]Community Seismic Network[END_REF][START_REF] Caro | Comparison of two lightweight protocols for smartphone-based sensing[END_REF][START_REF] Mehmeti | Performance analysis of "on-thespot" mobile data offloading[END_REF][START_REF] Lee | Mobile data offloading: How much can wifi deliver?[END_REF].

III. EXPERIMENTAL RESULTS

In [START_REF] Bouloukakis | Performance modeling of the middleware overlay infrastructure of mobile things[END_REF], we presented our simulator -MobileJINQS, which is an extension of the Java Implementation of a Networkof-Queues Simulation (JINQS) [START_REF] Field | Jinqs: An extensible library for simulating multiclass queueing networks[END_REF]. In this work we further extend MobileJINQS to implement the ON/OFF models described in Section II.

For our experimental setup, each different ON/OFF model is utilized to represent a mobile sender with different QoS settings when sending messages. Such a mobile peer generates approximately 7,500,000 messages to accurately calculate permessage mean response times and success rates. It connects and disconnects in the scale of seconds. To represent such behavior for any experiment, we set the ON/OFF model parameters as follows: i) the server remains in the ON and OFF states for exponentially distributed time periods T ON = T OFF = 20 sec, thus, the server changes its state every 20 sec; ii) messages are processed with a mean service demand D = 0.125 sec and iii) messages arrive to the queue with a mean arrival rate λ, varying from 0.05 to 3.9 messages per sec (λ max = 3.9 messages/sec). Then, we setup 4 different experiments which correspond to different ON/OFF models by applying the following QoS settings: 1) join probability (ζ) = 1, message availability (lifetime) = ∞ and buffer capacity

(K) = ∞; 2) ζ = 1, lifetime = 30 sec, K = ∞; 3) ζ = 0.75, lifetime = ∞, K = ∞; and 4) ζ = 1, lifetime = ∞, K = 100.

A. Response Times

For each one of the above (4) experiments we run its simulator model (based on the corresponding ON/OFF model) 10 times and average across these runs to derive the simulated curves depicted in Fig. 4. As expected, the highest level of mean response times (∼10-360 sec) are noticed in the 1st experiment where all messages must be buffered and served (without any message losses). At any λ rate, the lowest level of mean response times (∼5-10 sec) are noticed in the 4th experiment when applying lifetime periods at every message. When applying a deterministic buffer capacity (3rd experiment), we notice mean response times from 10 to 16 sec. Finally, the 4th experiment, provides low mean response times (∼8-14 sec) for λ = 0.05 -2 and higher response times (∼17-49 sec) for λ = 2.5 -3.9.

B. Success Rate vs. Response Time

In order to study the trade-off between response times and success rates, we present the rates of successful message service in Fig. 5 for the same set of experiments. As expected, there is a 100% of message success rate in the 1st experiment since there are not message losses. On the other hand, the lowest rates are noticed when applying lifetime periods, especially when λ rates increase. In the 3rd experiment, low λ rates provide high message success rates (∼90%-100%) and higher λ rates decrease them (∼87%-77%). Finally, a stable message success rate (75%) can be achieved through the 4th experiment.

IV. CONCLUSION

Modeling the performance of IoT applications is a tedious task. Messaging in IoT may result to high message delays or message losses due to several QoS settings, such as message availability, intermittent connectivity, resource constrained devices, etc. In this paper, we define simulation queueing models that represent such QoS settings. Using the simulator we have developed, we evaluate the trade-off between response times and message success rates for different queueing models. By relying on our models, system designers are able to create queueing networks representing IoT interactions and analyze their performance or even tune it based on our experimental results.

In our future work, we intend to identify analytical models for queueing networks with lifetime periods and finite capacity buffers, as well as for the ON/OFF probability and loss models. Then, we plan to model the performance of IoT applications for domains such as emergency response scenarios.

  ON/OFF loss queue.

Fig. 2 :

 2 Fig. 2: Probabilistic and loss ON/OFF queues.

Fig. 3 :

 3 Fig. 3: Queues with event expirations and finite capacity.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Response Time of different ON/OFF models.