Flávio Luiz Cardoso-Ribeiro
email: flaviocr@ita.br

Denis Matignon

Laurent Lefèvre

A partitioned finite element method for power-preserving discretization of open systems of conservation laws

Keywords: geometric spatial discretization, structure-preserving discretization, port-Hamiltonian systems, partitioned finite element method

This paper presents a structure-preserving spatial discretization method for distributed parameter port-Hamiltonian systems. The class of considered systems are hyperbolic systems of two conservation laws in arbitrary spatial dimension and geometries. For these systems, a partitioned finite element method (PFEM) is derived, based on the integration by parts of one of the two conservation laws written in weak form. The non-linear one-dimensional shallow-water equation (SWE) is first considered as a motivation example. Then, the method is investigated on the example of the non-linear two-dimensional SWE. Complete derivation of the reduced finite-dimensional port-Hamiltonian system (pHs) is provided and numerical experiments are performed. Extensions to curvilinear (polar) coordinate systems, spacevarying coefficients and higher-order pHs (Euler-Bernoulli beam equation) are provided.

Introduction

The port-Hamiltonian formalism has proven to be a powerful tool for the modelling and control of complex multiphysics systems. In many cases, spatio-temporal dynamics must be considered and infinite-dimensional port-Hamiltonian models are needed [START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF]. Classical academic examples such as the transmission line, the shallow water or the beam equations have been investigated in the port-Hamiltonian framework [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF].

Besides, two-and three-dimensional problems have been recently considered [START_REF] Trenchant | Structure preserving spatial discretization of 2D hyperbolic systems using staggered grids finite difference[END_REF][START_REF] Vu | A structured control model for the thermo-magnetohydrodynamics of plasmas in tokamaks[END_REF][START_REF] Wu | Power preserving model reduction of 2D vibroacoustic system: a port Hamiltonian approach[END_REF]. In many of these examples, e.g. those arising from fluid mechanics, systems of two balance equations are considered such as mass and momentum or volume and momentum balance equations.

In order to simulate and design control laws, obtaining a finite-dimensional approximation which preserves the port-Hamiltonian structure of the original system can be advantageous. It may serve as a design guide such as in control by interconnection or in interconnection and damping assignment passivity-based control (IDA-PBC). Besides, preserving the underlying Dirac interconnection structure results in energy conservation properties and associated dynamical properties (e.g. stability, controllability, etc.).

Mixed finite element methods were introduced a long time ago to perform structure-preserving spatial discretization of the Maxwell field equations [START_REF] Bossavit | Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism[END_REF][START_REF] Bossavit | Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements[END_REF]. [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF] applied these mixed finite element methods to open systems with boundary energy flows. More precisely, they considered a mixed finite element structure-preserving spatial discretization for one-dimensional hyperbolic systems of conservation laws, making use of distinct low-order Whitney basis functions to approximate the energy and co-energy variables, respectively. This idea was applied later for the discretization of a parabolic diffusion problem related to pressure-swing-adsorption columns (Baaiu et al., 2009b) and piezo-electric beams [START_REF] Cardoso-Ribeiro | Piezoelectric beam with distributed control ports: a power-preserving discretization using weak formulation[END_REF] and was connected to finite volume and staggered grid finite difference methods for one-dimensional problems [START_REF] Kotyczka | Finite volume structure-preserving discretization of 1D distributed-parameter port-Hamiltonian systems[END_REF]. They were also generalized for two-dimensional systems by [START_REF] Wu | Power preserving model reduction of 2D vibroacoustic system: a port Hamiltonian approach[END_REF] where they are applied to a vibro-accoustic systems. [START_REF] Trenchant | Structure preserving spatial discretization of 2D hyperbolic systems using staggered grids finite difference[END_REF] considered two-dimensional finite difference staggered grids schemes for the same vibro-accoustic system. Finally, these structure-preserving mixed finite element methods, applied to the spatial discretization of general port-Hamiltonian systems (pHs) with boundary energy flows, were stated in a geometry independent form making use of discrete exterior calculus results by [START_REF] Seslija | Explicit simplicial discretization of distributedparameter port-Hamiltonian systems[END_REF].

In these previous works, the central idea was to define different discretization bases for the energy and co-energy variables such that the strong form of the equations was exactly satisfied in the corresponding spanned finite-dimensional approximation spaces. This idea was extended to geometric pseudo-spectral methods using conjugate high-order polynomial bases [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF] or Bessel [START_REF] Vu | Symplectic spatial integration schemes for systems of balance equations[END_REF] basis functions, globally defined on the whole spatial domain. Farle and his co-authors also used different approximation bases for the discretization of the one-dimensional transmission line and three-dimensional Maxwell's equations (Farle et al., 2014a(Farle et al., ,b, 2013)). In these latter works, one of the balance laws is kept in strong form (with exact spatial derivation) while the other one is considered in the weak sense only and is being integrated by parts. A metric dependent Stokes-Dirac structure is introduced, making use of the Hodge star duality product, to make the element spaces compatible. As it was noticed by [START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF], defining these compatible spaces-with power-conjugate approximation bases for the energy and co-energy variables-is (relatively) straightforward for one-dimensional systems but seems to be cumbersome for higher spatial dimensions or higher-order methods. The discretization of the co-boundary (exterior derivative) and boundary operators may be performed independently such that the Stokes theorem still holds in the chosen approximation dual spaces. Then, for complex geometries (e.g. non-rectangular domain in N-dimensional spaces or non-smooth boundary, etc.) the approximated relations between indomain and boundary conjugate variables may become intricate, and the expression of the kernel rather involved, potentially leading to dimensionality problems. As suggested by [START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF] the discretization of the weak formulation of the considered pHs may be a practical solution to deal with these higher dimensional problems or more complex geometries. We propose in this paper to follow this approach but to perform integration by parts-which was used by [START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF] to get the weak formulation-on one of the two balance equations only, defining in this way a partitioned mixed finite element method. Doing so, the discretization in the chosen bases for the energy and co-energy variables (and the associated test functions) directly leads to a finite-dimensional Dirac interconnection structure and no further projection is required to get finite-dimensional pHs equations with reversible causality. Besides, boundary conditions are naturally handled, even in the case of higher-order finite element bases. Finally, the use of this weak-form formulation enables the use of standard finite-element software to perform the proposed structure-preserving discretization and conse-quently paves the way for further applications to more involved higher-order problems with complex geometries.

This paper starts with a motivation example detailed in Section 2: a structure-preserving spatial discretization for the one-dimensional shallow-water equations (1D SWEs). In Section 3, the approach is generalized to higher spatial dimensions, where the initial model is stated in vector calculus, independent of the specific geometry. The formulation of systems of two conservation laws with boundary energy flows is presented, as well as their weak form, making use of Cartesian coordinates. In Section 4, the proposed partitioned finite element method (PFEM) is applied to the two-dimensional shallow-water equation (2D SWE) example. Numerical experiments are presented in Section 5. Finally, extensions of the method to curvilinear (polar) coordinate systems, space-varying coefficients and higher-order pHs (Euler-Bernouilli beam equations) are provided in Section 6. The paper ends with conclusions and open questions which are discussed in Section 7.

An introductory example

The aim of this section is to present the general idea of this paper-partial integration by parts of the weak form for systems of two conservation laws and structure-preserving projections in finitedimensional approximation spaces-applied on a simple one-dimensional example, namely the 1D SWE written in the port-Hamiltonian formulation. First, the port-Hamiltonian formulation for these equations is recalled (Subsection 2.1). Then, the partial integration by parts idea is performed on the weak form for this 1D SWE example and a structure-preserving finite element spatial discretization method is applied to obtain the finite-dimensional pHs (Subsection 2.2). This general idea differs from previous works [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF][START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF] where the central idea was to define different discretization spaces for the energy and co-energy variables such that the strong form of the equations were exactly satisfied in these finite-dimensional spaces. Instead, we use a weak-form representation for the equations, where only one of the conservation laws is integrated by parts. This partitioned approach naturally leads to a skew-symmetric interconnection matrix between the flow and effort variables. Furthermore, the use of weak form enables the use of classical finite-element methods to perform the discretization.

Port-Hamiltonian strong formulation for the 1D SWE

The SWEs are sets of partial differential equations that can be used to represent an incompressible fluid with free-surface motion. These equations are typically used to model fluid motion in water channels [START_REF] Hamroun | Control by interconnection and energy-shaping methods of port Hamiltonian models. Application to the shallow water equations[END_REF], wave propagations in oceans and lakes and sloshing in fluid tanks [START_REF] Alemi Ardakani | A symplectic integrator for dynamic coupling between nonlinear vessel motion with variable cross-section and bottom topography and interior shallow-water sloshing[END_REF][START_REF] Cardoso-Ribeiro | A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system[END_REF], 2020c). When one considers the frictionless flow in a horizontal channel with uniform rectangular cross-section, the one-dimensional mass and momentum balance equations may be written as

∂ ∂t h = - ∂ ∂z (hu) , ∂ ∂t u = - ∂ ∂z u 2 2 + gh , (2.1)
where h(z, t) is the fluid height, u(z, t) is the fluid average velocity in a cross-section, z is the spatial coordinate, t is the time and g is the gravitational acceleration.

The total energy of the system inside the domain [0, L] is given by the sum of kinetic and potential (gravitational) energy:

H = 1 2 [0,L] ρb hu 2 + ρbg h 2 dz , (2.2)
where b is the width of the water channel (or tank) rectangular cross-section and ρ the water density. Defining the energy-variables q(z, t) := bh(z, t) and p(z, t) := ρu(z, t), the system Hamiltonian (total energy) is given by

H (q(z, t), p(z, t)) = 1 2 [0,L] qp 2 ρ + ρg b q 2 dz. (2.3)
Using these newly defined variables, (2.1) can be rewritten as

q(z, t) = - ∂ ∂z e p (z, t), ṗ(z, t) = - ∂ ∂z e q (z, t), (2.4)
where e q (z, t) and e p (z, t) are the co-energy variables (respectively, the total pressure and the water flow) which are defined as the variational derivatives 1 of the Hamiltonian with respect to q(z, t) and p(z, t):

e q = δH δq = p 2 2ρ + ρg b q = ρ u 2 2 + gh , e p = δH δp = qp ρ = bhu.
(2.5)

From the above definitions of energy and co-energy variables, using the SWE written in the canonical Hamiltonian form (2.4) and Stokes theorem, one obtains for the power balance equation

Ḣ(t) = [0,L] e q (z, t)q(z, t) + e p (z, t)ṗ(z, t) dz, = - [0,L] ∂ ∂z e q (z, t)e p (z, t) dz, = - ∂[0,L]
e q (z, t)e p (z, t),

= u T ∂ y ∂ , (2.6)
1 For the definition of variational derivative, see, e.g. (Olver, 1993, Definition 4.1) or (Duindam et al., 2009, Chapter 4).

where boundary port input variables, u ∂ , are defined as the values of the co-energy variables evaluated at the spatial domain boundary The time derivatives of the energy variables are usually called flow variables. We define f q (z, t) := -q(z, t) and f p (z, t) := -ṗ(z, t). The co-energy variables are also called effort variables. The power balance equation (2.6) defines a natural pairing or bilinear form

• |• : B → R , (e, f ) → e |f := [0,L] e q (z, t)f q (z, t) + e p (z, t)f p (z, t) dz + u T ∂ y ∂ , (2.9)
where the bond space B := E × F is defined as the product of the effort real vector space E := e := e q e p e ∂ T e q , e p ∈ H 1 (0, L) ; e ∂ ∈ R 2 , (2.10) and its power-conjugate flow vector space (2.11) with H 1 (0, L) and L 2 (0, L) denoting the Sobolev space of functions with square integrable derivatives on [0, L] and the usual Lebesgue space of square integrable functions on [0, L], respectively. Using the bilinear form (2.9), the power balance equation (2.6) simply reads

F := f := f q f p f ∂ T f q , f p ∈ L 2 (0, L) ; f ∂ ∈ R 2 ,
⎡ ⎣ e q (t, •) e p (t, •) e ∂ (t) ⎤ ⎦ ⎡ ⎣ -q(t, •) -ṗ(t, •) f ∂ (t) ⎤ ⎦ = 0, (2.12)
with e T ∂ (t) = e q (0, t) -e q (L, T) and f T ∂ (t) = e p (0, t) e p (L, T) (or the reverse). Besides, the pairing (2.9) may be symmetrized to obtain the associated indefinite bilinear form

• |• : B × B → R (e 1 , f 1 ), (e 2 , f 2 ) → (e 1 , f 1 ), (e 2 , f 2 ) := 1 2 e 1 f 2 + e 2 f 1 .
(2.13)

It may be shown that the Hamiltonian formulation (2.4) for the SWE, together with the boundary conditions (2.7) and output (2.8), may be equivalently implicitly defined as

δH δq , δH δp , u ∂ , - dq dt , - dp dt , y ∂ ∈ D , (2.14)
where D ⊂ B is the linear subspace which is maximally isotropic (i.e. D = D ⊥ ) with respect to the inner product (2.13) [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF]. In that sense, the natural pairing (2.9) fully describes the geometric structure of the pHs (2.4) with boundary values (2.7). Therefore, in this paper, structurepreserving (or symplectic) spatial discretization will be understood as approximations (projections) which preserve this power form (2.9). Symplecticity in that sense implies not only preservation of the power balance (2.6) or (2.12) (i.e. isotropy) but also preservation of the whole geometric structure of the system (e.g. the Poisson structure in the example of closed systems or the Dirac structure in the case of open systems with time-varying boundary conditions) [START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF][START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF]. Note that the particular input and output port variables chosen here above in (equations (2.7) and (2.8)) is only one among other possible choices. A description of all the possible choices of input/output variables which lead to well-posed problems (in the linear case) is described by Le [START_REF] Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF].

Partitioned weak-form and structure-preserving discretization for the 1D SWE

We will now introduce a weak formulation for the 1D SWE and then perform integration by parts on the first balance equation. Let v q (z) ∈ H 1 (0, L) and v p (z) ∈ L 2 (0, L) denote any arbitrary test functions, we may obtain from the strong formulation (2.4) the following weak form:

[0,L] v q (z)q(z, t)dz = - [0,L] v q (z) ∂ ∂z e p (z, t)dz, [0,L] v p (z)ṗ(z, t)dz = - [0,L] v p (z) ∂ ∂z e q (z, t)dz.
(2.15)

Integrating by parts the first balance equation only, we get the following partitioned weak form:

[0,L] v q (z)q(z)dz = [0,L] e p (z, t) ∂ ∂z v q (z)dz -v q (L)e p (L, t) + v q (0)e p (0, t), [0,L] v p (z)ṗ(z)dz = - [0,L] v p (z) ∂ ∂z e q (z, t)dz.
(2.16)

Remark 2.2 In the specific case where v q (z) = 1 and v p (z) = 1, we get [0,L] q(t, z)dz = e p (0, t)e p (L, t), [0,L] ṗ(t, z)dz = e q (0, t)e q (L, t),

(2.17) which show that the two conservation laws for the total mass and the total momentum in the spatial domain [0, L] are preserved in the weak formulation. When v q = e q (z, t) and v p = e p (z, t) are chosen, one gets

[0,L] e q (z, t)q(z)dz = [0,L]
e p (z, t) ∂ ∂z e q (z)dze q (L, t)e p (L, t) + e q (0, t)e p (0, t),

[0,L] e p (z, t)ṗ(z)dz = - [0,L]
e p (z) ∂ ∂z e q (z, t)dz.

(2.18)

Therefore, the power balance equation (2.6) reads

Ḣ = [0,L]
(e q (z, t)q(z) + e p (z, t)ṗ(z))dz = -e q (L, t)e p (L, t) + e q (0, t)e p (0, t), (2.19) which shows that the power balance is also preserved in the weak formulation.

We will now project the partitioned weak formulation (2.16) into finite-dimensional approximation spaces chosen in such a way as to preserve the total mass and momentum conservation laws, the power balance equation and the underlying Dirac structure of the original port-Hamiltonian model (2.4). Unlike in [START_REF] Moulla | Pseudo-spectral methods for the spatial symplectic reduction of open systems of conservation laws[END_REF] and [START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF] where different approximation bases are chosen for the energy and co-energy variables, we obtain the mass, momentum, power and structure-preservation by the selection of different approximation bases for the mass and momentum densities. This 'partitioned' choice for the approximation bases lead us to square skew-symmetric interconnection matrices.

Let us approximate the energy variables q(z, t) and p(z, t) as

q(z, t) ≈ q ap (z, t) := N q i=1 φ i q (z)q i (t) = φ T q (z)q(t), p(z, t) ≈ p ap (z, t) := N p j=1 φ j p (z)p j (t) = φ T p (z)p(t),
(2.20)

where φ i q (z) , i ∈ {1, . . . , N q } are the chosen approximation basis functions in H 1 (0, L), φ i p (z) , i ∈ {1, . . . , N p } the chosen approximation basis functions in L 2 (0, L), while q(t) and p(t) are the approximation coordinates for q ap (z, t) and p ap (z, t) in the approximation bases φ q (z) and φ p (z). The test functions v q (z) and v p (z) are approximated in the same bases as q(z, t) and p(z, t), respectively. From the substitution of the approximated variables (2.20) in the weak form (2.16), the following finitedimensional equations are obtained:

v q T [0,L] φ q (z)φ T q (z)dz q(t) = v q T [0,L]
dφ q dz (z)φ T p (z)dz e p (t)+ v T q φ q (0)e ap p (0, t)-v T q φ q (L)e ap p (L, t),

v p T [0,L] φ p (z)φ T p (z)dz ṗ(t) = -v p T [0,L] φ p (z) dφ T q dz (z)dz e q (t),
(2.21)

where the H 1 (0, L) effort functions e q (•, t) have been approximated in the {φ i q (z)} basis, while the L 2 (0, L) effort functions e p (•, t) have been approximated in the {φ i p (z)} basis. Similarly, the flow functions q(•, t) and ṗ(•, t) have been approximated in the {φ i q (z)} and {φ i p (z)} bases, respectively. Since these equations should remain valid for any choices of test functions coordinates v 1 and v 2 , one gets

M q q(t) = De p (t) + B e p (0, t) e p (L, t) , M p ṗ(t) = -D T e q (t), (2.22)
where M q and M p are square mass matrices (of size N q × N q and N p × N p , respectively) defined as

M q := [0,L] φ q (z)φ T q (z)dz , M p := [0,L] φ p (z)φ T p (z)dz. (2.23)
Matrix D is of size N q × N p and is defined as (2.24) and B := φ q (0) -φ q (L) is an N q × 2 matrix. Using the input-output conjugate boundary port variables as defined in (2.7) and (2.8), since the boundary values of e ap q (z, t) may be written as (2.25) the approximation (2.22) may be written using the following finite-dimensional Dirac structure representation:

D := [0,L] dφ q dz (z)φ T p (z)dz,
e ap q (0, t) -e ap q (L, t) = φ T q (0) -φ T q (L) e q = B T e q ,
M q 0 0 M p f q (t) f p (t) = 0 D -D T 0 e q (t) e p (t) + B 0 e p (0, t) e p (L, t) , e ap q (0, t) -e ap q (L, t) = B T 0 e q (t) e p (t) , (2.26)
where f q (t) and f p (t) denote the vector coordinates for the flow approximations in the φ i q (z) and

φ j p (z) approximation bases, that is, f ap q (z, t) = φ T q (z)f q (t) and f ap p (z, t) = φ T p (z)f p (t)
, respectively. We obtained the finite-dimensional Dirac structure representation (2.26) from the projection of the Stokes-Dirac structure on the chosen approximation spaces. We will now derive the corresponding approximation of the pHs dynamics (2.14) by restricting the Hamiltonian functional to the same approximation spaces. From the definition of co-energy variables as variational derivatives of the Hamiltonian with respect to q and p, we get

Ḣ(t) = [0,L]
e q (z, t)q(z, t) + e p (z, t)ṗ(z, t)) dz.

(2.27)

Using the approximations (2.20) for the energy and co-energy variables, this power balance may be approximated as

Ḣd (t) := [0,L]
e ap q (z, t)q ap (z, t) + e ap p (z, t)ṗ ap (z, t) dz, = e T q (t)M q q(t) + e T p (t)M p ṗ(t).

(2.28)

Therefore, in order to write the power balance (2.28) as the total time derivative of the discrete Hamiltonian written as a function of the finite-dimensional vector coordinates for the energy variables, the following relationship between the co-energy variables and the discrete Hamiltonian gradient must hold:

M q e q (t) = ∂H d ∂q , M p e p (t) = ∂H d ∂p , (2.29)
where

H d (q, p) := H φ T q (z)q(t), φ T p (z)p(t) .
(2.30)

We obtain the finite-dimensional pHs formulation for the proposed structure-preserving reduction scheme by combining equations (2.26) (for the linear finite-dimensional Dirac interconnection structure) and the non-linear constitutive equations (2.29) and (2.30).

In the 1D SWE example, the Hamiltonian function is neither quadratic nor separable. Nevertheless, an explicit form may be obtained for the constitutive equations (2.29). Since the Hamiltonian function restricted to the approximation spaces for q and p reads

H d (q, p) := 1 2 [0,L] ⎛ ⎜ ⎝ φ T q (z)q φ T p (z)p 2 ρ + ρg b φ T q (z)q 2 ⎞ ⎟ ⎠ dz.
(2.31)

One obtains for the effort variables the expressions

M q e q (t) = ∂H d ∂q = ρg b M q q(t) + 1 2ρ [0,L] φ q (z)p T (t)φ p (z)φ T p (z)dz p(t), M p e p (t) = ∂H d ∂p = q T (t) 1 ρ [0,L] φ q (z)φ p (z)φ T p (z)dz p(t).
(2.32)

Note that both constitutive equations exhibit non-linear terms. In order to compute them, the following procedure was used. The first equation can be written as

∂H d ∂q = ρg b M q q(t) + 1 2ρ [0,L] φ q (z)p T (t)φ p (z)φ T p (z)dz p(t) w(t) , (2.33)
where the components of w(t) can be computed as

w i (t) = p T (t) [0,L] φ q,i (z)φ p (z)φ T p (z)dz p(t).
(2.34)

Note that [0,L] φ q,i (z)φ p (z)φ T p (z)dz, for {i = 1, . . . , N q } are N q matrices of dimension N p × N p , which can be computed once and remain constant. Similarly, the second constitutive relationship can be written as a function of constant matrices.

Remark 2.3 Matrices M q := [0,L] φ q φ T q dz and M p := [0,L] φ p φ T p dz are the classical mass matrices of the FEM. Since the basis functions have very small support (the stencil), often non-overlapping, these matrices are diagonal for P0-Lagrange elements, tridiagonal for P1-Lagrange elements, pentadiagonal for P2-Lagrange elements or indeed sparse, even if an integral over the whole domain [0, L] appears. Similarly, in (2.34), the matrices obtained are even sparser, since the support of three different basis functions are involved in the computation of the integral. See the illustration on Fig. 6.

Remark 2.4 We may deduce from the pHs representation (2.41), the non-linear constitutive equations (2.46) and (2.48) that the power balance equation reads Ḣd = -e 1 (0, t)e 2 (0, t) + e 1 (L, t)e 2 (L, t).

(2.35)

Hence, the power balance (and the corresponding power product value) is preserved by the proposed partitioned spatial discretization scheme. In that sense, we call it a structure-preserving or symplectic scheme. When the basis functions satisfy

N q i=1 φ i q (z) = N p i=1 φ i p (z) = 1 , ∀z ∈ [0, L],
(2.36)

the mass and momentum conservation laws are also satisfied in the finite-dimensional approximation spaces:

[0,L] φ q (z) T dz q(t) = [0,L] qap (z, t)dz = e p (0, t) -e p (L, t), [0,L] φ p (z) T dz ṗ(t) = [0,L]
ṗap (z, t)dz = e q (0, t)e q (L, t).

(2.37)

In the beginning of this section, we motivated this work by the fact that previous work on structure-preserving spatial discretization that relies on exact satisfaction of the strong form of the equations usually lead to difficulties when generalizing to two-or three-dimensional systems. The following questions arise: can the proposed PFEM method be easily generalized to higher-dimensional problems (two-and three-dimensional)? Does it work with different coordinate systems? What about convergence? We will answer these questions in the following sections.

A general setting

In this section, we generalize the definition of pHs, given in the previous section for the 1D SWE example, for systems of two conservation laws in arbitrary dimension involving the (divergence, gradient) system in vectorial representation. This will allow us to define, independently from the particular spatial dimension, geometry or coordinate system, the class of problems which can be solved by using the structure-preserving spatial discretization scheme proposed in this paper.

In Section 3.1, the class of pHs in dimension N is defined. In Section 3.2, we formulate the 2D SWE example in this setting. In Section 3.3, we define the Stokes-Dirac interconnection structures associated to these pHs and the corresponding power pairings which will be preserved in the discretization. Finally, in Section 3.4, we will give the general weak form which will be used for the structure-preserving discretization scheme presented next in Section 4.

Port-Hamiltonian formulation of systems of two conservation laws with boundary energy flows: a vector calculus statement in higher space dimension

We will now extend the port-Hamiltonian formulation which has been presented in Section 2.1 only for the 1D SWE example. Hyperbolic systems of two conservation laws will be stated using vector calculus in R n . Let us consider the two conserved quantities, a scalar one α q and a vector-valued one α p ; the geometric domain of interest is a connected domain Ω ⊂ R n with Lipschitz boundary ∂Ω. These variables α q = α q (z, t) and α p = α p (z, t) are scalar-or vector-valued distributed energy state variables defined for any z ∈ Ω (z is the position vector) and time t 0. Let the Hamiltonian functional H be defined as

H(α q , α p ) := Ω H α q (z, t), α p (z, t), z d z, (3.1)
where H denotes the Hamiltonian density which is assumed to be a smooth function. The variational derivatives of H with respect to α q and α p are the unique scalar δ q H and vector-valued δ p H functions2 , such that

H(α q + εδα q , α p + εδα p ) = H(α q , α p ) + ε Ω δ q H δα q + δ p H • δα p d z + o(ε). (3.2)
Therefore, from the Hamiltonian defined in (3.1), we may define the co-energy variables (efforts):

e q := δ q H, e p := δ p H. (3.3)
In practice, when the canonical scalar product is being used, and when the functional H does not depend on the derivatives of the functions α q and α p , then the variational derivative amounts to the partial derivative δ q H = ∂ α q H and the gradient δ p H = grad α p H , respectively.

The Hamiltonian system of two canonically interacting conservation laws for α q and α p may be defined as

αq (z, t) αp (z, t) = 0 -div -grad 0 J e q (z, t) e p (z, t) , (3.4)
and since the formal adjoint of the divergence is the gradient, then the matrix-valued differential operator J remains formally skew-symmetric (that is, skew-symmetry assuming zero boundary conditions for the arguments). According to (3.2), the time derivative of the energy functional (power balance) reads

Ḣ = Ω δ q H αq + δ p H • αp d z = δ q H, αq L 2 (Ω) + δ p H, αp L 2 (Ω;R n ) . (3.5)
According the state equations (3.4) and to the celebrated Green's formula, this power balance may be written as

Ḣ = -e q , div e p L 2 (Ω) -e p , grade q L 2 (Ω;R n ) = - Ω div(e q e p ) d z = - ∂Ω e q ∂Ω e p ∂Ω • n ds. (3.6)
This latter formula suggests the following boundary port variables:

e ∂ f ∂ = γ 0 0 0 -γ ⊥ δ q H(z, t) δ p H(z, t) , (3.7)
involving the Dirichlet trace operator γ 0 : H 1 (Ω) → H 1 2 (∂Ω), defined by γ 0 (e q ) = e q ∂Ω , and the

normal trace operator γ ⊥ : H div (Ω) → H -1 2 (∂Ω), defined by γ ⊥ (e p ) = e p ∂Ω
• n; we refer the reader to (Tucsnak & Weiss, 2009, Chapter 10) for a precise definition of the functional spaces and both the boundary operators. These boundary port variables are defined in such a way that the power balance equation (3.5) and (3.6) may be written as

e q , f q L 2 (Ω) + e p , f p L 2 (Ω;R n ) + e ∂ , f ∂ ∂Ω = 0. (3.8)
The last term on the boundary is not a scalar product in general but stands for a duality bracket between H 1 2 (∂Ω) and H -1 2 (∂Ω); moreover, the flow variables f q and f p are defined as

f q (z, t) := -αq (z, t), f p (z, t) := -αp (z, t).
(3.9)

The power balance equation (3.8) simply states that the time derivative of the energy increase inside the domain Ω equals the power supplied through the boundary ∂Ω. As it may be seen, we have extended the systems of two conservation laws (3.4) with a boundary power supply and related boundary port variables, obtaining an open system of two conservation laws with boundary energy flows. The explicit definition of an open pHs of two canonically interacting conservation laws is given by the distributed state equation (3.4) together with the definition of the boundary port variables (3.7). It leads to the structural power balance equation (3.8) which is independent of the specific considered Hamiltonian function (i.e. from the effort constitutive equations (3.2)). Many one-, two-or threedimensional examples, either linear or non-linear, may be recast in that framework and satisfy this definition [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: an introductory overview[END_REF]. Even parabolic systems may be formally represented with the skew-symmetric operator J with an appropriate definition of the effort variables [START_REF] Vu | A structured control model for the thermo-magnetohydrodynamics of plasmas in tokamaks[END_REF].

The irrotational 2D SWE example

We will consider as a running example for this paper the irrotational 2D SWEs which describe the flow of an inviscid liquid where the horizontal components of the velocity field may be averaged on the water level and where the vertical velocity component may be omitted (low depth or shallowwater assumption). Besides, we will consider a 'non-rotating' flow. It is known that the corresponding 2D SWE then expresses the mass and momentum balance equations. Therefore, we will choose for the energy state variables the mass density (which is proportional to the water level h(z, t)) and the momentum density. For instance, using Cartesian coordinates, one would choose α q := h(z, t) (where h(z, t) denotes the water level) and α p := ρ u(z, t) v(z, t)

T where u(z, t) and v(z, t) denote the horizontal components of the fluid velocity, while ρ denotes the fluid mass density. Note that both α q and α p are defined in the two-dimensional (n = 2) horizontal spatial domain Ω of the flow. Using these energy state variables, one gets for the total (kinetic and potential) energy inside the domain Ω:

H(α q , α p ) := 1 2 Ω ρg (α q ) 2 + 1 ρ α q α p 2 d z. (3.10)
Therefore, the co-energy variables are defined as (3.11) which are the hydrodynamic pressure e q (an intensive variable in the 2D domain Ω) and the volume flow e p , respectively. For instance, using the same Cartesian coordinates as previously, one gets e q = ρgh + ρ 2 u 2 + v 2 and e p := h u(z, t) v(z, t) T . Using these co-energy variables, the Hamiltonian system of two canonically interacting conservation laws (3.4) reads .12) which are exactly the usual irrotational 2D SWE. For instance, (3.12) reads in Cartesian coordinates using the usual vector calculus notations The port-Hamiltonian model (3.4, 3.7) may then be implicitly defined as a linear subspace in the Bond space of effort and flow variables which embeds boundary effort and flow variables. In turn, this linear subspace may be geometrically defined as the linear subspace which is maximally isotropic with respect to some inner product associated to the natural power product-or power form in the Bond space-between effort and flow variables. Therefore, in the sequel, we aim at structure-preserving discretization which will preserve this power form in the approximation spaces. We will speak about symplectic discretization in the sense that this power form is preserved. In this section, we will define the Bond space, the power symplectic form, the associated inner product and the associated Stokes-Dirac structure which implicitly defines the port-Hamiltonian model (3.4, 3.7). Readers are referred to van der Schaft & Maschke ( 2002) for details about this representation.

e q := δ q H = ∂ α q H = ρgα q + 1 2ρ α p 2 , e p := δ p H = grad α p H = α q α p ρ ,
αq αp = 0 -div -grad 0 e q e p , ( 3
⎡ ⎣ ḣ ρ u v ⎤ ⎦ = 0 -div -grad 0 ⎡ ⎣ ρgh + ρ u 2 +v 2 2 h u v ⎤ ⎦ . ( 3
Let the Bond space of extended flow and effort variables be B := F × E with 3

F := L 2 (Ω) × L 2 (Ω; R n ) × H -1 2 (∂Ω), E := H 1 (Ω) × H div (Ω) × H 1 2 (∂Ω).
(3.14) 3 Let us recall that H 1 (Ω) is the space of scalar-valued functions with square integrable weak gradient, H div (Ω) is the space of vector-valued functions with square integrable weak divergence, and the duality bracket e ∂ , f

∂ ∂Ω := f ∂ , e ∂ H -1 2 (∂Ω), H 1 2 (∂Ω)
will be needed at the boundary, though for regular enough data it amounts to a scalar product in L 2 (∂Ω).

We may define on this Bond space the real power pairing or power form which maps any effort-flow vector (e, f ) ≡ (f q , f p , f ∂ ), (e q , e p , e ∂ ) ∈ B to (3.15) in such a way that every pair (e, f ) of extended effort and flow variables in the Bond space, satisfying the port-Hamiltonian equations (3.4) and (3.7), also satisfies the power balance equation e|f = 0. From the power pairing (3.15), we may define the following symmetric bilinear form:

e |f := e q , f q L 2 (Ω) + e p , f p L 2 (Ω;R n ) + e ∂ , f ∂ ∂Ω ,
•, • : B × B → R e 1 , f 1 , e 2 , f 2 := 1 2 e 1 f 2 + e 2 f 1 . (3.16)
With the help of this symmetric bilinear form (inner product on B ), we may define the Dirac structure associated to the power pairing (3.15) as the linear subspace D ⊂ B which is maximally isotropic, that is, such that D = D ⊥ where the orthogonality is defined with respect to the inner product •, • .

In particular, any (e, f ) ∈ D satisfies (e, f ) , (e, f ) = 0, hence the power balance e| f = 0. Dirac interconnection structure may be used to define implicitly the dynamics of pHs. In particular, the linear subspace D in the Bond space B := F × E , with F and E as in (3.15), which is defined by

D := ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (f q , f p , f ∂ ), (e q , e p , e ∂ ) ∈ B f q f p = 0 div grad 0 e q e p e ∂ f ∂ = γ 0 0 0 -γ ⊥ e q e p ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (3.17)
is a Dirac structure associated to the natural power pairing (3.15). This is proved by using the generalized Stokes theorem (van der [START_REF] Van Der Schaft | Hamiltonian formulation of distributed-parameter systems with boundary energy flow[END_REF]. Therefore, in this particular case, the interconnection structure is called a Stokes-Dirac structure. The dynamics (3.4), generated by the Hamiltonian function H(α q , α p ) (see Definition 3.1), with boundary energy flow and port boundary variables (3.7), may be implicitly defined by

(-αq , -αp , f ∂ ), (δ q H, δ p H, e ∂ ) ∈ D . (3.18)
In that sense, we will say that (3.18) defines a boundary pHs of two canonically interconnected conservation laws. In order to define a well-posed Cauchy problem, boundary conditions still need be chosen either for e ∂ or f ∂ .

Remark 3.1 The Stokes-Dirac structure (3.17), associated to the natural power bilinear form (3.15), may be used to represent hyperbolic systems, either linear (when the Hamiltonian density is quadratic) or non-linear (in the other cases). It may even be used to represent parabolic systems when the effort differential forms do not derive from the same Hamiltonian operator (Baaiu et al., 2009a). The choice of Dirichlet boundary conditions for e ∂ or f ∂ will lead to a well-posed system. There are, however, many other possible choices for admissible boundary conditions which lead to well-posed systems. For linear systems (quadratic Hamiltonian), these admissible boundary conditions have been parameterized

elegantly by Le [START_REF] Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF]. Many examples of physical systems have been represented using the port-Hamiltonian formulation and its implicit representation using Stokes-Dirac structures, including Maxwell field equations and Navier-Stokes flow problems (van der Schaft & Maschke, 2002), beam and membrane equations [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF], vibro-acoustic problems [START_REF] Trenchant | A port-Hamiltonian formulation of a 2D boundary controlled acoustic system[END_REF], shallow-water flow problems [START_REF] Cardoso-Ribeiro | Modeling of a fluid-structure coupled system using port-Hamiltonian formulation[END_REF][START_REF] Hamroun | Port-based modelling for open channel irrigation systems[END_REF][START_REF] Pasumarthy | Port-Hamiltonian formulation of shallow water equations with Coriolis force and topography[END_REF], advection-diffusion or adsorption problems (Baaiu et al., 2009a), Tokamak plasma MHD problems [START_REF] Vu | A structured control model for the thermo-magnetohydrodynamics of plasmas in tokamaks[END_REF], etc.

The partitioned weak form for pHs of two conservation laws

We will follow the approach presented in Section 2.2 for the 1D SWE and generalize it to the formulation (3.4) for pHs of two conservation laws. Let v q ∈ H 1 (Ω) and v p ∈ L 2 (Ω) denote any arbitrary test functions. We may obtain from the strong formulation (3.4) the following weak form:

v q , αq L 2 (Ω) = -v q , div e p L 2 (Ω) , v p , αp L 2 (Ω) = -v p , grade q L 2 (Ω) . (3.19)
Integrating by parts the first state equation only, and using the integration by parts formula, we get the following partitioned weak form:

v q , αq L 2 (Ω) = gradv q , e p L 2 (Ω) - ∂Ω v q ∂Ω e p ∂Ω • n ds, v p , αp L 2 (Ω) = -v p , grade q L 2 (Ω) .
(3.20)

Remark 3.2 With the particular choice v q = e q and v p = e p , one gets e q , αq L 2 (Ω) = grade q , e p L 2 (Ω) -∂Ω e q ∂Ω e p ∂Ω • n ds, e p , αp L 2 (Ω) =e p , grade q L 2 (Ω) .

(3.21)

Therefore, the power balance equation (2.6) reads

Ḣ = e q , αq L 2 (Ω) + e p , αp L 2 (Ω) = - ∂Ω e q ∂Ω e p ∂Ω • n ds, (3.22)
which shows that the power balance is also preserved in the weak formulation.

Remark 3.3 Instead of integrating the first equation by parts in (3.20), letting v q ∈ L 2 (Ω) and v p ∈ H div (Ω), we could integrate the second equation, which would also lead to another skew-symmetric structure (involving div instead of grad); the boundary inputs and outputs would not be the same either, more precisely their role would be switched.

In the case of the irrotational 2D SWE, written in Cartesian coordinates, the partitioned weak formulation (3.20) reads

Ω v q αq dx dy = Ω gradv q • e p dx dy - ∂Ω v q e p • n ds, Ω v p • αp dx dy = - Ω v p • grade q dx dy . (3.23)
In this latter equation, α q denotes the water level (proportional to the water surfacic density), while the bold vector notations α p := ρh(x, y, t) [u(x, y, t), v(x, y, t)] T , e p and v p denote the vectors of Cartesian coordinates for the momentum, water flow and test functions, respectively. With n as the outward unit normal to the boundary ∂Ω, the term -n • e p ∂Ω denotes the boundary port variable f ∂ (perpendicular water flow at the boundary) which will be chosen as the boundary input u ∂ ; thus, the collocated boundary output will be y ∂ := e ∂ = e q ∂Ω . In the following section, we discretize the partitioned weak form (3.23), and we show that the resulting system is a finite-dimensional pHs preserving the power balance of the original system at the discrete level.

PFEM for the 2D SWE

In this section, the partitioned weak-form representation for the SWEs, presented in Section 3.4, will be projected into finite-dimensional approximation spaces in such a way as to preserve the total mass and momentum conservation equations and the underlying Dirac structure of the original port-Hamiltonian model. This section is divided in two parts, firstly, the weak form is discretized in Section 4.1 and a finite-dimensional pHs is obtained. Secondly, we show how to obtain the discrete constitutive relationships in Section 4.2.

Structure-preserving finite element discretization

Let us approximate the energy variables α q (x, y, t) using the following basis with N q elements:

α q (x, y, t) ≈ α ap q (x, y, t) := N q i=1 φ i q (x, y)α i q (t) = φ q (x, y) T α q (t), (4.1)
where φ i q (x, y), i ∈ {1, . . . , N q } are the chosen approximation basis functions in H 1 (Ω), and α i q (t) are the approximation coordinates for α ap q (x, y, t). The test functions v q and the co-energy variables e q are also approximated using φ q (x, y).

Similarly, the vectorial variable α p (x, y, t) is approximated as .

α p (x, y, t) ≈ α p ap (x, y, t) := N p i=1 φ p i (x, y)α i p (t) = Φ p (x, y) T α p (t)
With this particular choice, we recover the case where the basis functions components of the vectorial variables are decoupled. We studied this case by [START_REF] Cardoso-Ribeiro | A structure-preserving partitioned finite element method for the 2D wave equation[END_REF].

Finally, the boundary input can be discretized using any one-dimensional set of basis functions, say

ψ(s) = [ψ i (s)]: u ∂ (s, t) ≈ u ap ∂ (s, t) := N ∂ i=1 ψ i (s)u i ∂ (t) = ψ(s) T u ∂ (t). (4.4)
Remark 4.2 In the sequel, in the implementation of our finite element method, we conveniently chose ψ(s) as φ q (x(s), y(s)) evaluated on the boundary. Other choices could be investigated.

From the substitution of the approximated variables, (4.1), (4.2) and (4.4), in the weak form (3.28), the following finite-dimensional equations are obtained:

v q T Ω φ q φ T q dx dy αq = v q T Ω ∂φ q ∂x ∂φ q ∂y Φ T p dx dy e p + -v q T ∂Ω φ q (x(s), y(s))Ψ T (s) ds u ∂ (t), v p T Ω Φ p Φ T p dx dy αp = -v p T ⎡ ⎣ Ω Φ p ⎡ ⎣ ∂φ T q ∂x ∂φ T q ∂y ⎤ ⎦ dx dy ⎤ ⎦ e q .
(4.5)

Since these equations should remain valid for any test functions coordinates v q and v p , one gets

M q αq (t) =De p (t) + Bu ∂ (t), M p αp (t) = -D T e q (t), (4.6)
where M q and M p are square mass matrices (of size N q × N q and N p × N p , respectively), defined as

M q := Ω φ q φ T q dx dy, M p := Ω Φ p Φ T p dx dy. (4.7)
The matrix D is of size N q × N p , defined as

D := Ω ∂φ q ∂x ∂φ q ∂y Φ T p dx dy, (4.8)
and B is an

N q × N ∂ matrix B := ∂Ω φ q (x(s), y(s))ψ T (s) ds. (4.9)
Defining y ∂ (t), the output conjugate to the input u ∂ (t) as y ∂ (t) := B T e q (t), (4.10)

the approximated system can be written using the following finite-dimensional Dirac structure representation:

M q 0 0 M p f q (t) f p (t) = 0 -D D T 0 e q (t) e p (t) + -B 0 u ∂ (t) y ∂ (t) = B T 0 e q (t) e p (t) , (4.11)
where f q (t) := -αq (t) and f p (t) := -αp (t) denote the vector coordinates for the flow approximations coordinates in the φ i q (x, y) and φ p i (x, y) approximation bases, that is, f ap q (x, y, t) = φ q (x, y) T f q (t) and

f ap p (x, y, t) = Φ p (x, y) T f p (t).
From the definition of the co-energy variables as the variational derivatives of the Hamiltonian with respect to α q (x, y, t) and α p (x, y, t), the time derivative of the continuous Hamiltonian is given by Ḣ = Ω αp (x, t) • e p (x, t) + αq (x, t)e q (x, t) dΩ.

(4.12)

Using the approximations for the energy and co-energy variables, this power balance can be approximated as

Ḣd =

Ω αp (t) T Φ p (x, y)Φ p (x, y) T e p (t) + αq (t) T φ q (x, y)φ q (x, y) T e q (t) dΩ, = αp (t) T M p e p (t) + αq (t) T M q e q (t).

(4.13)

The approximated equations (4.11), together with the power balance (4.13), provide a Dirac structure representation that is a projection of the Dirac-Stokes structure.

Therefore, in order to write the power balance as the total time derivative of the discrete Hamiltonian written as a function of the finite-dimensional vector coordinates for the energy variables, the following relationships must hold:

M q e q = ∂H d ∂α q , M p e p = ∂H d ∂α p , (4.14)
where the approximated Hamiltonian is defined as

H d (α q , α p ) := H[α q (x, t) = α T q (t)φ q (x), α p (x, t) = α T p (t)Φ p (x) .
(4.15)

From (4.11) and (4.13), the time derivative of the approximated Hamiltonian is given by

Ḣd (t) = αq (t) T M q e q (t) + αp (t) T M p e p (t), = e T p D T e q + u ∂ B T e q -e T q D e p (t) = y T ∂ u ∂ .
(4.16)

Note that u T ∂ y ∂ is the discrete analogue of the transferred power over the boundary (right side of power balance equation (2.19) in the continuous case). Furthermore, this power balance is exactly preserved in the finite-dimensional approximation spaces. From the definition of the B matrix (4.5), the definition of the approximated boundary input u ap ∂ (s, t) := ψ(s) T u ∂ (t) and approximated co-energy variable e ap p (x(s), y(s), t) := φ q (x(s), y(s)) T e q (t), we get

Ḣd = e T q ∂Ω φ q (x(s), y(s))ψ T (s) ds u ∂ , = ∂Ω e ap q (x(s), y(s), t)u ap ∂ (s, t) ds.
(4.17)

Remark 4.3 Note that using classical first-order finite-element one-dimensional discretization basis for the boundary input, the coordinates u ∂ (t) provide the values of the inflow (-n • e q ) at the boundary nodes. For instance, in the case of SWEs, these are the values of volumetric influx into the system. The conjugate output y ∂ is related with the curve integral of e q (x(s), y(s), t) along the elements. The coenergy variable e q (x(s), y(s), t) is the pressure; thus, the discretized outputs coordinates y ∂ (t) are related to the forces per unit length applied along the external boundary.

Note that we can also define an output that is given by the point-wise values of the co-energy variables e q (s, t) evaluated on the boundary. In this case, a convenient choice of basis function for the approximation would be the same as the input, i.e. (4.18) where the coordinates ŷ∂ (t) represents the values of pressure at the boundary nodes. The power balance through the boundary is computed as (4.19) where

y ap (s, t) = ψ T (s)ŷ ∂ (t),
Ḣd = ∂Ω y ap ∂ (s, t)u ap ∂ (s, t) ds = ŷT ∂ (t) ∂Ω ψ(s)ψ T (s) ds u ∂ = ŷT ∂ (t)M ψ u ∂ ,
M ψ = ∂Ω ψ(s)ψ T (s) ds (4.20) is a symmetric positive-definite N ∂ × N ∂ mass matrix.
Furthermore, since the power balances (4.16) and (4.19) must coincide, the following relationship between these two output definitions must hold:

M ψ ŷ∂ = y ∂ . (4.21)
Consequently, from (4.13) and (4.19), the following power product must hold (satisfying (4.11)):

⎡ ⎣ f p f p ŷ∂ ⎤ ⎦ ⎡ ⎣ e q e p u ∂ ⎤ ⎦ := f T p M p e p + f T q M q e q + u T ∂ M ψ ŷ∂ = 0. (4.22)
With the help of a symmetric bilinear form as (3.16), using the previous power product, a finitedimensional Dirac interconnection structure can be defined.

Obtaining the non-linear constitutive relationships: discretization of the Hamiltonian

In the previous section, a finite-dimensional Dirac structure was obtained for the 2D SWEs, relating the energy and co-energy variables as well as the boundary inputs and outputs. The next step is to obtain the constitutive relationships from the Hamiltonian.

The Hamiltonian of the 2D SWE (3.10) can be rewritten using the coordinate variables as

H[α q (x, y, t), α p (x, y, t)] := 1 2 Ω α q (x, y, t) α p 2 ρ + ρg(α q (x, y, t)) 2 dΩ . (4.23)
The energy variables are restricted to the approximation spaces for α q (x, y, t) and α p (x, y, t). From (4.15), the discretized Hamiltonian reads

H d (α q , α p ) := 1 2 Ω ⎛ ⎜ ⎝ φ T q (x, y)α q Φ T p (x, y)α p 2 ρ + ρg φ T q (x, y)α q 2 ⎞ ⎟ ⎠ dΩ . (4.24)
The reduced effort variables are obtained from the gradient of the discretized Hamiltonian:

∂H d ∂α q = ρgM q α q (t) + 1 2ρ Ω φ q (x, y)α T p (t)Φ p (x, y)Φ T p (x, y) dΩ α p (t) ∂H d ∂α p = α T q (t) 1 ρ Ω φ q (x, y)Φ p (z)Φ T p (z) dΩ α p (t).
(4.25)

Both constitutive equations exhibit non-linear terms. In order to compute them, the following procedure was used. The first equation can be written as

∂H ∂α q = ρgM q α q (t) + 1 2ρ Ω φ q (x, y)α T p (t)Φ p (x, y)Φ T p (x, y) dΩ α p (t) w(t) , (4.26)
where the components of w(t) can be computed as

w i (t) = α T p (t) Ω φ q,i (x, y)Φ p (x, y)Φ T p (x, y) dΩ α p (t). (4.27)
Note that Ω φ q,i (x, y)Φ p (x, y)Φ T p (x, y) dΩ, for {i = 1, . . . , N q } are N q matrices of dimension N p × N p , which can be computed once and remain constant. Similarly, the second constitutive relationship can be written as a function of constant matrices.

Numerical solution

The finite-dimensional Dirac structure (4.11) together with the constitutive relations (4.14) provides a finite-dimensional dynamical system in port-Hamiltonian form, detailed below:

M q 0 0 M p αq (t) αp (t) = 0 D -D T 0 e q (t) e p (t) + B 0 u ∂ (t) y ∂ (t) = B T 0 e q (t) e p (t) , (4.28)
where the co-energy variables are obtained from the constitutive relationships

M q e q (t) = ∂H d ∂α q α q (t), α p (t) , M p e p (t) = ∂H d ∂α p α q (t), α p (t) , (4.29)
and the initial conditions for the energy variables are: α q (t = 0) and α p (t = 0); some compatibility condition must be met relating the boundary values of the co-energy variable e p and the control input u ∂ at time t = 0. From a scientific computational point of view, three strategies are possible for solving the implicit equations (4.28) and (4.29) with sparse mass matrices. Let us give a short overview of these techniques.

1. Firstly, at each time-step, the co-energy variables can be obtained from (4.29) as a solution of the linear system. Then, (4.28) is solved using usual ODE numerical integration schemes in the time domain. This strategy has already been used on three-dimensional Maxwell's equations with 40.000 degrees of freedom (dof) [START_REF] Payen | Modelling and structure-preserving discretization of Maxwell's equations as port-Hamiltonian system[END_REF].

2. Secondly, both equations can be solved at the same time using implicit methods specific to DAEs (see, e.g. [START_REF] Egger | On structure-preserving model reduction for damped wave propagation in transport networks[END_REF][START_REF] Mehrmann | Structure-preserving discretization for port-Hamiltonian descriptor systems[END_REF], and the references therein).

3. Finally, when the number of dof is small, typically a few hundreds, it is also possible to find explicit pHs representations, involving the inversion of the mass matrices. For instance, defining new energy variables αp (t) := M p α p (t) and αq (t) := M q α q (t), it is possible to rewrite (4.28) and (4.29) as

αq αp = 0 D -D T 0 ∂H d ∂ αq ∂H d ∂ αp + B 0 u ∂ , y ∂ =B T ẽq , (4.30)
where the approximated Hamiltonian is given by

H d ( αq , αp ) := H α ap q (x, y, t) = M -1 q αq (t)φ q (x, y), α ap p (x, y, t) = M -1 p αp (t)Φ p (x, y) .
(4.31)

The code corresponding to this latter implementation on the two-dimensional shallow water can be found in Cardoso-Ribeiro et al. (2020a) with 1024 dof.

Numerical experiments

In this section, we present numerical experiments to test the PFEM. Firstly, results for the 1D SWE are presented in Section 5.1. Then, the two-dimensional case is presented in Section 5.2. The source codes for the numerical results presented in this section are available online at Cardoso-Ribeiro et al. (2020b).

One-dimensional SWEs

The example presented in Section 2 was implemented using finite elements with polynomial basis functions.

Firstly, a spectral convergence analysis of the numerical method was done. The eigenvalues obtained from the linearized 4 numerical model were compared with the exact eigenfrequencies of the linear wave equation with constant coefficients. The inputs u ∂ (t) of (2.26) were considered to be zero.

Recall that the variables v q (z, t), e q (z, t) and q(z, t) must be discretized in z with polynomials of order at least one (since they are derived once on (2.16)). Figure 1 shows the relative error of the first modal frequency for four different choices of polynomial approximations. P 1 P 0 stands for first-order polynomial for the variables related to q, and order zero for the variables related to p. P 1 P 1 uses firstorder polynomial for both all variables. P 1 P 2 uses first-order polynomial for the q variables and p variables. Finally, P 3 P 3 uses third-order polynomial for all variables.

Secondly, time-domain simulations for the SWE 4 Hereafter, the linearized partial differential equations are used for validation purposes. They are obtained assuming small amplitudes of variations of height h(t, z) and velocity ṽ(t, z) (such that h(t, z) = h + h(t, z) and v(t, z) = v + ṽ(t, z), with v = 0). Consequently, (2.1) can be approximated using only the first-order terms as

∂ ∂t h = - h ∂ ∂z ũ, ∂ ∂t ũ = -g ∂ ∂z h,
which is a linear wave equation with constant coefficients. The previous equation can be written as a pHs, using p = ρ ṽand q = b h, with exactly the same structure as (2.4). The only difference is that in the linear case, the Hamiltonian (2.3) becomes quadratic: were performed. The following initial conditions were considered:

H l = 1 2 [0,L] b h p2 ρ + ρg b q2 dz.
q(z, t = 0) = 1, 0 z L, p(z, t = 0) = 0, (5.1)
together with the following boundary conditions:

e p (0, t) = e p (L, t) = A t 1 + t sin(ωt), (5.2)
such that these conditions represent a harmonic influx through both boundaries. In order to avoid spurious oscillations, the initial condition of the control input must agree with the initial conditions within the domain (avoiding discontinuities). Furthermore, the term t/(1 + t) is included to smooth the control input at initial condition.

The simulations were performed using 20 elements of type P 1 P 1 . Trapezoidal method was used as time integration scheme.

Two different amplitudes A were used in the simulations. Snapshots of the simulations are presented in Figs 2 and3. The first figure shows snapshots for a small amplitude value for the input inflow. In the second figure, the amplitude is multiplied by 100. Non-linear phenomenon is observed in this case. Finally, Fig. 4 presents the evolution of the height in time, considering the non-linear simulation with large input amplitude.

Two-dimensional SWEs

The example presented in Section 4 was implemented using quadrilateral finite elements with polynomial basis functions for a square domain.

As we did for the 1D SWE in the previous subsection, firstly, a spectral convergence analysis of the numerical method was done.

For the 2D SWE, recall that the variables with index q (α q , e q , v q ) must be discretized with polynomials of order at least one (since they are derived once on (3.23)). Figure 5 shows the relative error of the first modal frequency for three different choices of polynomial approximations. P i P j P j stands for 'i-th' order polynomial for the variables with index q, and 'j-th' order for the two components of the vector variables of index p (α p , e p and v p ). The sparsity of the of the PFEM matrices M q , M p and D for elements of type P 1 P 0 P 0 are presented in Fig. 6. Time-domain simulations were performed using the discretized system under boundary-port excitation. The following initial conditions were considered:

α q (x, y, t = 0) = 1, [x, y] ∈ Ω, α p (x, y, t = 0) = 0,
(5.3) together with the following boundary conditions: (5.4) and u ∂ (x, y, t) = 0, [x, y] ∈ ∂Ω, t > 1s. The boundary is split in four sides: ∂Ω = ∂Ω up ∪ ∂Ω left ∪ ∂Ω down ∪ ∂Ω right . These conditions impose a harmonic inflow on one side of the boundary and the opposite condition on the other side. The simulations were performed using a 10 × 10 grid of P 1 P 0 elements. Again, trapezoidal method was used for time integration. Figure 11 shows how the approximated Hamiltonian (4.24) and V(t) = Ω α ap 1 (x, y, t) dΩ, the total volume of fluid, change with time. As expected, the Hamiltonian only changes during the first second of simulation, while the system is excited through the boundary ports. After that, since Ḣd = u T ∂ y ∂ = 0, the Hamiltonian is constant. Finally, the total volume is kept nearly constant and the changes are only due to numerical precision (of order 10 -14 ).

u ∂ (x, y, t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ A t 1 + t sin(π t) , [x, y] ∈ ∂Ω up , -A t 1 + t sin(π t) , [x, y] ∈ ∂Ω left , 0 , [x, y] ∈ ∂Ω down ∪ ∂Ω right . t 1s,

Extensions

In this final section, the extension of PFEM to more general configurations is addressed, in order to illustrate the flexibility of the proposed numerical method: in Section 6.1 the specific expression in polar coordinates is investigated, in Section 6.2 the case of heterogeneous medium with variable coefficients is presented and in Section 6.3 the method is applied to second-order differential operators, such as those involved in the port-Hamiltonian formulation of the Euler-Bernoulli beam in one-dimensional.

Polar coordinates

The goal of this subsection is to prove the applicability of PFEM in two-dimensional, when the chosen coordinate system is not Cartesian: polar coordinates are presented. In order to avoid unnecessary technicalities, the geometry of a disc has been chosen to illustrate this extension. Further control applications on this test case can be found in [START_REF] Cardoso-Ribeiro | Port-Hamiltonian modeling, discretization and feedback control of a circular water tank[END_REF].

The 2D SWEs in a disc, as a pHs in polar coordinates. Let us consider the disc Ω = D R of radius R > 0 with boundary ∂Ω = C R , the circle of radius R. Polar coordinates r and θ will be used. In vector calculus, the energy variables are represented by the scalar function α q = h, and the vector function α p := ρ [u r (t, r, θ), u θ (t, r, θ)] T . The Hamiltonian reads

H = 1 2 D R [ρg h 2 + ρh ((u r ) 2 + (u θ )
2 )] r dr dθ, (6.1) The effort or co-energy variables can be computed as e q := δ q H = ρg α q + 1 2ρ |α p | 2 and e p := δ p H = 1 ρ α q α p = h [u r (t, r, θ), u θ (t, r, θ)] T . With these notations and definitions, we get the same system as (3.13) for the strong form of the pHs, namely (6.3) with boundary control u ∂ (θ , t) := -e p •n = -e r p (R, θ , t) and collocated boundary observation y ∂ (θ , t) := e q (R, θ , t) at the boundary ∂Ω = C R . Let us conclude with the energy balance for this system,

= D R [ 1 2 ρg α 2 q + 1 2ρ α q |α p | 2 ] r dr dθ . (6.2)
⎡ ⎣ ḣ ρ ur uθ ⎤ ⎦ = 0 -div -grad 0 ⎡ ⎣ ρ(gh + (u r ) 2 +(u θ ) 2 2 ) h u r u θ ⎤ ⎦ ,
d dt 1 2 D R [ρg h 2 + ρh ((u r ) 2 + (u θ ) 2 )] r dr dθ = C R u ∂ (θ , t) y ∂ (θ , t) R dθ . (6.4)
The PFEM directly applies to the pHs in polar coordinates. Now, let us first rewrite the weak form (3.23) with test functions v q and v p .

D R v q αq r dr dθ = D R gradv q • e p r dr dθ - C R v q n • e p R dθ , D R v p • αp r dr dθ = - D R
v p • grade q r dr dθ .

(6.5) Let us approximate the scalar energy variables α q (r, θ , t) using the following basis with N q elements:

α q (r, θ , t) ≈ α ap q (r, θ , t) := N q i=1 φ i q (r, θ)α i q (t) = φ q (r, θ) T α q (t). (6.6)
The variables e q and v q are also approximated using φ q (r, θ).

Similarly, the vectorial energy variable α p is approximated as

α p (r, θ , t) ≈ α p ap (r, θ , t) := N p k=1 φ p k (r, θ)α k p (t) = Φ p (r, θ) T α p (t), (6.7)
where

φ p k (r, θ) = φ r,k p (r, θ) φ θ,k p (r, θ)
represents a two-dimensional vectorial basis function and, consequently, Φ p (r, θ) is an N p × 2 matrix. Remark 4.1 does apply here also. Furthermore, e p and v p are also approximated using Φ p (r, θ). Finally, the boundary input, localized on the circle of radius r = R can be discretized using any onedimensional set of basis functions, say ψ = [ψ m ], provided 2π -periodicity is ensured (trigonometric polynomials are a fair trial approximation basis, see, e.g. (Boyd, 2001, Chapter 18)):

u ∂ (θ , t) ≈ u ap ∂ (θ , t) := N ∂ m=1 ψ m (θ )u m ∂ (t) = ψ(θ ) T u ∂ (t). (6.8)
Introducing the notation ∂ r φ q := [∂ r φ i q ] and ∂ θ φ q := [∂ θ φ i q ] for the matrices of partial derivatives of the functions φ i q , we define matrix (6.9) where the apparent singularity at r = 0 has been removed. Then, with classical mass matrices M q := D R φ q φ T q r dr dθ, M p := D R Φ p Φ T p r dr dθ , together with the control matrix B := (6.10) where M q and M p are square matrices (of size N q × N q and N p × N p , respectively). D is an N q × N p matrix and B is an N q × N ∂ matrix.

D := D R ∂ r φ q r -1 ∂ θ φ q Φ T p r dr dθ = D R r ∂ r φ q ∂ θ φ q Φ T p dr dθ ,
C R φ q (R, θ) ψ T (θ ) R dθ , the finite-dimensional equations become M q αq =D e p + B u ∂ (t), M p αp = -D T e q ,
Defining y ∂ (t), the output conjugate to the input u ∂ (t) as y ∂ (t) := M ψ -1 B T e q (t), (6.11) with boundary mass matrix M ψ := C R ψψ T R dθ , the approximated system can be written using the finite-dimensional Dirac structure representation given by (4.11), and as found in Remark 4.3, the global energy balance reads (6.12) which mimicks that at the continuous level, namely (6.4).

d dt 1 2 (α T q M q α q + α T p M p α p ) = y T ∂ M ψ u ∂ ,

Heterogeneous case with variable coefficients

The goal of this subsection is to prove the applicability of PFEM when the coefficients are space varying.

In order to avoid unnecessary technicalities, the choice has been made to tackle the one-dimensional model, first derived in Section 2 as introductory example. Another fully worked out example can be found in Serhani et al. (2019b) on the anisotropic heterogeneous wave equation in two-dimensional.

The variable-coefficient physical model as a pHs. Let us consider the SWE in a water channel with a space-varying cross section, i.e. with z → b(z) the width of the channel, it is easy to understand that the energy given by (2.2) remains unchanged, with function b(z) instead of coefficient b, but then the balance equations (2.1) must be modified as follows:

∂ ∂t (b h) = - ∂ ∂z (bhu) , ∂ ∂t u = - ∂ ∂z u 2 2 + gh . (6.13)
With the appropriate choice of energy variables q := b h, and p := ρ b u, the Hamiltonian (2.2) now reads

H (q(z, t), p(z, t)) = 1 2 [0,L] qp 2 ρ b 2 + ρg b q 2 dz. (6.14)
The co-energy variables are found to be e q := δ q H = ρ u 2 2 + gh , and e p := δ p H = h u; thus, system (6.15) becomes in compact form

q(z, t) = - ∂ ∂z b(z) e p (z, t) , ṗ(z, t) = -b(z) ∂ ∂z e q (z, t), (6.15)
to be compared with (2.4) in the uniform case. Hence, the new interconnection operator J b reads The PFEM directly applies to the pHs with variable coefficients. Here, the same procedure as in Section 2.2 is being followed. We begin with a weak formulation of (6.15), then two complementary choices can be made.

J b := 0 -∂ z [b(z) .] -b(z) ∂ z 0 . Since [0,L] ϕ ∂ z (b(z) ψ) dz = -[0,L] b(z)(∂ z ϕ) ψ
If we choose to integrate by parts the mass balance equation only, i.e. the first line of the obtained weak form, we get exactly the same finite-dimensional pHs as (2.36), but with D a new N q × N p matrix, defined by

D := [0,L] b(z) dφ q dz (z) φ T p (z)dz, (6.16)
and a new N q × 2 control matrix B := b(0) φ q (0) -b(L) φ q (L) . The boundary control remains u ∂ (t) := e p (0, t) e p (L, t) T , and the new collocated boundary observation reads y ∂ (t) := B T e q (t) = b(0) e q (0, t) -b(L) e q (L, t) T .

If instead, we choose to integrate by parts the momentum balance equation only, i.e. the second line of the obtained weak form, we get the following finite-dimensional pHs: M q q(t) = De p (t) , M p ṗ(t) = -D T e q (t) + B e q (0, t) e q (L, t) , (6.17) but with D another N q × N p matrix, defined by Finally, note that in both the above cases, the following power balance is met:

D := - [0,L] φ q (z) d dz [b(z) φ T p (z)]dz, ( 6
Ḣd (t) := e T q (t)M q q(t) + e T p (t)M p ṗ(t) = y T ∂ (t) u ∂ (t).

Higher-order systems

In the previous sections, the PFEM was applied to first-order (one-and two-dimensional) formally skewsymmetric differential operators. Indeed, the method seems to be much more general and can be applied similarly to higher-order equations.

The Euler-Bernoulli beam equation can be written as a pHs of second order. The equations are given by (see, e.g. Cardoso-Ribeiro et al. ( 2016)):

ẋ1 (z, t) = - ∂ 2 ∂z 2 e 2 (z, t), ẋ2 (z, t) = ∂ 2 ∂z 2 e 1 (z, t), (6.19)
where e 1 and e 2 are obtained from the variational derivative of the Hamiltonian:

H = 1 2 L 0 x 1 2 + x 2 2 dz. (6.20)
From the definition of the variational derivatives, the time derivative of the Hamiltonian is computed as

Ḣ = L z=0 e 1 ẋ1 + e 2 ẋ2 dz, = L z=0 -e 1 ∂ 2 z 2 e 2 + e 2 ∂ 2 z 2 e 1 dz, = L z=0 ∂ z -e 1 ∂ z (e 2 ) + ∂ z (e 1 ) e 2 dz, = -e 1 (L, t) ∂ z (e 2 )(L, t) + ∂ z (e 1 )(L, t)e 2 (L, t)
+ e 1 (0, t) ∂ z (e 2 )(0, t) -∂ z (e 1 )(0, t)e 2 (0, t). (6.21)

Note that Ḣ depends only on the boundary values of e 1 (vertical speed), e 2 (moment), ∂ z e 1 (rotation speed) and ∂ z e 2 (force). This motivates the definition of the boundary ports, which allows writing the infinite-dimensional equations as pHs. From (6.21), one possible definition for the boundary ports is as follows: Integrating the first equation by parts twice, we get the following partitioned weak form: (6.27) and B is an N 1 × 4 matrix:

y ∂ := ⎡ ⎢ ⎢ ⎣ f 1∂ f 2∂ f 3∂ f 4∂ ⎤ ⎥ ⎥ ⎦ := ⎡ ⎢ ⎢ ⎣ ∂ z e 1 (L
L 0 v 1 (z)ẋ 1 (z, t)dz = - L 0 ∂ 2 ∂z 2 v 1 (z) e 2 (z, t)dz + ∂ z v 1 (L) -∂ z v 1 (0) -v 1 (L) v 1 (0) ⎡ ⎢ ⎢ ⎣ e 2 (L
D := L z=0 ∂ 2 φ 1 ∂z 2 (z) φ 2 (z) T dz ,
B := ∂φ 1 ∂z (L) - ∂φ 1
∂z (0) -φ 1 (L) φ 1 (0) . (6.28)

Finally, the conjugate-output can also be written in terms of the previous B matrix: Defining the flow variables as f 1 (t) := -ẋ 1 (t) and f 1 (t) := -ẋ 2 (t), we find the following finitedimensional Dirac structure representation: (6.30)

y ∂ = ⎡ ⎢ ⎢ ⎣
M 1 0 0 M 2 f 1 (t) f 2 (t) = 0 -D D T 0 e q (t
Following the same procedure presented in the previous sections for the 1D and 2D SWEs, from the discretization of the Hamiltonian using the energy variables approximation spaces, one gets the underlying port-Hamiltonian dynamics for the approximated Euler-Bernoulli beam equations.

The analogue of the Euler-Bernoulli beam in one-dimensional is the Kirchhoff plate in twodimensional, one can refer to Brugnoli et al. (2019b) for the modelling as a pHs using tensor calculus, and the application of PFEM to it, with various boundary controls; note that the analogue of the Timoshenko beam in one-dimensional is the Mindlin plate in two-dimensional, and PFEM can also be applied to this model-see Brugnoli et al. (2019a).

Conclusion and open questions

The partitioned finite element method (PFEM) provides a full-rank structure-preserving representation of pHs in two-and three-dimensional: a general setting has been proposed here, written in the language of vector calculus for common PDE applications. This method can be easily implemented thanks to ready to use FEM software to compute the matrices of the representation, which are all sparse. It applies to complex geometries, works in any coordinate systems and allows for space-varying coefficients; moreover, higher-order differential operators can also be tackled. Although PFEM has already been successfully applied to linear PDEs with quadratic Hamiltonian functionals, e.g. vibrating membranes and plates, here the methodology carries over to a non-linear PDE with non-quadratic and non-separable Hamiltonian functional, the irrotational SWE in two-dimensional.

The ongoing work include the mixed boundary control (possibly leading to differential algebraic problems as pHDAEs), see, e.g. Brugnoli et al. (2020a); the three-dimensional application of PFEM to Maxwell's equation, see, e.g. [START_REF] Payen | Modelling and structure-preserving discretization of Maxwell's equations as port-Hamiltonian system[END_REF]; the application of PFEM to the heat equation, see, e.g. Serhani et al. (2019a); the mathematical convergence analysis (choice of the finite element bases and theoretical rate of convergence), see, e.g. [START_REF] Haine | Numerical analysis of a structure-preserving space-discretization for an anisotropic and heterogeneous boundary controlled n-dimensional wave equation as port-Hamiltonian system[END_REF]; and the development of a unified computational framework, see, e.g. Brugnoli et al. (2020b).

Future work include the description of PFEM in a more general setting making use of exterior calculus on differential forms, following, e.g. [START_REF] Flanders | Differential Forms with Applications to the Physical Science[END_REF] or [START_REF] Frankel | The Geometry of Physics: An Introduction[END_REF] for their use in physical systems modelling, and, e.g. [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF] and [START_REF] Arnold | Spaces of finite element differential forms[END_REF] for the numerical approximations using finite element spaces of differential forms. Another topic of interest is the study of some one-dimensional test cases to analyse the case of appearance of a shock in finite time, even with regular initial data: it is known that the energy might not be preserved but transferred to the thermal domain thanks to an entropy production. Some worked-out two-dimensional test cases will be studied on coupled systems, e.g. fluid-structure interaction, or thermal-structure coupling. Lastly, structurepreserving model reduction techniques will be tested on the high-fidelity finite-dimensional systems obtained by PFEM; see, e.g. [START_REF] Egger | On structure-preserving model reduction for damped wave propagation in transport networks[END_REF] for pHs or [START_REF] Hauschild | Model reduction techniques for linear constant coefficient port-Hamiltonian differential-algebraic systems[END_REF] for pHDAEs. The reduced order system will then be most useful to apply dedicated control laws for pHs, like IDA-PBC, which do take advantage of the specific structure of these dynamical systems with collocated inputs and outputs; see, e.g. [START_REF] Ortega | Control by interconnection and standard passivity-based control of port-Hamiltonian systems[END_REF].

  The geometric structure of pHs In the previous section, we proposed a port-Hamiltonian formulation for open systems of two canonically interacting conservation laws (distributed state space equations (3.4) and the boundary equations (3.7)). Since boundary energy flows are considered, boundary port variables are needed to derive the power balance equation (3.8).

Fig. 1 .

 1 Fig. 1. Convergence of the first natural frequency for different polynomial interpolation of the basis functions.

Fig. 2 .

 2 Fig. 2. Four snapshots of a time-domain simulation representing the fluid height as a function of horizontal position z. This result uses a harmonic boundary excitation with very small amplitude. The non-linear and the linearized equations exhibit very similar results.

Fig. 3 .

 3 Fig. 3. Four of a time-domain simulation representing the fluid height as a function of horizontal position z. This result uses a harmonic boundary excitation with larger amplitudes, such that non-linear phenomenon is now observed.

Fig. 4 .

 4 Fig. 4. Time-domain simulation representing the fluid height as a function of horizontal position z and time t. This result uses a harmonic boundary excitation with large amplitude, as in the previous figure.

Fig. 5 .

 5 Fig. 5. Convergence of the first natural frequency of the two-dimensional linear SWEs.

Fig. 6 .

 6 Fig.6. Sparsity of the matrices M q , M p and D. The black dots represents the non-zero values. P1 elements were used for the q variables and P0 elements were used for the p variables.

Fig. 7 .

 7 Fig. 7. Snapshots of simulation for a harmonic inflow excitation at two of the boundaries of the domain. The variable α q (fluid height) is shown. Here, small inputs are considered, such that the non-linear and linearized time responses are almost equivalent.

Fig. 8 .

 8 Fig. 8. Snapshots of simulation for a harmonic inflow excitation at two of the boundaries of the domain. The variable α q (fluid height) is shown along a cross-section in the middle of the domain. Here, small inputs are considered, such that the non-linear and linearized time responses are almost equivalent.

Fig. 9 .

 9 Fig. 9. Snapshots of simulation for a harmonic inflow excitation at two of the boundaries of the domain. The variable α q (fluid height) is shown. Differences between the non-linear and linearized time responses are now observed.

Fig. 10 .

 10 Fig. 10. Snapshots of simulation for a harmonic inflow excitation at two of the boundaries of the domain. The variable α q (fluid height) is shown along a cross-section in the middle of the domain. Differences between the non-linear and linearized time responses are observed.

Fig. 11 .

 11 Fig. 11. (a) Hamiltonian as a function of time for the simulation for a harmonic inflow excitation at two of the boundaries of the domain. Hamiltonian only changes during the first second, while an external excitation is applied. (b) Total volume of fluid as a function of time for the simulation for a harmonic inflow excitation at two of the boundaries of the domain. The total volume is constant along all the simulation.

  dz for smooth scalar functions ϕ and ψ, a straightforward computation shows that (J b u , v) = -(u , J b v) for smooth vector-valued functions u and v vanishing at the ends of the interval, and with the standard scalar product in L 2 × L 2 . Hence, the unbounded matrix-valued differential operator J b proves skew-symmetric in L 2 × L 2 .

  .18) and a new N p × 2 control matrix B := b(0) φ p (0) -b(L) φ p (L) . The boundary control is now defined by u ∂ (t) := e q (0, t) e q (L, t) T , and the new collocated boundary observation reads y ∂ (t) := B T e p (t) = b(0) e p (0, t) b(L) e p (L, t) T .

  Remark 2.1 These boundary input and output variables are the classical hydrodynamic pressure and water flow which are indeed power-conjugate variables for lumped model hydraulic systems. Besides, this canonical choice of input-output variables satisfies the matrix well-posedness condition proposed by Le[START_REF] Gorrec | Dirac structures and boundary control systems associated with skew-symmetric differential operators[END_REF] for linear pHs. Indeed, with this choice, equation (2.6) precisely states that

	u ∂ :=	e p (0, t) e p (L, t)	,	(2.7)
	while the power-conjugate boundary output variables are defined as
	y ∂ :=	e q (0, t) -e q (L, t)	.	(2.8)

systems (2.5) and (2.7) are impedance passive, even impedance energy preserving. This property, which is sufficient in the linear case, should be useful to address the well-posedness issues in the non-linear case as well.

  in the Φ p (x, y) approximation space. The test functions v p and the co-energy variable e p are also approximated using Φ p (x, y).

	where φ p	i (x, y) =	φ x,i p (x, y) φ y,i p (x, y)	represents a two-dimensional vectorial basis function and, conse-
	quently, φ Remark 4.1 Note that a particular choice for Φ p is
				Φ p (x, y) =	φ x p (x, y) 0	φ	0 p (x, y) y	(4.3)
	such that we can decompose the variables with index p in their Cartesian components as α p =	α x p y α p
						,	(4.2)

p (x, y) is an N p ×2 matrix. The variables α i p (t) are the approximation coordinates for α ap p

(x, y, t)

form representation of Euler-Bernoulli beam equation.

  Let us use arbitrary test functions v 1 (z) and v 2 (z) and develop a weak form of (6.19):

		⎤ ⎥ ⎥ ⎦ , u ∂ =	⎡ ⎢ ⎢ ⎣	e 1∂ e 2∂ e 3∂	⎤ ⎥ ⎥ ⎦ =	⎡ ⎢ ⎢ ⎣	e 2 (L, t) e 2 (0, t) ∂ z e 2 (L, t)	⎤ ⎥ ⎥ ⎦ .	(6.22)
						e 4∂			∂ z e 2 (0, t)
	The final power balance ( Ḣ) can thus be written as
		Ḣ = y T ∂ u ∂ .					(6.23)
	Weak-L 0 L 0	v 1 (z)ẋ 1 (z, t)dz = -v 2 (z)ẋ 2 (z, t)dz = 0	L	L v 2 (z) 0 v 1 (z) ∂z 2 e 1 (z, t) dz, ∂ 2 ∂z 2 e 2 (z, t)dz, ∂ 2	(6.24)

, t) -∂ z e 1 (0, t) -e 1 (L, t) e 1 (0, t)

  Similarly to the development in Section 2.2, we chose finite-dimensional basis functions (2.20) as φ 1 (z) and φ 2 (z), for the variables with index 1 and 2, respectively. From the substitution of the approximation functions in the weak form (6.25), we find M 1 ẋ1 (t) = -De 2 (t) + B M 1 and M 2 are square mass matrices (of size N 1 × N 1 and N 2 × N 2 , respectively), equivalent to (2.23). The matrix D is of size N 1 × N 2 :

	⎡ ⎢ ⎢ ⎣	e 2 (L, t) e 2 (0, t) ∂ z e 2 (L, t) ∂ z e 2 (0, t)	⎤ ⎥ ⎥ ⎦ ,	(6.26)
	M 2 ẋ2 (t) = D T e 1 (t)
			⎤
			⎥ ⎥ ⎦ .

, t) e 2 (0, t) ∂ z e 2 (L, t) ∂ z e 2 (0, t) L 0 v 2 (z)ẋ 2 (z, t)dz = L 0 v 2 (z) ∂ 2

∂z 2 e 1 (z, t)dz, (6.25)

Finite-dimensional pHs.

where

In this definition, Riesz identification theorem has been used to identify continuous linear functionals with functions, relying on the use of the canonical scalar product, but in general, care must be taken of the chosen scalar product in the pivot Hilbert space.

Funding

ANR-16-CE92-0028 (entitled Interconnected Infinite-Dimensional Systems for Heterogeneous Media, financed by the Agence Nationale de la Recherche and the Deutsche Forschungsgemeinschaft). Further information is available at: https://websites.isae-supaero.fr/infidhem/the-project.


Laura Cyr

Bernard Laura

Pedinielli Christophe Cutarella

Vincent Bréjard

come

Context

Craving is defined as a "strong desire or urge to consume a substance" (APA, 2013). It is a central symptom in substance use disorders. Researches on substance cues reactivity in addiction showed the influence of situational and affective cues-elicited craving (e.g., Preston et al., 2009) and highlighted strong links between negative affects (NA) and craving reactivity (e.g., Cooney et al., 1997). NA and craving are both know to be associated with substance use and to interfere with intervention outcomes. However, it is still not known how these variables relate to each other and if some others affective variables, for instance positive affects (PA) and alexithymia, may play an important role.

AFFECTIVITY, ALEXITHYMIA AND CRAVING IN SUBSTANCE USE DISORDERS :

A SYSTEMATIC REVIEW L. Cyr 1 2 , L. Bernard 1 2 , C. Cutarella 2 , J.-L. Pedinielli 1 , V. Bréjard 1

Objectives

This review aims to synthesize independent and dependent affective variables ( dispositional trait or transient state) with craving and identify gaps in the literature.

Method

Data base used: PsycINFO, PubMed/Medline, Web of Science (+ grey literature).

Research algorithm:

Alexithymia OR emotion regulation OR emotion dysregulation OR negative affect OR positive affect AND craving OR drug craving OR alcohol craving OR urges.

Negative affectivity (trait and state) and craving

• Risk factor : mood/anxiety disorders high associated with craving. • Causal link (n=7).

• Moderating effect (n=6) of personality variables (self-consciousness, anxiety sensitivity, emotional. intelligence) which act like protective or risk factors.

Positive affectivity (trait and state) and craving

• Protective or risk factor.

• Mediating effect of PA between personality variables (impulsivity, hedonic capacities) and craving. • Majority of studies on tobacco (n=3).

Alexithymia (trait) , affectivity and craving

• One study on the link between alexithymia, affectivity and craving. • Positive correlation between alexithymia and craving (n=3). • Alexithymia moderated the association between depression on craving in dfferent ways according to gender (protective and risk factor).

Study selection criterions

• Period: 1990 to 2019.

• Type of studies: 12 longitudinal (until 6 months), 9 cross sectional.

• Type of assessments: All quantitative studies included self assessment (with satisfactory validity).

• Size and age of samples: 30 to 362 subjects, 18 to 66 years old.

• Type of substances: Alcohol (n=9), Tobacco (n=9), Opioids (n=5), Cannabis (n=2), Methamphetamine (n=1), Benzodiazepine (n=1).

Conclusions

• This systematic review supports the role of affective individual differences in reactivity to craving among substance dependence. • Overall, a strong association between affective disturbances (depression, anxiety) and craving was found including mediating or moderating effects of personality and level of positive affectivity despite of little data on negative affectivity like a trait. • There is a lack of consideration on the thymic-dependent role of alexithymia in craving whereas it's strong association with depression has been demonstrated for substance use disorders (Thorberg et al., 2009). • Future researches should included longitudinal design to determine directional pathways from affectivity to craving and clarify its moderating and/or mediating effects according to substances types.

1

  Aix-Marseille Université, LPCPP, Aix-en-Provence, France, 2 Clinique Saint-Barnabé, Ramsay Santé group, Marseille, France Results of the 21 studies selected after the selection process of psychiatry -EPA 2020 (Madrid) Sources • APA (2013). Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition. Washington, DC : American Psychiatric Association. • Cooney, N.L., Litt, M.D, Morse, P.A, Bauer, L.O, & Gaupp L. (1007). Alcohol cue reactivity, negative-mood reactivity, and relapse in treated alcoholic men. Journal of abnormal Psychology, 106, 243-250. • Preston, K. L., Vahabzadeh, M., Schmittner, J., Lin, J. L., Gorelick, D. A. & Epstein, D. H. (2009). Cocaine craving and use during daily life. Psychopharmacology, 207(2), 291-301. doi: 10.1007/ s00213-009-1655-8 • Thorberg, F.A., Young, R.M., & Sullivan, K.A. (2009). Alexithymia and alcohol use disorders: a critical review. Addictive Behaviors, 34(3), 237-245. doi: 10.1016/j.addbeh.2008.10.016