Amel Belaggoun
email: amel.belaggoun@cea.fr

Ansgar Radermacher
email: ansgar.radermacher@cea.fr

Valerie Issarny
email: valerie.issarny@inria.fr

ASLA: Adaptive System Level in AUTOSAR

This paper presents an initial approach towards making AU-TOSAR dynamic starting from the application layer down to the operating system level (task model and RTE)(i.e. extending AUTOSAR ECU Software architecture) and describes ASLA, which is the framework that provides taskslevel adaptation techniques in AUTOSAR.

INTRODUCTION AND MOTIVATION

Nowadays, the complexity of the next generation of automotive embedded systems such as Fully Electric Vehicles (FEVs) is increasing due to the growing number of functionality (more than 2.500 functions such as power trains, steering or braking systems, "X-by-wire" systems etc.). These systems are by nature, real-time. Moreover, most of them work under several resource constraints, due to cost, space and energy limitations. In addition, they are running in highly dynamic environments. Combining real-time features in tasks with dynamic behavior, together with cost and resource constraints may create new problems to be addressed in the design of such systems. Using the classical design approaches adopted in hard real-time systems, such as WCET analysis to guarantee timeliness, for example, is no longer acceptable in highly dynamic environments because it would waste resources and increase costs [START_REF] Buttazzo | Adaptive mechanisms for component-based real-time systems[END_REF]. Instead of allocating resources for the worst case, there is a need for smarter techniques to sense the current state of environment and react accordingly, which means, to cope with dynamic environments, automotive systems need to be adaptive; that is, they must be capable of changing their structure and/or their behavior to better reflect their current situation or in order to keep the system requirements at a desired level; if this is not possible, degrade it in a controlled way. Besides dealing with problems generated by adaptation in terms of meeting real-time constraints despite the system evolution; system reliability remains an important factor to be investigated within embedded systems due to their interactions with the physical word. Although, there has been much research on building reliable distributed real-time systems [START_REF] Balasubramanian | Middleware for resource-aware deployment and configuration of fault-tolerant real-time systems[END_REF][7], a trend towards more complex features in a system with cost and resource constraints poses a major challenges in developing such a system. In order to address the above challenges we need to study the feasibility of applying adaptation techniques in real-time embedded systems. Specifically, in automotive systems it would be very useful to support adaptation in order to increase the availability and reliability of software-based applications without additional hardware costs. Currently AUTOSAR [1]-The standard architecture for automotive systems-has no support for runtime adap-tation. Making it adaptive requires specific support at different levels of the software architecture. The most important component affecting adaptivity is the Operating System (OS), but some flexibility can also be introduced in the runtime environment(RTE). Therefore, we propose a layer called ASLA, which is used for task mapping, bandwidth allocation and adaptation for mixed-criticality distributed systems. The novelty of our approach is the capability to support the adaptation of applications with soft and hard real-time requirements (mixed) while respecting timing and safety requirements in AUTOSAR.

Overcoming AUTOSAR's limitations with ASLA: ASLA provides the ability to dynamically adapt the system, such as adding a new application or moving an existing application to a different ECU. In AUTOSAR, the system configuration is by design static: the AUTOSAR RTE is configured at design-time for specific ECUs and partly generated based on the requirements of the Software Components(SWC). A reconfiguration of the system, such as adding an application or moving an application from one ECU to another, cannot be done dynamically at runtime. ASLA extends AUTOSAR in two ways. It changes the scheduling policy from a fixed-priority (assigned to tasks at design time) to a dynamic preemptive policy based on EDF and CBS schedulers to guarantee mixed critical requirements. It also contains RTE extension that supports the runtime application migration between ECUs in response to both anticipated changes caused by the environment, such as network connectivity, as well as unexpected failures in both software and hardware. The deployment and re-configuration of an application onto ASLA is based on the work of [START_REF] Saraswat | Task mapping and bandwidth reservation for mixed hard/soft fault-tolerant embedded systems[END_REF], which uses a Tabu Search based meta-heuristic to search design space exploration to give the best task allocations and bandwidth allocation. Unlike the existing approaches [START_REF] Martorell | Improving adaptiveness of autosar embedded applications[END_REF][START_REF] Zeller | Timing constraints for runtime adaptation in real-time, networked embedded systems[END_REF], ASLA explicitly introduces the concept of runtime adaptation in mixed-criticality applications in automotive systems.

Contribution. This paper presents the design of a realtime adaptive distributed architecture, ASLA, to provide task-level adaptation techniques in AUTOSAR.

The remainder of this paper is organized as follows. Section 2 summarizes the background, the scope and requirements needed to understand the proposed approach. Section 3 presents our approach and the architecture of ASLA. Section 4 presents the related work. We state our conclusion and future work in Section 5.

BACKGROUND AND DEFINITIONS

Before delving into the approach, it is important to clarify the scope and the requirements that need to be addressed with our solution. We use Kiviat diagrams [START_REF]software empiricist[END_REF] to show visually the characteristics of our solution (Fig. 1) and the requirements we consider for adaptive automotive systems (Fig. 2). These diagrams provide developers of automotive software an easy way of viewing the characteristics of their applications. The dimensions represent axes of the Kiviat diagrams and characteristics of the dimensions represent the set of properties to be met by our solution (i.e. the red bullets).

The scope of our Research

The adaptation Model. The adaptation may typically be synchronous and/or asynchronous with respect to the execution of applications. In the synchronous case, the applications synchronize their execution, and new tasks (actions) are introduced only after all applications have finished performing the actions specified in the initial configuration. The schedulability analysis of the adaptation is thus not required, because there is no interference between tasks before and after the adaptation (i.e. in the new system configuration). On the other hand, in the case of asynchronous adaptation, all the applications start changing their configuration as soon as they receive an adaptation trigger without considering the behavior of the other applications. As a result, actions of the initial system configuration run concurrently with the new ones during the transition, which calls for schedulability analysis.

Promptness. Adaptation is well suited for systems that require reactive behavior, there is no need to wait until an idle period or slack before performing such a change as in Tindell's model [START_REF] Tindell | Mode changes in priority preemptively scheduled systems[END_REF] or at the end of cycles like in some approaches based on cyclic executive scheduling [START_REF] Real | Mode change protocols for real-time systems: A survey and a new proposal[END_REF]. The adaptation may be classified in three different categories in the time domain: (i) time-based adaptation: These are adaptations where we know the arrival time of the adaptation request in advance. (ii) event-driven adaptation: These are adaptations triggered by events rather than time. We don't know exactly when they happen.(iii) irregular eventdriven adaptation: These are adaptations when no prediction can be made about the arrival time of the adaptation request. An example of this type is a system fault. The system may change its configuration and migrate to a degraded one where not all the functionalities are provided. In our work, we consider all three types of adaptation.

Scheduling policy. The use of dynamic priorities scheme suits better with systems running in highly dynamic environments [START_REF] Real | Mode change protocols for real-time systems: A survey and a new proposal[END_REF] [START_REF] Abeni | Stochastic analysis of a reservation based system[END_REF].

Scheduling algorithm. Combining two scheduling mechanisms CBS and EDF has proven its efficiency to solve the problem of temporal isolation due to the integration of soft and hard real-time applications on the same platform [START_REF] Abeni | Stochastic analysis of a reservation based system[END_REF] [START_REF] Saraswat | Task mapping and bandwidth reservation for mixed hard/soft fault-tolerant embedded systems[END_REF].

Architecture. Distributed architecture. In our work a real-time distributed system is defined to be a system with multiple autonomous processing units (ECUs) cooperating together to achieve a common goal. We use the term distributed architecture to refer to loosely coupled architectures where message passing is required (full connectivity is assumed between ECUs i.e. each of the ECUs is connected to each other).

Timing requirement criticality. In our solution we are dealing with mixed-criticality systems, to the best of our knowledge no one has applied runtime adaptation considering both soft and hard real-time constraints in AUTOSAR.

Definition of our Adaptation Requirements

(R1) Timeliness: the timeliness requirements of an adaptation characterize the time constraints under which the adaptation is executed. Hard real-time constraints require the execution of adaptation within a firm deadline. Adaptations executed under soft real-time constraints minimize the adaptation execution and blackout time (which is the time the application is unavailable due to state transfer and reconfiguration). Unbounded adaptations are executed without any time bound [START_REF] Fritsch | Time-bounded adaptation for automotive system software[END_REF].

(R2) Consistency: preserving the system consistency and leaving the system under change in a correct state after adaptation are two major requirements that must be ensured when performing adaptation of the running system. Many adaptation approaches freeze the entities to be reconfigured into an adaptation safe state called quiescent state [START_REF] Kramer | Self-managed systems: An architectural challenge[END_REF].

(R3) Flexibility: The dynamic behavior of real-time systems requires executing applications with certain flexibility requirements in which the temporal properties and the number of applications vary during runtime. That means, the tasks of flexible application provide implementations that can adapt their execution to the available processing resources. These tasks have variable period and/or may demand variable WCET (i.e. stochastic [START_REF] Abeni | Stochastic analysis of a reservation based system[END_REF]). So executing flexible applications, prevents the use of an efficient static temporal partitioning of the processing time. A static temporal partitioning that lasts over the entire lifetime of a system would result in an oversized system. Hence, in order to efficiently use the processing time, the flexibility of applications and possible demand changes during runtime have to be considered during runtime analysis.

(R4) Adaptation trigger: adaptation can be triggered either internally due to the monitoring infrastructure or externally; requested from an outside entity, for example the user. Furthermore , adaptation triggers may arrive synchronously (i.e. at specific time) or asynchronously (i.e. with unknown arrival pattern). Timelines (R5) Type of adaptation: the type of adaptation defines what is being reconfigured. This type can be either, Resource adaptation(i.e.allocation of resources dynamically based on observed conditions), or Software adaptation which includes three categories: (i) structural, (ii) behavioral and (iii) parameter. The first one (i) structural adaptation, changes the actual architectural parts of an application, e.g., by removing a SWC, introducing a new one or replacing /updating an existing SWC with another newer version. The second category (ii) behavioral adaptation allows changes of the behavior of the application and the last category (iii) parameter adaptation involves modifying variable values that determine program behavior.

Type of adaptation

Flexibility

TOWARDS ADAPTIVE AUTOSAR

We present an approach that provides an adaptive solution for AUTOSAR [1]. The approach is based on monitoring various applications distributed on different ECUs to detect the need for adaptation at the application and the systemlevel ( i.e. the Basic SoftWare(BSW)) while maintaining the system consistency [START_REF] Kramer | Self-managed systems: An architectural challenge[END_REF], which means the system under adaptation must be left in a correct state after adaptation.

The ASLA Architecture

Fig. 3 provides an overview of ASLA's overall design. Every ECU that supports adaptation through ASLA layer consists of a real-time OS with EDF and CBS scheduling policies, an Adaptive SWC that is responsible for reconfiguring applications running on the system, RTE, and application layers. The real-time OS is responsible for HW abstraction, communication, scheduling and executing tasks in real-time. We assume that the underlying HW is a fail-silent system and the communication network is fault-tolerant. Our application layer consists of a set of SWCs (similar to AU-TOSAR's) and the new Adaptive SWC which can be distributed over several ECUs. The RTE provides a communication abstraction to SWCs. Unlike AUTOSAR, our RTE extension contains functions to support adaptation. These functions are managed by the Adaptive SWC (more precisely by the Reconfiguration Manager (RM)) which also communicates with the others Adaptive SWCs running on the different ECUs to make one of the adaptation actions such as: adding, deleting or updating application. When a new application is being added, the mapping between the application's SWCs and the ECUs is given to the RM, then each RM analyzes the mapping and renews the RTE's functions. The ASLA layer is composed of an Adaptive SWC (one on each ECU) and plugin offering a task execution container. This plugin enables any task launched on ASLA layer to be periodically executed. The adaptive component has a coordination-based architecture. One Adaptive SWC acts as a coordinator of the other Adaptive SWCs which are responsible for handling tasks on each ECU and monitoring a health vector. The latter contains all Non-Functional Requirements(NFRs) needed for the adaptation such as the ECU's processor utilization, resources, QoS, HW NFRs..etc. All operational ECUs compute their resources and processor utilization in form of a health vector at a fixed time period and share their health vector with each other. This provides each ECU a consistent view of the available resources and utilization on the other nodes. Since our Adaptive SWC has a coordination based-architecture, we define a management protocol between the different Adaptive SWCs inspired by [START_REF] Mitzlaff | Enabling mode changes in a distributed automotive system[END_REF]. However, we differ from them in the sense that our protocol is much more simpler and it is used for managing the process of an adaptation in distributed real-time systems. In our protocol, all the Adaptive SWCs including the coordinator broadcast messages to each other. The coordinator can detect the failure of the others Adaptive SWCs by a lack of heartbeat messages. The major components of ASLA are described below.

A. The Adaptive SWC. As illustrated in Fig 3, an Adaptive SWC is composed of a monitor,a Mapping Manager (MM) and a Reconfiguration Manager(RM). The monitor is responsible for monitoring events that trigger the adaptation. The MM offers a dynamic deployment of tasks on the ECUs and the RM can automatically reconfigure tasks inside/or between the different ECUs :

• The Monitor. The monitor periodically sends messages to other ECUs in the system via the network 1 . The monitor allows ASLA to agree on the availability of each ECU. Any adaptation trigger received by the application during its execution may invoke the monitor which sends a message to the RM in order to adapt the application. The loss of a message for two consecutive cycles means that the ECU is no longer alive and the adaptation needs to be triggered to accommodate the desired changes.

• The Mapping Manager. The MM offers an automatic deployment of tasks on ECUs. We use TSMBA (Tabu search Mapping and Bandwidth Allocation) [START_REF] Saraswat | Task mapping and bandwidth reservation for mixed hard/soft fault-tolerant embedded systems[END_REF] as a base line algorithm for our work to do the allocation which provides a comprehensive solution that allocates mixed critical application to a distributed heterogeneous architecture and reserves processor bandwidth for guaranteeing timing requirements. We propose the extension O-TSMBA (Operational chains-TSMBA) a variant of TSMBA that supports task dependencies (i.e. pipeline tasks model). The MM takes as input the application description (an initial system configuration file) and changes the current mapping when it's necessary to do so. Changes of the allocation can occur due to the adaptation or in case of one or several ECUs failures.

• The Reconfiguration Manager. The RM is a sporadic task that gets triggered upon the reception of an adaptation trigger (requests for adding new tasks, requests for migrating failed tasks/and or failed ECUs, replacement of tasks with an improved version and removing tasks). B. ASLA Plugins. All applications will run on the top of ASLA plugins. ASLA plugins support the mechanisms for task reconfiguration and bandwidth allocation (i.e.TSeRBA algorithm "Tabu Search Reconfiguration and Bandwidth Allocation") and also enables tasks to have guaranteed and protected access to required processing resources during reconfiguration in a timely manner. 1 The Network is beyond the scope of this paper. We assume a synchronous communication network

ASLA Development Process

Our process for developing automotive software in compliance with AUTOSAR standard is shown in Fig. 4. It starts with an application description, in terms of a SWC architecture, dependencies between SWCs, real-time constraints, HW resource requirements and other information needed in the vehicle. In our approach we are interested in mixed critical applications with soft and hard real-time requirements. We consider that each SWC contains one runnable and is represented by one AUTOSAR task. We construct an operational chain OP which is composed of periodically executing runnables generating data and events regularly that flow through multiple runnables (i.e. they are connected by data flow or/and control flow). OP correspond to the AUTOSAR execution model (see Fig4-w g 1). We distinguish chains with soft and hard timing requirements (OP Soft and OP Hard ). We assume that the initial mapping of runnables to AUTOSAR tasks is given at design time similarly to AU-TOSAR and all runnables are executed periodically within the context of an AUTOSAR Task (see Fig4-w g 2 ). At this step we will have an initial solution which is not necessary schedulable and it is given as input to the task mapping algorithm O-TSMBA which we designed for the operational chain model. The output of O-TSBMA algorithm will be used as input for TSeRBA (i.e. a configuration schedulable and tagged optimized (see Fig4-w g 3 )). TSeRBA algorithm is used for task mapping, bandwidth allocation and reconfiguration for mixed-criticality distributed systems. The runtime support to AUTOSAR will be realized by TSeRBA as it allows the dynamic allocation of SWCs with both hard and soft real-time constraints as well as supports the insertion, the deletion and the migration of SWCs at runtime (see Fig4-w g 4).

RELATED-WORK

Runtime adaptation in real-time systems have been extensively studied in [START_REF] Becker | Model-based extension of autosar for architectural online reconfiguration[END_REF][START_REF] Zeller | Timing constraints for runtime adaptation in real-time, networked embedded systems[END_REF][START_REF] Kim | SAFER: system-level architecture for failure evasion in real-time applications[END_REF]. Unlike existing works, ASLA provides a framework to support runtime adaptation with taking into account schedulability analysis and task allocation for mixed-criticality applications in AUTOSAR. Except SAFER [START_REF] Kim | SAFER: system-level architecture for failure evasion in real-time applications[END_REF] which can be used as a complementary solution for our proposal for the dependability point of view. Another well researched topic in real-time systems-mode change concept, which has a close relationship with the adaptation. In order to guarantee timing requirements in the presence of changes in the system several approaches have been proposed and focused on mode change protocols (see the survey of [START_REF] Real | Mode change protocols for real-time systems: A survey and a new proposal[END_REF] for more details), [START_REF] Nelis | Global-edf scheduling of multimode real-time systems considering mode independent tasks[END_REF], which can be potentially used as an extension to ASLA. The authors of [START_REF] Nolte | A synchronization protocol for temporal isolation of software components in vehicular systems[END_REF] proposed SIRAP, a protocol for synchronization in a hierarchical realtime scheduling framework. Their approach is relevant to our case because we are using a server-based technology to schedule AUTOSAR SWCs. However, we have a different focus; we are tackling the challenge of making runtime adaptation in AUTOSAR.

CONCLUSION AND FUTURE WORKS

Runtime adaptation in embedded real-time systems is a topic expected to grow in the coming years, gaining particular moment in the context of designing FEVs. Yet there is a lack of techniques and tools for performing such adaptation. We presented ASLA, a novel framework that supports task-level reconfiguration features in AUTOSAR. We have built an experimental platform using three ARM-based STM32F4Discovery boards, an Open-source AUTOSAR implementation, ERIKA-OS [2] served as the BSW. A CAN bus communication was established between the three ECUs and currently we are focusing on the implementation of the algorithms behind ASLA framework to demonstrate the theoretical ideas. As future work, ASLA will be validated by means of a real case study from the SafeAdapt project. We will measure the overhead of the adaptation mechanisms, in particular the impact on timeliness.

Figure 1 :

 1 Figure 1: The scope of our Research.

Figure 2 :

 2 Figure 2: Requirements for Runtime Adaptation. (R5) Type of adaptation: the type of adaptation defines what is being reconfigured. This type can be either, Resource adaptation(i.e.allocation of resources dynamically based on observed conditions), or Software adaptation which

Figure 3 :

 3 Figure 3: The ASLA architecture.The ASLA layer is composed of an Adaptive SWC (one on each ECU) and plugin offering a task execution container. This plugin enables any task launched on ASLA layer to be periodically executed. The adaptive component has a coordination-based architecture. One Adaptive SWC acts as a coordinator of the other Adaptive SWCs which are responsible for handling tasks on each ECU and monitoring a health vector. The latter contains all Non-Functional Requirements(NFRs) needed for the adaptation such as the ECU's processor utilization, resources, QoS, HW NFRs..etc. All operational ECUs compute their resources and processor

Figure 4 :

 4 Figure 4: ASLA Development Process

An initial configuration which is schedulable & tagged <optimized> i.e. AUTOSAR Tasks i.e. AUTOSAR execution model (runnables and messages) CBS CBS RTE: The proposed tasks adaptation mechanisms … EDF CBS CBS CBS ECU HW RTE Adaptive SWC SWC (app1) SWC (app2) EDF CBS CBS CBS ECU HW RTE: The proposed tasks adaptation mechanisms RTE Adaptive SWC SWC (app1) SWC (app2) EDF CBS CBS CBS ECU HW RTE Adaptive SWC SWC (app1) SWC (app2) Monitor Mapping Manager Reconfiguration Manager ... ... ... Plugin Plugin Plugin Plugin Plugin Plugin Plugin Plugin Plugin Network abstraction Adaptive SWC Plugin RTE: The proposed tasks adaptation mechanisms ASLA layer ASLA layer ASLA layer RTE: The proposed tasks adaptation mechanisms

			Application
			Description
			Construct
	Level 1 mapping	1	operational
	@ design time		chains
		2	Construct tasks
		3	Deployment
	Level 2 mapping and
	reconfiguration	4	TSeRBA
	@ run-time
		5	Reallocation & reconfiguration

Acknowledgments: We would like to thank Robert Davis for the interesting discussions about this work.


Yoann Roth
email: yoann.roth@cea.fr

Jean-Baptiste Doré
email: jean-baptiste.dore@cea.fr

Laurent Ros
email: laurent.ros@gipsa-lab.grenoble-inp.fr

Vincent Berg
email: vincent.berg@cea.fr

A Comparison of Physical Layers for Low Power Wide Area Networks

Keywords: Low Power Wide Area, LPWA, Internet of Things, PHYlayer, FSK, Turbo Code, Low Rate, Low SNR, Machine-to-Machine

In Low Power Wide Area networks, terminals are expected to be low cost, low power and able to achieve successful communication at long range. Communication should be low rate compared to cellular mobile networks. Most of the current technologies dedicated to Machine-to-Machine communication rely on the use of a spreading factor to achieve low levels of sensitivity. To contain power consumption at the terminal side, the cost of complexity should be paid by the receiving side. We propose to use turbo processing schemes as potential physical layers. We compare these schemes to a standard and an industrial solution, and show that a significant gain in sensitivity can be achieved, and very energy efficient scheme can be designed by mixing turbo processing and orthogonal modulation concepts.

Introduction

The Internet-of-Things (IoT) is a vast concept, where every object is expected to be connected. Amongst the foreseeable applications for these Machine-to-Machine (M2M) communications, many require the design of a new Low Power Wide Area (LPWA) network [START_REF] Rebbeck | Low-powered wireless solutions have the potential to increase the M2M market by over 3 billion connections[END_REF]. While recent advances (3G, 4G) sought for higher data rate and spectral efficiency, this kind of network is expected to be low rate, low power and long range. Packet sizes should be from a few bytes to hundreds maximum. Terminal power consumption and sensitivity working levels are critical issues, relying on the choice of an efficient and robust physical layer. We consider a star network topology, where each node is connected to a base-station. We focus on the uplink communication.

To achieve low sensitivity levels, a widely used technique is the low complexity repetition code. Indeed, repeating the information by a Spreading Factor (SF) λ at the transmitter can offer a gain of 10 log 10 (λ) dB in sensitivity at the receiver side. This technique offers a good compromise between sensitivity and complexity and many systems rely on it, as for instance in Direct Sequence Spread Sprectrum (DSSS). Another natural choice for low data rate and low power applications is the use of M -ary orthogonal modulation. Increasing the size of the alphabet gives a gain in energy efficiency, while reducing the spectral efficiency [START_REF] Proakis | Digital Communications 3rd Edition, ser. Communications and signal processing[END_REF]. This property can be extended to Orthogonal Sequence Spread Spectrum (OSSS), and using high orders of alphabet can bring low levels of sensitivity.

Some solutions have been proposed to answer the need of a new LPWA network. A recent proposal [START_REF]LoRa Alliance[END_REF][START_REF] Seller | Low power long range transmitter[END_REF] is based on the combination of OSSS with a low complexity block code. Also, efforts have been made to design the standard 802.15.4k using the DSSS technique combined with a convolutional code [5], which can be considered for LPWA networks. These two schemes are low complexity and straightforward solutions. However, since low power consumption is expected mainly at the terminal side, the cost of complexity can be paid by the receiver side, i.e the base station. More elaborate decoding schemes can then be used, such as turbo processing [START_REF] Berrou | Near Shannon limit error-correcting coding and decoding: Turbo-codes. 1[END_REF], and lead to a better energy efficiency of the system. We consider the widely known Turbo Code (TC) standardized by the Universal Mobile Telecommunications System (UMTS) [START_REF]LTE Evolved Universal Terrestrial Radio Access (E-UTRA): Multiplexing and Channel Coding[END_REF][START_REF] Valenti | The UMTS Turbo Code and an Efficient Decoder Implementation Suitable for Software Defined Radios[END_REF] combined with a common modulation and a SF as a potential physical layer for LPWA. Our recent work [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF] pointed out that if the repetition scheme allows a decrease of the required Signal-to-Noise Ratio (SNR) of a specific modulation, it does not increase the energy efficiency of the system (which can be expressed in terms of energy per bit to noise spectral density ratio, denoted E b /N 0 ). This is a consequence of the information rate reduction implied by the use of the SF. Redundancy should be used more efficiently by the receiver, by means of turbo processing. The scheme proposed in [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF], dedicated to LPWA networks, combines orthogonal modulation with turbo processing, and achieves very promising level of performance.

In this paper, we propose a comparison of four different schemes: the 802.15.4k standard [5], OSSS with a Hamming code, the TC UMTS from [START_REF]LTE Evolved Universal Terrestrial Radio Access (E-UTRA): Multiplexing and Channel Coding[END_REF][START_REF] Valenti | The UMTS Turbo Code and an Efficient Decoder Implementation Suitable for Software Defined Radios[END_REF] with a common modulation (Binary Phase Shift Keying (BPSK) or M -ary orthogonal modulation), and the Turbo-FSK scheme proposed in [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF]. All of them perform at a low spectral efficiency, aim a high energy efficiency, and would be suitable for uplink communication in the context of LPWA networks. We focus on the physical layer, and assume perfect synchronization and Additive White Gaussian Noise (AWGN). All architectures are described, and simulations are performed. To the best of our knowledge, comparison of these potential LPWA networks have never been done. Also, the use of the UMTS Turbo Code and a modulation have never been considered before as a potential uplink scheme for LPWA.

The paper is organized as follows. The different schemes and their architectures are introduced Section 2. Performance comparison are presented Section 3. In Section 4, results are discussed and Section 5 concludes the paper.

System Model

This section is dedicated to the description of the four conceivable physical layers previously mentioned. For each one, a short introduction is followed by the presentation of the block diagram of the transmitter. As receivers is rarely standardized, a receiver scheme has been proposed for each physical layer.

IEEE 802.15.4k

The IEEE 802.15.4k is a standard for local and metropolitan area networks, and is part of the Low-Rate Wireless Personal Area Networks (LR-WPANs) [5]. It aims at low energy critical infrastructure monitoring networks. This standard supports three physical layer (PHY) modes: DSSS with BPSK or Offset-Quadrature Phase Shift Keying (O-QPSK), or Frequency Shift Keying (FSK). DSSS with BPSK modulation is adapted to more constrained situations, and will be presented here. Standard specifications allow the use of a spreading factor value λ from 16 to 32768.

Block diagrams for transmitter and receiver are given in Figure 1. The transmitter is composed of a Forward Error Correction (FEC) block, defined to be the convolutional code of rate 1/2, with generators polynomial [171 133] (in octal) and constraint length k = 7. After encoding, interleaving is done to ensure diversity at the reception side. Differential encoding is applied, and encoded data is "repeated" by the use of a binary direct sequence of size λ. BPSK modulation is then used to transmit the signal through the channel. The normalized spectral efficiency of the physical layer, expressed in bits/s/Hz, is

η 1 = 1 2λ . (1)
The receiver executes the reverse operations of the transmitting side. After despreading the signal (executing the mean weighted by the elements of the binary direct sequence), we use a soft differential BPSK decoder. Knowing how data was interleaved at the transmitter side, the deinterleaving operation allows the soft input Viterbi decoder to decode the data and to finally retrieve the information bits.

Orthogonal modulation with block code

It is a known fact that orthogonal modulations reach channel capacity for large size of alphabet [START_REF] Proakis | Digital Communications 3rd Edition, ser. Communications and signal processing[END_REF]. An original solution relies on the use of orthogonal modulations with a simple channel code [START_REF] Seller | Low power long range transmitter[END_REF][START_REF]SX1272 from Semtech, datasheet[END_REF]. Spreading operation is replaced by increasing the alphabet size, thus lowering the spectral efficiency but also increasing the energy efficiency. As we assumed ideal synchronization and AWGN channel, performance does not depend on the orthogonal alphabet choice (e.g FSK, OSSS, Pulse Position Modulation, . . . ).

The transmitter and receiver architectures are given in Figure 2. After applying a simple Hamming (7, 4, 3) block code, coded bits are interleaved and then mapped to an orthogonal alphabet of size λ. The spectral efficiency is defined as

η 2 = log 2 (λ) λ 4 7 . (2)
At the receiving side, soft demodulation is performed. The output is the log likelihood ratio (LLR) of the bits. It is defined as

L(b n ) = log p(b n = 1) p(b n = 0) . (3)
The principle of soft decoding is to compute the A Posteriori Probabilities (APP) of the codewords, i.e for each codeword the probability of having a given codeword knowing the observation, and then compute the LLR of the information bits, knowing the APP for all the possible codewords. The elements of the orthogonal alphabet are denoted c i with i ∈ [0, λ -1] and

c i = [c i 0 , . . . c i λ-1 ]. For each c i , the associated information word is denoted b i = [b i 0 , . . . b i log 2 (λ)-1
]. We consider one noisy received codeword y, and d its associated decoded information word. The APP can be expressed as

p(c i |y) = p(y|c i ) p(c i ) p(y) , (4)
after applying Bayes' law. p(c i ) is the a priori probability of having the codeword, which here is equal to 1/λ. As y and c i are vector with λ elements, the APP becomes:

p(c i |y) = 1 λp(y) λ-1 m=0 p(y m |c i m ). (5)
Considering the complex AWGN case,

p(y m |c i m ) = 1 2πσ 2 exp - 1 2σ 2 y m -c i m 2 , ( 6
)
where σ is the noise variance and . is the Euclidean norm. y m and c i m are complex numbers, hence (5) can be expressed p(

c i |y) = C λp(y) exp 1 σ 2 λ-1 m=0 Re y m .c i m , ( 7
)
where C is a constant that will be canceled out in further computations, and c i m the complex conjugate of c i m . LLR of the information bits can be computed using the expression

L(d n |y) = log λ-1 i=0, b i n =1 p(c i |y) -log λ-1 i=0, b i n =0 p(c i |y), (8)
where the index b i n = 1 (resp. b i n = 0) for the first sum (resp. the second sum) signifies that the sum is done over the codewords c i that encode an information word for which the value of bit b i n is 1 (resp. 0). The factor before the exponential in ( 7) is suppressed at this step.

Once LLR have been computed, deinterleaving is done. Soft Hamming decoding can then be executed to get the decoded information bits.

Standard modulation with turbo coding

3G and 4G standards both rely on a turbo code for forward error correction [START_REF]LTE Evolved Universal Terrestrial Radio Access (E-UTRA): Multiplexing and Channel Coding[END_REF]. The TC used is a parallel concatenation of two identical recursive systematic codes of rate 1/2, making the code rate equal to 1/3 (or slightly less when the termination bits needed to close the trellis are considered). The constituent code has a constraint length of k = 4, and its feedback and generator polynomials are respectively [1315] (in octal). This channel code allows any system to achieve an energy efficiency gain relatively high, thanks to iterative soft decoding at the receiver side. Complexity at the transmitter side is rather low. We consider this code combined with a modulation (BPSK or M -ary FSK) as a potential physical layer for LPWA network, or at least as an interesting element of comparison.

The architecture of the transmitter is depicted Figure 3. The Parallel-to-Serial Conversion has three inputs : the information bits, the parity bits of the first Recursive Systematic Coder (RSC), and the parity bits of the second RSC which input has been interleaved. After serial conversion, a SF of value λ is applied, and then encoded data is modulated. The spectral efficiency of this physical layer is

η 3 = η mod 3λ , (9)
where η mod is the normalized spectral efficiency of the modulation used, defined as

η mod =    1 if BPSK is used, log 2 (M ) M
if M -ary orthogonal modulation is used.

The receiver suggested is described Figure 4. After the noisy channel, the SF Recombination module averages λ consecutive symbols. Signal is then demodulated, i.e LLR are computed. Serial-to-Parallel Conversion will give the systematic bits (i.e uncoded information bits), and parity bits of each encoder. Each decoder is then fed by the systematic bits (which are interleaved for the second decoder) and parity bits corresponding to each RSC. Maximum A Posteriori (MAP) rule is used to decode the trellis, using the Bahl, Cocke, Jelinek and Raviv (BCJR) algorithm [START_REF] Bahl | Optimal decoding of linear codes for minimizing symbol error rate (corresp.)[END_REF]. The max-log approximation can be used, implying a small loss in performance for an interesting reduction of complexity. For each decoder, the output of the other decoder is used as a priori input after interleaving. Several iterations of this process are performed, and a decision can be made at the end of each iteration.

Turbo-FSK

A recent scheme has been proposed, using orthogonal modulations and turbo decoding at the receiver side [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF]. This scheme allows for a better use of redundancy than in pure repetition scheme, and achieves an interesting energy efficiency gain. The transmitter low complexity makes it suitable for uplink communication in the LPWA context. Transmitter architecture is given in Figure 5. The structure is composed of λ stages, each one encoding an interleaved version of the input bits. Information bits are gathered into P groups of r bits, and Convolutional-FSK encoding is applied, as described Figure 6. For every group of r bits, parity is computed and accumulated in the memory. The r + 1 bits are then mapped to a codeword of the FSK alphabet, which size is M = 2 r+1 . Thanks to the accumulator, every consecutive symbol is linked to the previous one. The output of each encoder is then a set of P FSK codewords: this scheme mixes coding and modulation in the same process. A Parallel-to-Serial Conversion is then done to send the FSK codewords through the channel. The spectral efficiency is defined by

η 4 = log 2 (M ) -1 λM , (10)
The receiving side, depicted in Figure 7, consists of Serial-to-Parallel Conversion, to reconstruct the λ stages emitted at the first place. A Soft FSK Detector is used to determine the probabilities of each possible codewords. This step can be done using the Fast Fourier Transform (FFT) algorithm. These probabilities are then fed to the APP decoder, which will use them as channel observation, while output of the other decoders will be used as a priori information. BCJR algorithm is used to decode the trellis. Orthogonality of the transmitted codewords is one of the key feature of this scheme, which offers a great performance gain compared to simple repetition schemes. Detailed explanations about the transmitter and the receiver's computations can be found in [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF].

Simulation Results

To perform a comparison between the schemes presented in Section 2, we perform simulations for the AWGN channel, with coherent reception. Random interleavers are used, except for the TC UMTS where the interleaver specified by the 3G standard is considered. Packet size is set to 1024 bits (or 128 bytes), which is rather small but appropriate for the LPWA context. For iterative schemes, 10 decoding iterations are performed, since performance improvement is rather small for larger numbers of iterations. MAP algorithm is used.

For the sake of fairness, the different schemes should be compared with an equal normalized spectral efficiency η. This physical feature can be adjusted thanks to the definition of the parameter λ, which is related to either a repetition factor or an alphabet size. The normalized spectral efficiency of a configuration of the Turbo-FSK (32-FSK with λ = 4 stages, with η 4 1/32 = 3.125 • 10 -2 ) is taken as reference. The adjusted spectral efficiencies, the parameters λ, the modulations, the FEC chosen and the binary code rate are summarized in Table 1. For the third scheme, two different modulations are chosen: BPSK (physical layer 3) and 32-FSK (physical layer 3bis), hence the total of 5 physical layers. To obtain the SNR, we consider the straightforward conversion

SNR = η E b N 0 . (11)
Figure 8 shows the Bit Error Rate (BER) performance versus SNR for the 5 physical layers. All schemes have interesting performance at a very low SNR, and will work at low sensitivity levels, which was the goal to achieve. Scheme 1 shows the typical performance of a convolutional code. Scheme 2 uses the Hamming code, which is less powerful than the convolutional code of scheme 1, but offers better performance when combined with a relatively large size of alphabet (128). Performance clearly shows the gain of turbo processing. Scheme 3, corresponding to the Turbo Code UMTS with BPSK modulation and λ = 11, achieves a BER of 10 -5 for a SNR of -14.2dB, hence outperforming scheme 2 by 4.5dB. Scheme 3bis, for which 32-FSK modulation is used, suffers from a performance loss compared to 3. Scheme 4 outperforms all the other schemes, showing the beneficial impact of mixing convolutional encoding and FSK modulation as proposed in [START_REF] Roth | Turbo-FSK: A New Uplink Scheme for Low Power Wide Area Networks[END_REF]. It should be noted that both schemes 3bis and 4 use the same size of orthogonal alphabet, with yet a gain greater than 1.3dB at a BER of 10 -4 for the Turbo-FSK scheme. Packet Error Rate (PER) performance versus the SNR is depicted Figure 9. Schemes involving turbo processing are again shown to be far more efficient than the other schemes, with a gain superior to 3.5dB for a PER of 10 -3 . The ultimate Shannon's limit [START_REF] Shannon | A mathematical theory of communication[END_REF] is the maximal transmission rate with arbitrarily small bit-error probability, for a given SNR and a given bandwidth. Comparing schemes to this limit thus shows how efficiently the channel resource is used. The limit, expressed in terms of normalized spectral efficiency versus E b /N 0 , is represented in Figure 10, along with the performance of the 5 physical layers that are here presented. M -ary orthogonal modulation is also added for reference. Since we choose a specific value of normalized spectral efficiency, all performance corresponding to the schemes are approximately on the same horizontal line. Scheme 4 is the closest scheme to the E b /N 0 Shannon's limit as it is only 2.3dB away, showing it to be the scheme with the best energy efficiency. We also give the performance of scheme 3 but without any repetition (thus λ = 1). As repetition does not influence performance in E b /N 0 , both schemes have the same energy efficiency (same vertical line). Comparing scheme 4 to the uncoded M -ary orthogonal modulation, performance shows that turbo processing is a

Fig. 9. PER performance comparison in SNR. Packet size is 1024 bits. For schemes using turbo processing (plain curves), 10 iterations have been performed and the block size is set to 1024 bits. Spectral efficiency is roughly equal to 3.113 • 10 -2 bits/s/Hz.

way to get closer to Shannon's limit with a reasonable size of alphabet, hence a reasonably low spectral efficiency. The given E b /N0 values are the ones for which the BER is 10 -5 . For schemes using turbo processing, 10 iterations have been performed and the block size is set to 1024 bits.

Discussion

As presented in Section 3, schemes involving turbo processing have clearly the potential to reach low sensitivities. All gains presented in the previous section can be associated to the gain in sensitivity: using scheme 4 instead of scheme 1 with the same bandwidth will lead to the same data rate but sensitivity will be 5dB lower. Same performance can thus be achieved with a transmitted power divided by 3.1 (in linear scale) or distance between transmitter and receiver increased by a factor 1.7 (considering a free space path loss). In an other way, for a same level of sensitivity, schemes based on repetition will have to use a higher SF, thus a longer signal for an equivalent bandwidth. This can imply harder constraints for synchronization, in terms of phase noise or channel coherence time.

The performance improvements are done at the expense of an increased complexity. However, the cost of complexity of the turbo processing-based schemes is paid by the receiver side, i.e the base station, which can be considered to have unlimited resource. This is not often suggested for the LPWA context, but more advanced processing can be considered, leading to better performance.

As mentioned in Section 1, we assumed perfect synchronization. Since every communication system needs to be synchronized to work properly, this aspect must be taken into account. This is one of the major issues concerning low SNR communication, that still needs to be solved.

Conclusion

The design of a new LPWA network is a major issue addressed in the context of the IoT. The uplink scheme is critical, as terminals must have low power consumption and high reliability. We compared several schemes as potential candidates, and have shown that if repeating the information is a good way to lower the sensitivity, more complex systems can allow for higher energy efficiency gain. Simulation results show that turbo processing is an efficient way to achieve a gain in sensitivity. Amongst the studied turbo schemes, Turbo-FSK is the most efficient, and demonstrates the positive impact of mixing turbo processing, Mary orthogonal modulation and low spectral efficiency. However, synchronization for this ranges of SNR and complexity remain open questions.

Fig. 1 .

 1 Fig. 1. Standard IEEE 802.15.4k: (a) transmitter architecture, (b) receiver architecture.

Fig. 2 .

 2 Fig. 2. Orthogonal modulation and Hamming code: (a) transmitter architecture, (b) receiver architecture.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Modulation and Turbo Code UMTS: transmitter side.ChannelSpreading Factor Recombination

Fig. 5 .Fig. 6 .

 56 Fig. 5. The Turbo-FSK transmitter.

Fig. 7 .

 7 Fig. 7. The Turbo-FSK receiver.

Fig. 8 .

 8 Fig.8. BER performance comparison in SNR. For schemes using turbo processing (plain curves), 10 iterations have been performed and the block size is set to 1024 bits. Spectral efficiency is roughly equal to 3.113 • 10 -2 bits/s/Hz.

Fig. 10 .

 10 Fig. 10. Performance comparison considering spectral efficiency versus E b /N0.The given E b /N0 values are the ones for which the BER is 10 -5 . For schemes using turbo processing, 10 iterations have been performed and the block size is set to 1024 bits.

Table 1 .

 1 Parameters used for comparison. Termination bits for convolutional codes are considered when computing the normalized spectral efficiency values, implying a lower value than expected. For Turbo-FSK, no binary code rate can be defined, since coding and modulation are mixed.

	PHY-layer	1	2	3	3bis	4
	Modulation	DBPSK	128-Orthog	BPSK	32-FSK	32-FSK
	FEC	CC [171 133] Hamming	TC [13 15] TC [13 15] Turbo-FSK
	Binary code-rate	1/2	4/7	1/3	1/3	-
	λ	16	128	11	2	4
	η (•10 -2 )	3.125	3.120	3.019	2.594	3.113


Numerical experiments Artificial validation case

Scattering of a plane wave by a Au-nanowire (1) Ω ∈ R d , d = 2, 3, t f is the final physical time. E, H, and J respectively denote the electric field, the magnetic field and the dipolar current vector. 2 DGTD formulations Weak variation Space-discrete formulation 3 Reduced-order DG formulations POD-DGTD formulations Krylov subspace based solver Krylov-subspace based operator exponential approach Krylov-subspace based Padé approximation approach 4 Numerical experiments Artificial validation case Scattering of a plane wave by a Au-nanowire 1D: d i = p + 1. 2D: d i = 1 2 (p + 1)(p + 2). 3D: d i = 1 6 (p + 1)(p + 2)(p + 3).

1D:

Thus, it becomes an important problem how to get a low dimensional system that has the same response characteristics as the original system, while sufficiently accurate numerical solutions.

far less storage requirements and much lower evaluation time.

Krylov-subspace based operator exponential approach Krylov-subspace based Padé approximation approach

Numerical experiments Artificial validation case

Scattering of a plane wave by a Au-nanowire Summary Notations K i : an element of the mesh.

N Ω : the set of indices of the mesh elements. F I h : the union of all interior interfaces. F B h : the union of all boundary interfaces. F ij : the common face of the K i and K j , K j is a adjacent element. ν i := {j|K i ∩ K j = ∅}: the set of indices of the elements which are neighbors of K i .

Notations

K i : an element of the mesh. N Ω : the set of indices of the mesh elements. F I h : the union of all interior interfaces. F B h : the union of all boundary interfaces. F ij : the common face of the K i and K j , K j is a adjacent element.

Introduction Model Problem

Time-domain Maxwell-Drude equations

               ∂H ∂t + curlE = 0 in Ω × (0, t f ], ε ∞ ∂E ∂t -curlH + J = 0 in Ω × (0, t f ], ∂J ∂t + γ d J -ω 2 d E = 0 in Ω s × (0, t f ], E(x; 0) = E 0 (x), H(x; 0) = H 0 (x) in Ω, J(x; 0) = J 0 (x) in Ω s .
(1)

Time-domain Maxwell-Drude equations

               ∂H ∂t + curlE = 0 in Ω × (0, t f ], ε ∞ ∂E ∂t -curlH + J = 0 in Ω × (0, t f ], ∂J ∂t + γ d J -ω 2 d E = 0 in Ω s × (0, t f ], E(x; 0) = E 0 (x), H(x; 0) = H 0 (x) in Ω, J(x; 0) = J 0 (x) in Ω s .
(1)

Time-domain Maxwell-Drude equations

               ∂H ∂t + curlE = 0 in Ω × (0, t f ], ε ∞ ∂E ∂t -curlH + J = 0 in Ω × (0, t f ], ∂J ∂t + γ d J -ω 2 d E = 0
in Ω s × (0, t f ], E(x; 0) = E 0 (x), H(x; 0) = H 0 (x) in Ω, J(x; 0) = J 0 (x) in Ω s .

Introduction Model Problem

Boundary condition

Hard-wall condition [Schmitt N, Scheid C, Lanteri S, et al, JCP, 2016.]:

n • J = 0 on ∂Ω s .
(2)

Perfect electric conductor condition (PEC):

n × E = 0 on Γ m .

(3)

First order Silver-Müller absorbing boundary condition (ABC):

n × E + zn × (n × H) = n × E inc + zn × (n × H inc ) on Γ a . (4)

Boundary condition

Hard-wall condition [Schmitt N, Scheid C, Lanteri S, et al, JCP, 2016.]:

n • J = 0 on ∂Ω s . (2)
Perfect electric conductor condition (PEC):

n × E = 0 on Γ m .

(3)

First order Silver-Müller absorbing boundary condition (ABC):

n × E + zn × (n × H) = n × E inc + zn × (n × H inc ) on Γ a . (4) Γ m ∩ Γ a = ∅, Γ m ∪ Γ a = ∂Ω, ∅ is an empty set.
n denotes the unit normal vector pointing outward to ∂Ω.

E inc and H inc are the incident fields.

z = µ/ε.

Motivation

At present, most of numerical solvers for (1) with ( 2)-( 4) can be summarized as: FDTD, FETD, FVTD, DGTD, etc. The DGTD methods has many advantages compared with the previous numerical methods: the use of unstructured possibly non-conforming meshes.

easily obtained high order accuracy.

hp -adaptivity. natural parallelism.

...

Introduction

Motivation

Motivation

DOFs of the DGTD method: N dof = nt i=1 d i . n t is the number of elements. d i is the number of DOFs on element i.

Introduction

Motivation

Motivation

DOFs of the DGTD method:

N dof = nt i=1 d i .
n t is the number of elements. d i is the number of DOFs on element i. {v}

Outline

Introduction

F ij = 1 2 (v j + v i )
where v j and v i respectively stand for the traces of v on F ij from the interior of K j and K i . The finite dimensional subspace can be defined by:

V p h := {v ∈ (L 2 (Ω)) 3 | v |K i ∈ (P p (K i )) 3 , ∀i ∈ N Ω }.
The restrictions of V p h in Ω s is denoted by V

Weak variation

The DG method seeks an approximate solution (E h , H h , J h ) in the space

V p h × V p h × V p h that should satisfy for all K i in Ω h :      ( ∂H h ∂t , v) K i + (curlE h , v) K i = 0, ∀v ∈ V p h , (ε ∞ ∂E h ∂t , v) K i -(curlH h , v) K i + (J h , v) K i = 0, ∀v ∈ V p h , ( ∂J h ∂t , v) K i + (γ d J h , v) K i -(ω 2 d E h , v) K i = 0, ∀v ∈ V p h . (5)
Using appropriate Green's formulas and replacing the boundary terms by numerical traces E * h and H * h in (5), which are also known as numerical fluxes, we have:

     ( ∂H h ∂t , v) K i + (E h , curlv) K i -E * h × n, v ∂K i = 0 ∀v ∈ V p h , (ε ∞ ∂E h ∂t , v) K i -(H h , curlv) K i + (J h , v) K i + H * h × n, v ∂K i = 0 ∀v ∈ V p h , ( ∂J h ∂t , v) K i + (γ d J h , v) K i -(ω 2 d E h , v) K i = 0 ∀v ∈ V p h . (6)

Space-discrete formulation

Proper choice of numerical traces E * h and H * h is essential for the correctness and the convergence of the scheme. For the sake of simplicity, we choose to use a fully centered numerical flux, i.e.

E * h = {E h } F ij , H * h = {H h } F ij , ∀i ∈ N Ω , ∀j ∈ ν i .
Replacing surface integrals using centered fluxes in ( 6) and, re-integrating by parts yields:

         ( ∂H h ∂t , v) K i = -1 2 ((curlv, E h ) K i + (curlE h , v) K i ) + 1 2 l∈ν i E l × n, v F il , (ε ∞ ∂E h ∂t , v) K i = 1 2 ((curlv, H h ) K i + (curlH h , v) K i ) -1 2 l∈ν i H l × n, v F il , ( ∂J h ∂t , v) K i = -(γ d J h , v) K i + (ω 2 d E h , v) K i . (7)

Space-discrete formulation

We now define a set of scalar basis functions ϕ ik (1 ≤ k ≤ d i ) on the element K i , and then vectorize the basis functions (as well as the test functions), so in a 3D setting we have ∀i ∈ N Ω , 0 ≤ k ≤ d i , the vectorial basis functions

Φ 1 ik = (ϕ ik , 0, 0) T , Φ 2 ik = (0, ϕ ik , 0) T and Φ 3 ik = (0, 0, ϕ ik ) T .
The test function v can be chosen as the 3d i vectors Φ m ik (1 ≤ k ≤ d i , 1 ≤ m ≤ 3). Then the problem (7) can be denoted by an equivalent matrix form:

         M i dH i dt = -K i E i + l∈ν i S il E l , ε ∞ M i dE i dt = K i H i -M i J i - l∈ν i S il H l , dJ i dt = -γ d J i + ω 2 d E i . (8)

Space-discrete formulation

Here, the 3 × 3 block diagonal mass matrix M i , the 3 × 3 block diagonal stiffness matrix K i , and the 3 × 3 block diagonal surface matrix S il are defined by their respective diagonal blocks:

     ( M m i ) jk = (Φ m ij , Φ m ik ) K i , of size d i × d i , 1 ≤ m ≤ 3, ( K m i ) jk = 1 2 ((Φ m ij , curlΦ m ik ) K i + (curlΦ m ij , Φ m ik ) K i ), of size d i × d i , 1 ≤ m ≤ 3, ( S m il ) jk = Φ m ij , Φ m lk × n il F il , of size d i × d l , l ∈ ν i , 1 ≤ m ≤ 3,
here, d l denotes the number of DOFs on face l ∈ ν i .

Space-discrete formulation

Remark 1. For a metallic boundary face

F ij ∈ Γ m ∩ F B h , the metallic boundary condition is E h|F ij = 0. We therefore define E j|F ij = -E i|F ij , H j|F ij = H i|F ij . ( 9
)
For an absorbing boundary face

F ij ∈ Γ a ∩ F B h , we can define E j|F ij = z i N ij H i|F ij + E inc i|F ij -z i N ij H inc i|F ij , H j|F ij = -z -1 i N ij E i|F ij + H inc i|F ij + z -1 i N ij E inc i|F ij , (10)
where

N ij =   0 n ijz -n ijy -n ijz 0 n ijx n ijy -n ijx 0   , and E inc i|F ij , H inc i|F ij are given incident fields.

Space-discrete formulation

For all element, we can get a descriptor system of first order in time for (8) as following:

     M dH h dt = -KE h + S i E h + S e E h + B h u h (t), ε ∞ M dE dt = KH h -S i H h -S h H h -MJ h + B e u e (t), dJ h dt = -γ d J h + ω 2 d E h . (11)
The descriptor system of first order in time can then be written as:

M ẋ(t) = Kx(t) + bu(t), y (t) = p T x(t), (12

POD-DGTD formulations

Then considering the approximation of the state vector x(t) by another state vector, constrained to stay in the subspace spanned by the columns of U n , namely,

x(t) ≈ U n z(t), z(t) ∈ R n ,
yields an linear system with respect to the state variable z(t):

MU n ż(t) = KU n z(t) + bu(t), y (t) = p T U n z(t).

After left-multiplying the first equation by U n ,we have

U T n MU n ż(t) = U T n KU n z(t) + U T n bu(t), y (t) = p T U n z(t).
Then an n th reduced-order model of the linear system (12) in the time domain is naturally defined, and the quadruples (

M n , K n , b n , p T n ) can be simply expressed as (U T n MU n , U T n KU n , U T n b, p T U n ).

POD-DGTD formulations

Remark 2. For the discrete scheme ( 8), the number of DOFs (9N dof ) is very large, but the number DOFs n of the POD-DGTD formulations ( 22) is very small as we will show from the numerical results. Remark 3. The choice of the dimension d of the POD basis is not obvious. We denote by ε(ψ

1 , ψ 2 , • • • , ψ d ) =
r j=d+1 λ j the POD energy or error in the POD basis. If we want the POD energy to be less than some prescribed tolerance ρ, i.e., that ...

ε(ψ 1 , ψ 2 , • • • , ψ n ) ≤ ρ,

Padé approximation and moment-matching

Taking a Laplace transform for the descriptor system (12), assuming zero initial conditional, and then we obtain the following frequency domain formulation of ( 12):

sMX (s) = KX (s) + bU(s), Y (s) = p T X (s) ( 14
)
where X (s), Y (s), and U(s) respectively denotes the Laplace transform of x(x), y (x), and u(x). Eliminating X (s) in ( 14) leads to a transfer function of the original system ( 12) as

H(s) = Y (s) U(s) = p T (sM -K ) -1 b, ( 15
)
where s ∈ C, and H(s) ∈ R N-1,N is a scalar rational function.

Padé approximation and moment-matching

The Taylor series expansion of H(s) about

s i ∈ C (i = 1, 2, • • • , l) is
given by

H(s) = ∞ j=0 m j (s i )(s -s i ) j , ( 16
)
where m j (s i ) = p T [-(s i M -K ) -1 M] j (s i M -K ) -1 b are called moments about s i , l is the number of expansion points. We seek a reduced order system such that the Taylor expansion of the corresponding transfer function at s i has the form

H n (s) = ∞ j=0 m j (s i )(s -s i ) j , ( 17
)
where m j (s i ) are called moments about s i for a reduced order system.

Padé approximation and moment-matching

The function H n (s) ∈ R n-1,n is said to be an nth multi-point Padé approximation of H(s) about the expansion point

s i (i = 1, 2, • • • , l) if it
matches with the moments of H(s) as far as possible, i.e.

m j (s i ) = m j (s i ), j = 0, 1, 2, • • • , j i -1, i = 1, 2, • • • , l,
where j i > 0 denotes the number of the moments corresponding at s i , and n = l i=1 j i . Then,

H(s) = H n (s) + ∞ j=j i (m j (s i ) -m j (s i ))((s -s i ) j ), ( 18
)
where the last summation term denotes the error between the frequency response of the original system and that of the reduced-order system. If l = 1, ( 17) is a standard nth Padé approximant for H(s).

Padé approximation and moment-matching

If the transfer function H n (s) of the reduced-order system is expressed as a scalar rational function

H n (s) = P n-1 (s) Q n (s) , (19)
where the numerator P n-1 (s) and the denominator Q n (s) are defined respectively

P n-1 (s) = n-1 j=0 a j (s -s) j , Q n (s) = 1 + n j=1 b j (s -s) j ,
and the coefficients {a j } and {b j } are unknown, we can obtain a total of 2n linear equations by ( 18) and ( 19) by collecting the coefficients of the polynomials. Unfortunately, the coefficient matrix of the linear equations, which is called the Hanker matrix, usually is ill-conditioned.

Padé via Lanczos (l = 1, s0)

For the matrix B = -(K -s 0 M) -1 M, the vector r = (Ks 0 M)b, and a left staring vector p, the nth Krylov subspace K n (B, r ) is spanned by a sequence of n column vectors:

K n (B, r ) = span{r , Br , B 2 r , • • • , B n-1 r }.
This is sometimes called the right Krylov subspace. When the matrix B is nonsymmetric, there is a left Krylov subspace generated by B T and a starting vector p defined by

K n (B T , p) = span{p, B T p, (B T ) 2 p, • • • , (B T ) n-1 p}.
Note that the first 2n moments m i of H(s) are connected with Krylov subspaces through computing the inner products between the left and right Krylov sequences:

m 2i = ((B T ) i p) T • (B i b) T , m 2i+1 = ((B T ) i p) T • (B i+1 b) T , for i = 1, 2, • • • , n -1.

Padé via Lanczos

The left and right Krylov subspaces contain the desired information of moments, but the vectors B i r and (B T ) i p are unsuitable as basis vectors. The remedy is to construct more suitable basis vectors {v 1 , v 2 , • • • , v n } and {w 1 , w 2 , • • • , w n } such that they span the same desired Krylov subspaces. It is well-known that the Lanczos process is an elegant way to generate the desired basis vectors.

Lanczos process

Give a matrix B, a right staring vector r and a left staring vector p, the Lanczos process generates the basis vector {v i } and {w i }, known as the Lanczos vector. Moreover, these Lanczos vector are constructed to be biorthogonal (w j ) T v k = 0, for all j = k.

The Lanczos vectors can be generated by two three-term recurrences.

There recurrences can be stated compactly in matrix form as follows

BV n = V n T n + ρ n+1 v n+1 e T n , B T W n = V n T n + η n+1 w n+1 e T n ,
where T n and T n are the tridiagonal matrices.

Information about test problems:

The computational domain: the square Ω

2 = (x, y ) ∈ [0, l] × [0, l].
The boundary condition: PEC.

The Drude parameters: ε ∞ = 3.7362, ω d = 1.3871e7 GHz, γ d = 4.41544e4 GHz.

Codes: MATLAB 2013b (Memory: 8Gb).

Method: DG methods with RK 4 , PVL-DGTD methods. It can be seen that the convergence order of numerical simulation is to satisfy the theoretical convergence order. the difference between the error calculated by the DGTD method and that calculated by the PVL-DGTD is small. while the dimension of the reduced order model is much smaller than that of the original system. It is seen that the PVL-DGTD method saves the required CPU time to yield a uniform accuracy level.

The point with coordinates: (4.560e-8, 6.625e-7); Methods: DGTD-P 2 with M2, PVL-DGTD-P 2 with M2. We again find that PVL-DGTD solutions coincide with the DGTD solutions and the analytical solutions.

  DG formulations POD-DGTD formulations Krylov subspace based solver MOR Proper orthogonal decomposition (POD) methods [Holmes P, Lumley J L, Berkooz G, 1998]. Krylov subspaces methods. Padé-via-Lanczos (PVL) [R Freund, P Feldmann, 1993]. the Arnoldi method. the multi-point rational Krylov subspaces [K Gallivan, F ij ∈ F I h , we define average values {•}:

3

  Reduced-order DG formulations POD-DGTD formulations Krylov subspace based solver Krylov-subspace based operator exponential approach Krylov-subspace based Padé approximation approach 4 Numerical experiments Artificial validation case Scattering of a plane wave by a Au-nanowire

  ) with initial condition x(0) = x 0 = [E 0 (x), H 0 (x), J 0 (x)] T , where M, K ∈ R N×N , p, b ∈ R N , and N = 9N dof . DG formulations POD-DGTD formulations Krylov subspace based solver Krylov-subspace based operator exponential approach Krylov-subspace based Padé approximation approach Numerical experiments Artificial validation case Scattering of a plane wave by a Au-nanowire Summary

  we can choose d to be the minimum of positive integer such that πScattering of a plane wave by a Au-nanowire Summary Krylov-subspace based operator exponential approach Introducing the system matrix A asA = -M -1 K ,stable field approximations can be computed as [Zimmerling J, Wei L, Urbach P, et al, JCP, 2016]x(t) = u(t) * 2η(t)Re[η(A)exp(-At)M -1 b] for t > 0,where η(t) is the Heaviside unit step function. The central approximation of the Krylov subspace techniques described then reads asexp(-At)M -1 b = V n exp(-H n t)e 1 ,e 1 being a first unit basis vector, n being the dimensional of the reduced order system, and V n , H n can be obtained by Arnoldi-process or rational Lanczos-process. The references for Krylov-subspace based operator exponential approach: Zimmerling J, Wei L, Urbach P, et al. A Lanczos model-order reduction technique to efficiently simulate electromagnetic wave propagation in dispersive media. Journal of Computational Physics, 2016, 315: 348-362. Botchev M A. Krylov subspace exponential time domain solution of Maxwells equations in photonic crystal modeling. Journal of computational and applied mathematics, 2016, 293: 20-34. Druskin V, Lieberman C, Zaslavsky M. On adaptive choice of shifts in rational Krylov subspace reduction of evolutionary problems. SIAM Journal on Scientific Computing, 2010, 32(5): 2485-2496.

Figure : Figure :

 :: Figure: Comparisons between the usual DGTD methods and the PVL-DGTD methods. (a)∼(c): the usual DGTD solutions for Ez , Hx , Hy . (d)∼(f): the PVL-DGTD solutions for Ez , Hx , Hy .

Field

  distributions: E x . Time: t = (6*(1e-16)).

  Figure: 1D distribution along the x for the modulus of Ex in the Fourier domain.

Table :

 : Convergence results for the dispersive cavity problem.

			DGTD-P 1		DGTD-P 2		DGTD-P 3		DGTD-P 4
	h	L 2 -error	r	L 2 -error	r	L 2 -error	r	L 2 -error	r
	3.678e-08	3.482e-08	-	9.364e-10	-	4.161e-11	-	6.735e-13	-
	1.839e-08	1.618e-08	1.106	1.965e-10	2.253	5.140e-12	3.017	4.091e-14	4.042
	9.194e-09	7.931e-09	1.028	4.678e-11	2.070	6.448e-13	2.995	2.585e-15	3.984
	4.597e-09	3.946e-09	1.007	1.161e-11	2.011	8.099e-14	2.993	1.644e-16	3.975

Table :

 : Computational performance of the artificial problem.

		DGTD		PVL-DGTD
	Mesh	N	L 2 -error	n	s0	L 2 -error
			P1
	M1	3744 3.482e-08	25	1300f 3.489e-08
	M2	14976 1.618e-08	25	2100f 1.618e-08
	M3	59904 7.931e-09	35	3500f 7.930e-09
			P2
	M1	7488 9.364e-10	35	1700f 9.366e-10
	M2	29952 1.965e-10	35	3700f 1.964e-10
	M3	119808 4.678e-11	90 12000f 4.684e-11
			P3
	M1	12480 4.161e-11	45	1500f 4.162e-11
	M2	49920 5.140e-12	75	3100f 5.799e-12
	M3	199680 6.448e-13	100 15000f 6.518e-13
			P4
	M1	18720 6.735e-13	52	2100f 6.907e-13
	M2	74880 4.092e-14	90 10000f 4.292e-14
	M3	299520 2.585e-15	150 20000f 3.473e-15

Table :

 : Computational performance of the artificial problem.

		DGTD			PVL-DGTD
	Mesh Tconstruction(s) T solution (s)	Tconstruction(s) TPVL(s) T solution (s) Speed-up
			P1
	M1	0.566	0.405	0.566	0.148	0.026 57.037%
	M2	2.535	3.135	2.535	0.990	0.091 65.518%
	M3	14.499	29.041	14.499	9.137	1.550 63.200%
			P2
	M1	0.601	3.109	0.601	0.613	0.059 78.385%
	M2	2.873	28.020	2.873	3.725	0.437 85.146%
	M3	22.994	285.670	22.994	58.094	8.376 76.732%
			P3
	M1	0.638	10.680	0.638	2.623	0.147 74.064%
	M2	3.765	115.400	3.765	23.464	2.232 77.733%
	M3	52.804	1135.300	52.504	224.765	28.818 77.664%
			P4
	M1	0.896	36.593	0.896	7.689	0.360 78.004%
	M2	7.155	334.265	7.155	92.095	4.832 71.003%
	M3	126.198	3148.400	126.198	783.482	74.960 72.734%

Space-discrete formulation

POD-DGTD formulations

be a collection of snapshots of the solution of ( 12). We compute the singular value decomposition and truncate depending on how fast the singular values decay:

Here, r = rank(W ), U and V are N × N and × unitary matrices. Let

Krylov subspace based solver

Reduced Order Model

In general, the application of model order reduction to (12) yields a reduced order model by using the orthogonal projection

where M n , K n ∈ R n×n , and b n , p n ∈ R n have a significantly smaller dimension n N than the original system (12), and

Why is the simple choice l = 1, s 1 = s 0 ?

), e.g, the rational Arnoldi method, the selection of suitable expansion points (also called shifts) is crucial for the quality of the resulting reduced order model.

If l = 1, s 1 = infinit, e.g, PVL with expansion at infinity is fast while the accuracy of the method is very low. The computational domain: a gold disk of radius r = 2nm.

The boundary condition: ABC, and propagating along the -y axis.

The incident plane wave (a sinusoidal modulated Gaussian polse):

The Drude parameters:

Conclusions and Future works

Conclusions:

Conclusions and Future works

Conclusions:

The reduced-order DG method is developed for saving time;

The reduced-order DG method has less degrees of freedom compared to the usual DG method;

Numerical results show their advantages.

Conclusions and Future works

Conclusions:

The reduced-order DG method is developed for saving time;

The reduced-order DG method has less degrees of freedom compared to the usual DG method;

Numerical results show their advantages.

Future works:

Conclusions and Future works

Conclusions:

The reduced-order DG method is developed for saving time;

The reduced-order DG method has less degrees of freedom compared to the usual DG method;

Numerical results show their advantages.

Future works:

More complex 3d realistic applications should be investigated;

The error bounds between the usual DG solutions and the reducedorder DG solutions will be derived;

The extension of this reduce-order method to more general physical model, e.g., Hydrodynamic Drude and generalized non-local modes.

other MOR methods, such as, rational Krylov subspace technique.